/*Copyright (C) 1992, 1995 by Thomas Glen Smith.  All Rights Reserved.*/
/* matinv APL2 V1.0.0 **************************************************
* Produces the matrix inverse.  Implements APL monadic DOMINO.         *
* An argument of rank greater than two is rejected.  A scalar is       *
* treated as a 1-by-1 matrix, and a vector as a one-column matrix.     *
***********************************************************************/
#define INCLUDES APLCB
#include "includes.h"
Aplcb matinv(rite)
Aplcb rite;
{
	Errinit; Errstop; Getcb; Invert; Mdivarg; Real; Regress;
	extern int aplerr;
	int cols,*dimptr,i,rows;
	double *dataptr,f,*inv;
	Aplcb out=NULL;

	if (errinit())
		return(errstop(0,NULL,rite,out));
	if (ERR == mdivarg(rite,&rows,&cols))
		return(errstop(0,NULL,rite,out));
	if (rows < cols)
		return(errstop(6,NULL,rite,out));
	if (!(rite->aplflags & APLNUMB)) {
		rite = real(rite);
		if (aplerr) return(NULL);
	}
	if (rows == cols) { /* square matrix? */
		dataptr = invert(rite->aplptr.apldata,rows);
		if (NULL == dataptr)
			return(errstop(8,NULL,rite,out)); /* not invertible */
	}
	else { /* non-square matrix */
		dataptr = regress(NULL,rite,rows,rows,rows,cols);
		if (aplerr) return(NULL);
	}
	out=getcb(dataptr,rite->aplcount,APLTEMP+APLNUMB,rite->aplrank,NULL);
	if (rite->aplrank > 1) { /* build new dimensions? */
		dimptr=out->apldim;
		*dimptr++=*(rite->apldim+1);
		*dimptr=*(rite->apldim);
	}
	return(errstop(0,NULL,rite,out));
}