/*Copyright (C) 1992, 1995 by Thomas Glen Smith.  All Rights Reserved.*/
/* mdivide APL2 V1.0.0 *************************************************
* mdivide does matrix division with LEFT being the numerator and RIGHT *
* the denominator.  Arguments of rank greater than two are rejected.   *
* A scalar is treated as a 1-by-1 matrix, and a vector as a one-column *
* matrix.                                                              *
***********************************************************************/
#define INCLUDES APLCB+APLMEM
#include "includes.h"
Aplcb mdivide(left,rite)
Aplcb left,rite;
{
	Errinit; Errstop; Getcb; Invert; Matchok; Matmult; Regress; Mdivarg;
	extern int aplerr;
	int datacnt,*dimptr,i,lcols,lrows,rank,rcols,rrows;
	double *dataptr,f,*inv;
	Aplcb out=NULL;

	if (errinit())
		return(errstop(0,left,rite,out));
	if (ERR == mdivarg(rite,&rrows,&rcols))
		return(errstop(0,left,rite,out));
	if (rrows < rcols)
		return(errstop(6,left,rite,out));
	if (ERR == mdivarg(left,&lrows,&lcols))
		return(errstop(0,left,rite,out));
	if (lrows != rrows)
		return(errstop(7,left,rite,out));
	if ((left->aplflags & APLCPLX) || (rite->aplflags & APLCPLX))
		return(errstop(995,left,rite,out)); /* Can't do complex numbers. */
	if (!matchok(&left,&rite,APLNUMB))
		return(out);
	if (rrows == rcols) { /* square matrix? */
		inv = invert(rite->aplptr.apldata,rrows);
		if (NULL == inv)
			return(errstop(8,left,rite,out)); /* not invertible */
		dataptr = matmult(inv,left->aplptr.apldata,rrows,rcols,lcols);
		free(inv);
	}
	else { /* non-square matrix */
		dataptr = regress(left,rite,lrows,lcols,rrows,rcols);
		if (aplerr)
			return(out);
	}
	rank = 0;
	if (left->aplrank>1) rank++;
	if (rite->aplrank>1) rank++;
	datacnt = lcols*rcols;
	out = getcb(dataptr,datacnt,APLTEMP+APLNUMB,rank,NULL);
	if (rank > 1) { /* Build new dimensions? */
		dimptr = out->apldim;
		if (rite->aplrank > 1) *dimptr++ = *(rite->apldim+1);
		if (left->aplrank > 1) *dimptr = *(left->apldim+1);
	}
	return(errstop(0,left,rite,out));
}