MatrixImpl:
CEDAR PROGRAM
IMPORTS
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;
j ← 0; --always use column 0 for now
FOR i
IN [0..nrows)
DO
MakeAij[a,Aij,i,j];
det ← det + a[i][j]*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];
};
Multiply:
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[nra,ncb];
FOR j:
INTEGER
IN [0..nra)
DO
FOR i:
INTEGER
IN [0..ncb)
DO
c[j][i] ← 0;
FOR k:
INTEGER
IN [0..nca)
DO
c[j][i] ← c[j][i] + a[j][k]*b[k][i];
ENDLOOP;
ENDLOOP;
ENDLOOP;
};
MultiplyVec:
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.