MatricesImpl.mesa
Last Edited by: Arnon, June 10, 1985 4:19:22 pm PDT
DIRECTORY
Rope,
Basics,
Atom,
Ascii,
Convert,
IO,
AlgebraClasses,
Matrices;
MatricesImpl: CEDAR PROGRAM
IMPORTS Rope, Atom, Convert, IO
EXPORTS Matrices =
BEGIN OPEN AC: AlgebraClasses, Matrices;
Errors
SyntaxError: PUBLIC ERROR [reason: ATOM] = CODE;
BadElementStructure: PUBLIC ERROR [elementStructure: AC.Structure] = CODE;
CantInvertError: PUBLIC ERROR [elementStructure: AC.Structure] = CODE;
Classes for Matrix Rings and Matrix Algebras
ClassPrintName: AC.PrintNameProc = {
data: MatrixStructureData ← NARROW[structure.instanceData];
RETURN[Rope.Cat[
Convert.RopeFromCard[data.size],
" x ",
Convert.RopeFromCard[data.size],
" Matrices over ",
data.elementStructure.class.printName[data.elementStructure]
] ];
};
ClassCharacteristic: AC.StructureRankOp = {
data: MatrixStructureData ← NARROW[structure.instanceData];
RETURN[ data.elementStructure.class.characteristic[data.elementStructure] ]
};
ClassDimension: AC.StructureRankOp = {
data: MatrixStructureData ← NARROW[structure.instanceData];
RETURN[data.size * data.size]
};
ClassLegalFirstChar: AC.LegalFirstCharOp = {
SELECT char FROM
'[ => RETURN[TRUE];
ENDCASE;
RETURN[FALSE];
};
ClassRead: AC.ReadOp = {
data: MatrixStructureData ← NARROW[structure.instanceData];
RETURN[ReadMatrix[in, data.size, data.elementStructure] ];
};
ClassFromRope: AC.FromRopeOp = {
stream: IO.STREAM ← IO.RIS[in];
RETURN[ ClassRead[stream, structure] ];
};
ClassToRope: AC.ToRopeOp = {
matrix: Matrix ← NARROW[in];
data: MatrixStructureData ← NARROW[structure.instanceData];
RETURN[ MatrixToRope[matrix, data.elementStructure] ] -- use defaults
};
ClassWrite: AC.WriteOp = {
matrix: Matrix ← NARROW[in];
IO.PutRope[stream, ClassToRope[matrix, structure] ]
};
ClassZero: AC.NullaryOp = {
data: MatrixStructureData ← NARROW[structure.instanceData];
RETURN[ DiagMatrix[data.elementStructure.class.zero[data.elementStructure], data.size, data.elementStructure] ]
};
ClassOne: AC.NullaryOp = {
data: MatrixStructureData ← NARROW[structure.instanceData];
RETURN[ DiagMatrix[data.elementStructure.class.one[data.elementStructure], data.size, data.elementStructure] ]
};
ClassAdd: AC.BinaryOp = {
firstMatrix: Matrix ← NARROW[firstArg];
secondMatrix: Matrix ← NARROW[secondArg];
data: MatrixStructureData ← NARROW[structure.instanceData];
RETURN[ Add[firstMatrix, secondMatrix, data.elementStructure] ]
};
ClassNegate: AC.UnaryOp = {
matrix: Matrix ← NARROW[arg];
data: MatrixStructureData ← NARROW[structure.instanceData];
RETURN[ Negate[matrix, data.elementStructure] ]
};
ClassSubtract: AC.BinaryOp = {
firstMatrix: Matrix ← NARROW[firstArg];
secondMatrix: Matrix ← NARROW[secondArg];
data: MatrixStructureData ← NARROW[structure.instanceData];
RETURN[ Subtract[firstMatrix, secondMatrix, data.elementStructure] ]
};
ClassMultiply: AC.BinaryOp = {
firstMatrix: Matrix ← NARROW[firstArg];
secondMatrix: Matrix ← NARROW[secondArg];
data: MatrixStructureData ← NARROW[structure.instanceData];
RETURN[ Multiply[firstMatrix, secondMatrix, data.elementStructure] ]
};
ClassEqual: AC.EqualityOp = {
firstMatrix: Matrix ← NARROW[firstArg];
secondMatrix: Matrix ← NARROW[secondArg];
data: MatrixStructureData ← NARROW[structure.instanceData];
RETURN[ Equal[firstMatrix, secondMatrix, data.elementStructure] ]
};
ClassTranspose: AC.UnaryOp = {
matrix: Matrix ← NARROW[arg];
RETURN[ Transp[matrix] ]
};
ClassDeterminant: AC.StructuredToGroundOp = {
matrix: Matrix ← NARROW[structuredElt];
data: MatrixStructureData ← NARROW[structure.instanceData];
RETURN[ Det[matrix, data.elementStructure] ]
};
matrixRingOps: MatrixOps ← NEW[MatrixOpsRec ← [
diagonalMatrix: DiagMatrix,
transpose: ClassTranspose,
determinant: ClassDeterminant
] ];
matrixProp: Atom.DottedPair ← NEW[Atom.DottedPairNode← [$MatrixRing, matrixRingOps]]; -- having the $MatrixRing property actually means being a matrix structure but not an algebra
matrixRingClass: PUBLIC AC.StructureClass ← NEW[AC.StructureClassRec ← [
flavor: ring,
printName: ClassPrintName,
characteristic: ClassCharacteristic,
legalFirstChar: ClassLegalFirstChar,
read: ClassRead,
fromRope: ClassFromRope,
toRope: ClassToRope,
write: ClassWrite,
add: ClassAdd,
negate: ClassNegate,
subtract: ClassSubtract,
zero: ClassZero,
multiply: ClassMultiply,
commutative: FALSE,
one: ClassOne,
equal: ClassEqual,
integralDomain: FALSE,
gcdDomain: FALSE,
euclideanDomain: FALSE,
propList: LIST[matrixProp]
] ];
ClassInvert: AC.UnaryOp = {
matrix: Matrix ← NARROW[arg];
data: MatrixStructureData ← NARROW[structure.instanceData];
RETURN[ Invert[matrix, data.elementStructure] ]
};
ClassDivide: AC.BinaryOp = {
firstMatrix: Matrix ← NARROW[firstArg];
secondMatrix: Matrix ← NARROW[secondArg];
data: MatrixStructureData ← NARROW[structure.instanceData];
RETURN[ Divide[firstMatrix, secondMatrix, data.elementStructure] ]
};
ClassScalarMultiply: AC.BinaryOp = {
scalar: REF ← firstArg;
inMatrix: Matrix ← NARROW[secondArg];
data: MatrixStructureData ← NARROW[structure.instanceData];
RETURN[ ScalarMultiply[scalar, inMatrix, data.elementStructure] ]
};
matrixAlgebraClass: PUBLIC AC.StructureClass ← NEW[AC.StructureClassRec ← [
flavor: algebra,
printName: ClassPrintName,
characteristic: ClassCharacteristic,
dimension: ClassDimension,
legalFirstChar: ClassLegalFirstChar,
read: ClassRead,
fromRope: ClassFromRope,
toRope: ClassToRope,
write: ClassWrite,
add: ClassAdd,
negate: ClassNegate,
subtract: ClassSubtract,
zero: ClassZero,
multiply: ClassMultiply,
commutative: FALSE,
invert: ClassInvert, -- it is up to the user not to try to invert a singular matrix
divide: ClassDivide,
one: ClassOne,
scalarMultiply: ClassScalarMultiply,
equal: ClassEqual,
integralDomain: FALSE,
gcdDomain: FALSE,
euclideanDomain: FALSE,
propList: LIST[matrixProp]
] ];
Matrix Ring and Algebra Constructors
MakeMatrixRing: PUBLIC PROC [elementStructure: AC.Structure, size: NAT] RETURNS [matrixRing: AC.Structure] ~ {
matrixRingData: MatrixStructureData ← NEW[MatrixStructureDataRec ← [
elementStructure: elementStructure,
size: size
] ];
IF elementStructure.class.flavor#ring AND elementStructure.class.flavor#algebra THEN ERROR BadElementStructure[elementStructure];
RETURN[ NEW[AC.StructureRec ← [
class: matrixRingClass,
instanceData: matrixRingData
] ] ];
};
MakeMatrixAlgebra: PUBLIC PROC [elementField: AC.Structure, size: NAT] RETURNS [matrixAlgebra: AC.Structure] ~ {
matrixAlgebraData: MatrixStructureData ← NEW[MatrixStructureDataRec ← [
elementStructure: elementField,
size: size
] ];
IF elementField.class.flavor#field THEN ERROR BadElementStructure[elementField];
RETURN[ NEW[AC.StructureRec ← [
class: matrixAlgebraClass,
instanceData: matrixAlgebraData
] ] ];
};
Extract Matrix Operations from Class Record Property Lists
IsMatrixRing: PUBLIC PROC [structure: AC.Structure] RETURNS [BOOL] ~ {
IF Atom.GetPropFromList[structure.class.propList, $MatrixRing] # NIL THEN RETURN[TRUE] ELSE RETURN[FALSE];
};
Transpose: PUBLIC PROC [structure: AC.Structure] RETURNS [AC.UnaryOp] ~ {
IF IsMatrixRing[structure] THEN {
matrixOps: MatrixOps ← NARROW[ Atom.GetPropFromList[structure.class.propList, $MatrixRing] ];
RETURN[matrixOps.transpose];
}
ELSE ERROR;
};
Determinant: PUBLIC PROC [structure: AC.Structure] RETURNS [AC.StructuredToGroundOp] ~ {
IF IsMatrixRing[structure] THEN {
matrixOps: MatrixOps ← NARROW[ Atom.GetPropFromList[structure.class.propList, $MatrixRing] ];
RETURN[matrixOps.determinant];
}
ELSE ERROR;
};
Constructors
DiagMatrix: PUBLIC PROC [element: REF, size: NAT, elementStructure: AC.Structure] RETURNS [out: Matrix] ~ {
rows: RowSeq ← NEW[RowSeqRec[size]];
FOR i:NAT IN [1..size] DO
rows[i] ← MostlyZeroRow[element, i, size, elementStructure];
ENDLOOP;
out ← NEW[MatrixRec ← [size: size, rows: rows] ];
};
MostlyZeroRow: PROC [element: REF, position, size: NAT, elementStructure: AC.Structure] RETURNS [row: Row] ~ {
zeros everywhere except element in position
zero: REF ← elementStructure.class.zero[elementStructure];
row ← NEW[RowRec[size]];
FOR i:NAT IN [1..size] DO
row[i] ← IF i = position THEN element ELSE zero;
ENDLOOP;
};
Conversion and IO
ReadMatrix: PUBLIC PROC [in: IO.STREAM, size: NAT, elementStructure: AC.Structure] RETURNS [matrix: Matrix] = {
puncChar: CHAR;
rows: RowSeq ← NEW[RowSeqRec[size]];
[]← in.SkipWhitespace[];
puncChar ← in.GetChar[ ];
IF puncChar # '[ THEN ERROR SyntaxError[$LeftBracketExpected];
FOR i:NAT IN [1..size] DO
row: Row ← NEW[RowRec[size]];
[]← in.SkipWhitespace[];
puncChar ← in.GetChar[ ];
IF puncChar # '[ THEN ERROR SyntaxError[$LeftBracketExpected];
FOR j:NAT IN [1..size] DO
row[j] ← elementStructure.class.read[in, elementStructure];
[]← in.SkipWhitespace[];
puncChar ← in.GetChar[];
[]← in.SkipWhitespace[];
IF j < size THEN IF puncChar # ', THEN ERROR SyntaxError[$CommaExpected];
ENDLOOP;
IF puncChar # '] THEN ERROR SyntaxError[$RightBracketExpected];
rows[i] ← row;
puncChar ← in.GetChar[];
[]← in.SkipWhitespace[];
IF i < size THEN IF puncChar # ', THEN SyntaxError[$CommaExpected];
ENDLOOP;
IF puncChar # '] THEN ERROR SyntaxError[$RightBracketExpected];
matrix ← NEW[MatrixRec ← [size: size, rows: rows] ];
};
MatrixFromRope: PUBLIC PROC [in: Rope.ROPE, size: NAT, elementStructure: AC.Structure] RETURNS [out: Matrix] = {
stream: IO.STREAM ← IO.RIS[in];
out ← ReadMatrix[stream, size, elementStructure];
};
MatrixToRope: PUBLIC PROC [in: Matrix, elementStructure: AC.Structure] RETURNS [out: Rope.ROPE] = {
out ← "[ ";
FOR i:NAT IN [1..in.size] DO
out ← out.Concat["[ "];
FOR j:NAT IN [1..in.size] DO
out ← out.Concat[elementStructure.class.toRope[in.rows[i][j], elementStructure] ];
IF j < in.size THEN out ← out.Concat[", "];
ENDLOOP;
out ← out.Concat["] "];
IF i < in.size THEN out ← out.Concat[", "];
ENDLOOP;
out ← out.Concat["] "];
};
WriteMatrix: PUBLIC PROC [in: Matrix, elementStructure: AC.Structure, out: IO.STREAM] = {
out.PutF["\n %g \n", IO.rope[MatrixToRope[in, elementStructure]] ];
};
Arithmetic
Create: PROC [size: NAT] RETURNS [a: Matrix] = {
a ← NEW[MatrixRec ← [
size: size,
rows: NEW[RowSeqRec[size]]
] ];
FOR i: NAT IN [1..size] DO a.rows[i] ← NEW[RowRec[size]]; ENDLOOP;
};
Add: PUBLIC PROC [in1, in2: Matrix, elementStructure: AC.Structure] RETURNS [out: Matrix] ~ {
out ← Create[in1.size];
FOR i:NAT IN [1..in1.size] DO
FOR j:NAT IN [1..in1.size] DO
out.rows[i][j] ← elementStructure.class.add[in1.rows[i][j], in2.rows[i][j], elementStructure];
ENDLOOP;
ENDLOOP;
};
Negate: PUBLIC PROC [in: Matrix, elementStructure: AC.Structure] RETURNS [out: Matrix] ~ {
out ← Create[in.size];
FOR i:NAT IN [1..in.size] DO
FOR j:NAT IN [1..in.size] DO
out.rows[i][j] ← elementStructure.class.negate[in.rows[i][j], elementStructure];
ENDLOOP;
ENDLOOP;
};
Subtract: PUBLIC PROC [in1, in2: Matrix, elementStructure: AC.Structure] RETURNS [Matrix] ~ {
RETURN[ Add[ in1, Negate[ in2, elementStructure], elementStructure ] ];
};
Multiply: PUBLIC PROC [in1, in2: Matrix, elementStructure: AC.Structure] RETURNS [out: Matrix] ~ {
zero: REF ← elementStructure.class.zero[elementStructure];
out ← Create[in1.size];
FOR i:NAT IN [1..in1.size] DO
FOR j:NAT IN [1..in1.size] DO
out.rows[i][j] ← zero;
FOR k:NAT IN [1..in1.size] DO
out.rows[i][j] ← elementStructure.class.add[
out.rows[i][j],
elementStructure.class.multiply[in1.rows[i][k], in2.rows[k][j], elementStructure],
elementStructure]
ENDLOOP;
ENDLOOP;
ENDLOOP;
};
Invert: PUBLIC PROC [in: Matrix, elementField: AC.Structure] RETURNS [out: Matrix] ~ {
det: REF ← Det[in, elementField];
tempVal: REF;
Aij: Matrix ← Create[in.size-1];
IF elementField.class.flavor # field THEN ERROR CantInvertError[elementField];
out ← Create[in.size];
FOR i: NAT IN [1..in.size] DO
FOR j: NAT IN [1..in.size] DO
MakeAij[in, Aij, i, j];
tempVal ← elementField.class.divide[Det[Aij, elementField], det, elementField]; -- error will be generated here if Det[Aij] = 0
IF NOT Even[i + j] THEN tempVal ← elementField.class.negate[tempVal, elementField];
out.rows[j][i] ← tempVal;
ENDLOOP;
ENDLOOP;
RETURN[out];
};
Divide: PUBLIC PROC [in1, in2: Matrix, elementField: AC.Structure] RETURNS [out: Matrix]~ {
RETURN[ Multiply[ in1, Invert[ in2, elementField], elementField ] ];
};
ScalarMultiply: PUBLIC PROC [scalar: REF, inMatrix: Matrix, elementField: AC.Structure] RETURNS [outMatrix: Matrix]~ {
outMatrix ← Create[inMatrix.size];
FOR i:NAT IN [1..inMatrix.size] DO
FOR j:NAT IN [1..inMatrix.size] DO
outMatrix.rows[i][j] ← elementField.class.multiply[scalar, inMatrix.rows[i][j], elementField];
ENDLOOP;
ENDLOOP;
};
Matrix Operations
Equal: PUBLIC PROC [in1, in2: Matrix, elementStructure: AC.Structure] RETURNS [BOOL ← TRUE] ~ {
diff: Matrix ← Subtract[in1, in2, elementStructure];
FOR i: INTEGER IN [1..diff.size] DO
FOR j: INTEGER IN [1..diff.size] DO
IF NOT elementStructure.class.equal[diff.rows[i][j], elementStructure.class.zero[elementStructure], elementStructure] THEN RETURN[FALSE];
ENDLOOP;
ENDLOOP;
};
Transp: PUBLIC PROC [in: Matrix] RETURNS [out: Matrix] ~ {
out ← Create[in.size];
FOR i: INTEGER IN [1..in.size] DO
FOR j: INTEGER IN [1..in.size] DO
out.rows[j][i] ← in.rows[i][j];
ENDLOOP;
ENDLOOP;
RETURN[out];
};
MakeAij: PROC[a, Aij: Matrix, i,j: INTEGER] = { --assume Aij is well formed
n, m: NAT ← 1; --row index, column index for new matrix
FOR row: NAT IN [1..a.size] DO
IF row = i THEN LOOP; --row index for original matrix
m ← 1;
FOR col: NAT IN [1..a.size] DO --column index for original matrix
IF col = j THEN LOOP;
Aij.rows[n][m] ← a.rows[row][col];
m ← m+1;
ENDLOOP;
n ← n+1;
ENDLOOP;
};
Even: PROC[I: NAT] RETURNS[BOOL] = {RETURN[I MOD 2 = 0]};
Det: PUBLIC PROC [in: Matrix, elementStructure: AC.Structure] RETURNS [value: REF] ~ {
Expansion by minors.
IF in.size = 1 THEN value ← in.rows[1][1]
ELSE IF in.size=2 THEN value ← elementStructure.class.subtract[
elementStructure.class.multiply[in.rows[1][1], in.rows[2][2], elementStructure],
elementStructure.class.multiply[in.rows[1][2], in.rows[2][1], elementStructure],
elementStructure]
ELSE {
i, j: NAT;
zero: REF ← elementStructure.class.zero[elementStructure];
Aij: Matrix ← Create[in.size-1];
value ← zero;
j ← 1; --always use column 1 for now
FOR i IN [1..in.size] DO
tempVal: REF;
MakeAij[in, Aij, i, j];
tempVal ← elementStructure.class.multiply[in.rows[i][j], Det[Aij, elementStructure], elementStructure];
IF NOT Even[i + j] THEN tempVal ← elementStructure.class.negate[tempVal, elementStructure];
value ← elementStructure.class.add[value, tempVal, elementStructure];
ENDLOOP;
};
RETURN[value];
};
END.