/*Copyright (C) 1992, 1996 by Thomas Glen Smith.  All Rights Reserved.*/
/* invert APL2 V1.0.0 *************************************************
* Called by matinv and mdivide. Returns the matrix inverse of m.      *
**********************************************************************/
#define INCLUDES APLMEM
#include "includes.h"
double *invert(m,rows)
double *m;
int rows;
{
	int i,j,k,n,p;
	double w,x,*z,*v,*t,*s;

	/* initialize z as an identity matrix */
	z=malloc(rows*rows*sizeof(x));
	t=z+rows*rows;
	for (v=z;v<t;v++)
		*v=0.0; /* all of z = 0 */
	for (v=z;v<t;v+=(rows+1))
		*v=1.0; /* diagonal of z = 1 */
	/* main loop of inverse */
	for (i=0;i<rows;i++) {
		/* loop once for each row */
		n=i*rows; /* common expression */
		/* rearrange rows to give maximum precision */
		#include "inverta.h"
		/* Make the leftmost nonzero element 1 */
		#include "invertb.h"
		/* now make all elements in this column below the diagonal 0 */
		#include "invertc.h"
		/* now make all elements in this column above the diagonal 0 */
		#include "invertd.h"
	}
	return(z); /* z is the inverse of m */
}