--MatrixImpl.mesa
--Maureen Stone July 10, 1983 5:55 pm
DIRECTORY
Matrix;

MatrixImpl: CEDAR PROGRAM
EXPORTS Matrix
= BEGIN OPEN Matrix;
InvalidMatrix: PUBLIC SIGNAL = CODE;
InvalidOperation: PUBLIC SIGNAL = CODE;

ValidMatrix: PROCEDURE [a: MatrixN] RETURNS [nrows,ncols: INTEGER] = {
IF a=NIL OR a.nrows=0 THEN SIGNAL InvalidMatrix;
nrows _ a.nrows;
ncols _ a[0].ncols;
FOR i: INTEGER IN [0..nrows) DO
IF a[i].ncols#ncols THEN SIGNAL InvalidMatrix;
ENDLOOP;
RETURN[nrows,ncols];
};
Sign: PROC[i,j: INTEGER] RETURNS[REAL] = {RETURN[IF (i+j) MOD 2 = 0 THEN 1 ELSE -1]};
MakeAij: PROC[a, Aij: MatrixN, i,j: INTEGER] = {    --assume Aij is well formed
n,m: INTEGER _ 0;    --row index, column index for new matrix
FOR row: INTEGER IN [0..a.nrows) DO
IF row=i THEN LOOP;    --row index for original matrix
m _ 0;
FOR col: INTEGER IN [0..a[0].ncols) DO    --column index for new matrix
IF col=j THEN LOOP;
Aij[n][m] _ a[row][col];    --column 0 of old matrix supplies the aij values
m _ m+1;
ENDLOOP;
n _ n+1;
ENDLOOP;
};


Invert: PUBLIC PROCEDURE [a: MatrixN] RETURNS [ai: MatrixN] = {
nrows,ncols: INTEGER;
det: REAL;
Aij: MatrixN;
[nrows,ncols] _ ValidMatrix[a];
IF nrows#ncols THEN SIGNAL InvalidOperation;
det _ Determinant[a];
ai _ Create[nrows,ncols];
Aij _ Create[nrows-1,ncols-1];
FOR i: INTEGER IN [0..nrows) DO
FOR j: INTEGER IN [0..ncols) DO
MakeAij[a,Aij,i,j];
ai[j][i] _ Sign[i,j]*Determinant[Aij]/det;
ENDLOOP;
ENDLOOP;
RETURN[ai];
};
Determinant: PUBLIC PROC[a: MatrixN] RETURNS [det: REAL] = {
nrows,ncols: INTEGER;
[nrows,ncols] _ ValidMatrix[a];
IF nrows#ncols THEN SIGNAL InvalidOperation;
IF nrows=1 THEN det _ a[0][0]
ELSE IF nrows=2 THEN det _ a[0][0]*a[1][1]-a[0][1]*a[1][0]
ELSE {
i,j: INTEGER;
Aij: MatrixN _ Create[nrows-1,ncols-1];
det _ 0;
FOR i IN [0..nrows) DO
MakeAij[a,Aij,i,0];    --always use column 0
det _ det + a[i][0]*Sign[i,j]*Determinant[Aij];
ENDLOOP;
};
RETURN[det];
};
Create: PUBLIC PROC [nrows, ncols: INTEGER] RETURNS [a: MatrixN] = {
a _ NEW[MatrixSeq[nrows]];
FOR i: INTEGER IN [0..nrows) DO a[i] _ NEW[VecSeq[ncols]]; ENDLOOP;
};
Transpose: PUBLIC PROCEDURE [a: MatrixN] RETURNS [transpose: MatrixN] = {
nrows,ncols: INTEGER;
[nrows,ncols] _ ValidMatrix[a];
transpose _ Create[ncols,nrows];
FOR i: INTEGER IN [0..nrows) DO
FOR j: INTEGER IN [0..ncols) DO
transpose[j][i] _ a[i][j];
ENDLOOP;
ENDLOOP;
RETURN[transpose];
};
Multipy: PUBLIC PROCEDURE [a: MatrixN, b: MatrixN] RETURNS [c: MatrixN] = {
nra,nca: INTEGER;
nrb,ncb: INTEGER;
[nra,nca] _ ValidMatrix[a];
[nrb,ncb] _ ValidMatrix[b];
IF nca#nrb THEN SIGNAL InvalidOperation;
c _ Create[ncb,nra];
FOR j: INTEGER IN [0..nra) DO
FOR i: INTEGER IN [0..ncb) DO
c[i][j] _ 0;
FOR k: INTEGER IN [0..nca) DO
c[i][j] _ c[i][j] + a[j][k]*b[k][i];
ENDLOOP;
ENDLOOP;
ENDLOOP;
};

MultipyVec: PUBLIC PROC[a: MatrixN, v: ColumnN] RETURNS [c: RowN] = {
nrows,ncols,nels: INTEGER;
[nrows,ncols] _ ValidMatrix[a];
IF v=NIL THEN SIGNAL InvalidMatrix;
nels _ v.ncols;    --actually nrows for this case.  historical
IF nels#ncols THEN SIGNAL InvalidOperation;
c _ NEW[VecSeq[nrows]];
FOR i: INTEGER IN [0..nrows) DO
c[i] _ 0;
FOR j: INTEGER IN [0..nels) DO
c[i] _ c[i]+a[i][j]*v[j];
ENDLOOP;
ENDLOOP;
};

END.