MODULE MtxInt;
IMPORT SYSTEM, NbrInt, ArrayXdBytes, ArrayXd := ArrayXdInt, Array1d:= Array1dInt, DataErrors, Vec := VecInt, DataIO, NbrRe;
CONST
VERSION* = 1;
TYPE
Value* = ArrayXd.Value; Index* = LONGINT; Array* = ArrayXd.Array2; Map* = ArrayXd.Map;
Matrix* = OBJECT (ArrayXd.Array)
VAR lenx-, leny-: LONGINT;
rows-, cols-: LONGINT;
VAR ox-, oy-: LONGINT;
Get-: PROCEDURE {DELEGATE} ( x, y: Index ): Value;
PROCEDURE AlikeX( ): ArrayXdBytes.Array;
VAR copy: Matrix;
BEGIN
NEW( copy, origin[0], len[0], origin[1], len[1] ); RETURN copy;
END AlikeX;
PROCEDURE NewRangeX( neworigin, newlen: ArrayXdBytes.IndexArray; copydata: BOOLEAN );
BEGIN
IF LEN( newlen ) # 2 THEN HALT( 1001 ) END;
NewRangeX^( neworigin, newlen, copydata );
END NewRangeX;
PROCEDURE ValidateCache;
BEGIN
ValidateCache^;
IF dim # 2 THEN HALT( 100 ) END;
lenx := len[0]; leny := len[1]; ox := origin[0]; oy := origin[1]; rows := lenx; cols := leny;
END ValidateCache;
PROCEDURE SetBoundaryCondition*( c: SHORTINT );
BEGIN
SetBoundaryCondition^( c );
CASE c OF
ArrayXd.StrictBoundaryC:
Get := Get2;
| ArrayXd.AbsorbingBoundaryC:
Get := Get2BAbsorbing;
| ArrayXd.PeriodicBoundaryC:
Get := Get2BPeriodic;
| ArrayXd.SymmetricOnBoundaryC:
Get := Get2BSymmetricOnB
| ArrayXd.SymmetricOffBoundaryC:
Get := Get2BSymmetricOffB
| ArrayXd.AntisymmetricOnBoundaryC:
Get := Get2BAntisymmetricOnB
| ArrayXd.AntisymmetricOffBoundaryC:
Get := Get2BAntisymmetricOffB
END;
END SetBoundaryCondition;
PROCEDURE & New*( ox, rowsORw, oy, colsORh: LONGINT );
BEGIN
NewXdB( ArrayXdBytes.Array2( ox, oy ), ArrayXdBytes.Array2( rowsORw, colsORh ) );
END New;
PROCEDURE Alike*( ): Matrix;
VAR copy: ArrayXdBytes.Array;
BEGIN
copy := AlikeX(); RETURN copy( Matrix );
END Alike;
PROCEDURE NewRange*( ox, rowsORw, oy, colsORh: LONGINT; copydata: BOOLEAN );
BEGIN
IF (rowsORw # len[0]) OR (colsORh # len[1]) OR (ox # origin[0]) OR (oy # origin[1]) THEN
NewRangeX^( ArrayXdBytes.Array2( ox, oy ), ArrayXdBytes.Array2( rowsORw, colsORh ), copydata )
END;
END NewRange;
PROCEDURE Copy*( ): Matrix;
VAR res: ArrayXdBytes.Array;
BEGIN
res := CopyX(); RETURN res( Matrix );
END Copy;
PROCEDURE Set*( rowORx , colORy : Index; v: Value );
BEGIN
ArrayXdBytes.Set2( SELF, rowORx, colORy, v );
END Set;
PROCEDURE SwapRows*( row1, row2: Index );
BEGIN
ToggleElements( 0, row1, row2 )
END SwapRows;
PROCEDURE SwapColumns*( col1, col2: Index );
BEGIN
ToggleElements( 1, col1, col2 )
END SwapColumns;
PROCEDURE Transpose*;
BEGIN
PermuteDimensions( ArrayXdBytes.Array2( 1, 0 ), TRUE );
END Transpose;
PROCEDURE Dot*( x: Matrix );
VAR res: Matrix;
BEGIN
IF x # NIL THEN
IF x.len[0] = x.len[1] THEN
IF len[0] # len[1] THEN
IF len[0] = x.len[0] THEN
res := SELF * x; ArrayXdBytes.CopyArrayPartToArrayPart( res, SELF, res.origin, res.len, origin, len );
ELSE DataErrors.Error( "The two matrices were not compatible" )
END;
ELSE DataErrors.Error( "This matrix in not square." )
END
ELSE DataErrors.Error( "The supplied matrix was not square." )
END
ELSE DataErrors.Error( "The supplied matrix to be doted was NIL." )
END
END Dot;
PROCEDURE LeftDot*( x: Matrix );
BEGIN
Transpose; Dot( x );
END LeftDot;
PROCEDURE RightDot*( x: Matrix );
BEGIN
x.Transpose; Dot( x ); x.Transpose;
END RightDot;
PROCEDURE Row*( row: Index ): Vec.Vector;
VAR v: Vec.Vector;
BEGIN
NEW( v, 0, len[1] ); ArrayXd.CopyMtxToVec( SELF, v, 1, row, 0, 0, len[1] ); RETURN v;
END Row;
PROCEDURE InsertRow*( at: Index );
BEGIN
InsertElements( 0, at, 1 )
END InsertRow;
PROCEDURE DeleteRow*( x: Index );
BEGIN
DeleteElements( 0, x, 1 );
END DeleteRow;
PROCEDURE Col*( col: Index ): Vec.Vector;
VAR v: Vec.Vector;
BEGIN
NEW( v, 0, len[0] ); ArrayXd.CopyMtxToVec( SELF, v, 0, 0, col, 0, len[0] ); RETURN v;
END Col;
PROCEDURE InsertCol*( at: Index );
BEGIN
InsertElements( 1, at, 1 )
END InsertCol;
PROCEDURE DeleteCol*( x: Index );
BEGIN
DeleteElements( 1, x, 1 );
END DeleteCol;
PROCEDURE CopyToVec*( dest: ArrayXd.Array; dim: Index; srcx, srcy, destx, len: Index );
VAR slen: ArrayXdBytes.IndexArray;
BEGIN
IF (dest.dim # 1) THEN HALT( 1002 ) END;
slen := ArrayXdBytes.Index2( 1, 1 ); slen[dim] := len;
CopyToArray( dest, ArrayXdBytes.Index2( srcx, srcy ), slen, ArrayXdBytes.Index1( destx ), ArrayXdBytes.Index1( len ) );
END CopyToVec;
PROCEDURE CopyToMtx*( dest: ArrayXd.Array; srcx, srcy, destx, desty, lenx, leny: Index );
VAR slen: ArrayXdBytes.IndexArray;
BEGIN
IF (dest.dim # 2) THEN HALT( 1005 ) END;
slen := ArrayXdBytes.Index2( lenx, leny );
CopyToArray( dest, ArrayXdBytes.Index2( srcx, srcy ), slen, ArrayXdBytes.Index2( destx, desty ), slen );
END CopyToMtx;
PROCEDURE CopyToCube*( dest: ArrayXd.Array; dimx, dimy: Index; srcx, srcy, destx, desty, destz, lenx, leny: Index );
VAR slen: ArrayXdBytes.IndexArray;
BEGIN
IF (dest.dim # 3) OR (dimx >= dimy) THEN HALT( 1005 ) END;
slen := ArrayXdBytes.Index3( 1, 1, 1 ); slen[dimx] := lenx; slen[dimy] := leny;
CopyToArray( dest, ArrayXdBytes.Index2( srcx, srcy ), ArrayXdBytes.Index2( lenx, leny ),
ArrayXdBytes.Index3( destx, desty, destz ), slen );
END CopyToCube;
PROCEDURE CopyToHCube*( dest: ArrayXd.Array; dimx, dimy: Index;
srcx, srcy, destx, desty, destz, destt, lenx, leny: Index );
VAR slen: ArrayXdBytes.IndexArray;
BEGIN
IF (dest.dim # 4) OR (dimx >= dimy) THEN HALT( 1005 ) END;
slen := ArrayXdBytes.Index4( 1, 1, 1, 1 ); slen[dimx] := lenx; slen[dimy] := leny;
CopyToArray( dest, ArrayXdBytes.Index2( srcx, srcy ), ArrayXdBytes.Index2( lenx, leny ),
ArrayXdBytes.Index4( destx, desty, destz, destt ), slen );
END CopyToHCube;
PROCEDURE CopyTo1dArray*( VAR dest: ARRAY OF Value; sx, sy, slenx, sleny: Index; dpos, dlen: LONGINT );
VAR destm: ArrayXdBytes.ArrayMemoryStructure;
BEGIN
destm :=
ArrayXdBytes.MakeMemoryStructure( 1, ArrayXdBytes.Index1( 0 ), ArrayXdBytes.Index1( LEN( dest ) ), SYSTEM.SIZEOF( Value ),
SYSTEM.ADR( dest[0] ) );
ArrayXd.CopyArrayToArrayPartB( SELF, destm, bc, ArrayXdBytes.Index2( sx, sy ), ArrayXdBytes.Index2( slenx, sleny ),
ArrayXdBytes.Index1( dpos ), ArrayXdBytes.Index1( dlen ) );
END CopyTo1dArray;
PROCEDURE CopyTo2dArray*( VAR dest: ARRAY OF ARRAY OF Value; sx, sy, slenx, sleny: Index; dposx, dposy, dlenx, dleny: LONGINT );
VAR destm: ArrayXdBytes.ArrayMemoryStructure;
BEGIN
destm :=
ArrayXdBytes.MakeMemoryStructure( 2, ArrayXdBytes.Index2( 0, 0 ), ArrayXdBytes.Index2( LEN( dest, 1 ), LEN( dest, 0 ) ),
SYSTEM.SIZEOF( Value ), SYSTEM.ADR( dest[0, 0] ) );
ArrayXd.CopyArrayToArrayPartB( SELF, destm, bc, ArrayXdBytes.Index2( sx, sy ), ArrayXdBytes.Index2( slenx, sleny ),
ArrayXdBytes.Index2( dposx, dposy ), ArrayXdBytes.Index2( dlenx, dleny ) );
END CopyTo2dArray;
PROCEDURE CopyTo3dArray*( VAR dest: ARRAY OF ARRAY OF ARRAY OF Value; sx, sy, slenx, sleny: Index;
dposx, dposy, dposz, dlenx, dleny, dlenz: LONGINT );
VAR destm: ArrayXdBytes.ArrayMemoryStructure;
BEGIN
destm :=
ArrayXdBytes.MakeMemoryStructure( 3, ArrayXdBytes.Index3( 0, 0, 0 ),
ArrayXdBytes.Index3( LEN( dest, 2 ), LEN( dest, 1 ), LEN( dest, 0 ) ), SYSTEM.SIZEOF( Value ),
SYSTEM.ADR( dest[0, 0, 0] ) );
ArrayXd.CopyArrayToArrayPartB( SELF, destm, bc, ArrayXdBytes.Index2( sx, sy ), ArrayXdBytes.Index2( slenx, sleny ),
ArrayXdBytes.Index3( dposx, dposy, dposz ),
ArrayXdBytes.Index3( dlenx, dleny, dlenz ) );
END CopyTo3dArray;
PROCEDURE CopyTo4dArray*( VAR dest: ARRAY OF ARRAY OF ARRAY OF ARRAY OF Value; sx, sy, slenx, sleny: Index;
dposx, dposy, dposz, dpost, dlenx, dleny, dlenz, dlent: LONGINT );
VAR destm: ArrayXdBytes.ArrayMemoryStructure;
BEGIN
destm :=
ArrayXdBytes.MakeMemoryStructure( 4, ArrayXdBytes.Index4( 0, 0, 0, 0 ),
ArrayXdBytes.Index4( LEN( dest, 3 ), LEN( dest, 2 ), LEN( dest, 1 ), LEN( dest, 0 ) ), SYSTEM.SIZEOF( Value ),
SYSTEM.ADR( dest[0, 0, 0, 0] ) );
ArrayXd.CopyArrayToArrayPartB( SELF, destm, bc, ArrayXdBytes.Index2( sx, sy ), ArrayXdBytes.Index2( slenx, sleny ),
ArrayXdBytes.Index4( dposx, dposy, dposz, dpost ),
ArrayXdBytes.Index4( dlenx, dleny, dlenz, dlent ) );
END CopyTo4dArray;
PROCEDURE CopyFrom1dArray*( VAR src: ARRAY OF Value; spos, slen: Index; dx, dy, dlenx, dleny: Index );
VAR srcm: ArrayXdBytes.ArrayMemoryStructure;
BEGIN
srcm :=
ArrayXdBytes.MakeMemoryStructure( 1, ArrayXdBytes.Index1( 0 ), ArrayXdBytes.Index1( LEN( src ) ), SYSTEM.SIZEOF( Value ),
SYSTEM.ADR( src[0] ) );
ArrayXdBytes.CopyArrayPartToArrayPart( srcm, SELF, ArrayXdBytes.Index1( spos ), ArrayXdBytes.Index1( slen ),
ArrayXdBytes.Index2( dx, dy ), ArrayXdBytes.Index2( dlenx, dleny ) );
END CopyFrom1dArray;
PROCEDURE CopyFrom2dArray*( VAR src: ARRAY OF ARRAY OF Value; sposx, spoxy, slenx, sleny: Index;
dx, dy, dlenx, dleny: Index );
VAR srcm: ArrayXdBytes.ArrayMemoryStructure;
BEGIN
srcm :=
ArrayXdBytes.MakeMemoryStructure( 2, ArrayXdBytes.Index2( 0, 0 ), ArrayXdBytes.Index2( LEN( src, 1 ), LEN( src, 0 ) ),
SYSTEM.SIZEOF( Value ), SYSTEM.ADR( src[0, 0] ) );
ArrayXdBytes.CopyArrayPartToArrayPart( srcm, SELF, ArrayXdBytes.Index2( sposx, spoxy ),
ArrayXdBytes.Index2( slenx, sleny ), ArrayXdBytes.Index2( dx, dy ),
ArrayXdBytes.Index2( dlenx, dleny ) );
END CopyFrom2dArray;
PROCEDURE CopyFrom3dArray*( VAR src: ARRAY OF ARRAY OF ARRAY OF Value; sposx, spoxy, sposz, slenx, sleny, slenz: Index;
dx, dy, dlenx, dleny: Index );
VAR srcm: ArrayXdBytes.ArrayMemoryStructure;
BEGIN
srcm :=
ArrayXdBytes.MakeMemoryStructure( 3, ArrayXdBytes.Index3( 0, 0, 0 ),
ArrayXdBytes.Index3( LEN( src, 2 ), LEN( src, 1 ), LEN( src, 0 ) ), SYSTEM.SIZEOF( Value ),
SYSTEM.ADR( src[0, 0, 0] ) );
ArrayXdBytes.CopyArrayPartToArrayPart( srcm, SELF, ArrayXdBytes.Index3( sposx, spoxy, sposz ),
ArrayXdBytes.Index3( slenx, sleny, slenz ), ArrayXdBytes.Index2( dx, dy ),
ArrayXdBytes.Index2( dlenx, dleny ) );
END CopyFrom3dArray;
PROCEDURE CopyFrom4dArray*( VAR src: ARRAY OF ARRAY OF ARRAY OF ARRAY OF Value;
sposx, spoxy, sposz, spost, slenx, sleny, slenz, slent: Index; dx, dy, dlenx, dleny: Index );
VAR srcm: ArrayXdBytes.ArrayMemoryStructure;
BEGIN
srcm :=
ArrayXdBytes.MakeMemoryStructure( 4, ArrayXdBytes.Index4( 0, 0, 0, 0 ),
ArrayXdBytes.Index4( LEN( src, 3 ), LEN( src, 2 ), LEN( src, 1 ), LEN( src, 0 ) ), SYSTEM.SIZEOF( Value ),
SYSTEM.ADR( src[0, 0, 0, 0] ) );
ArrayXdBytes.CopyArrayPartToArrayPart( srcm, SELF, ArrayXdBytes.Index4( sposx, spoxy, sposz, spost ),
ArrayXdBytes.Index4( slenx, sleny, slenz, slent ),
ArrayXdBytes.Index2( dx, dy ), ArrayXdBytes.Index2( dlenx, dleny ) );
END CopyFrom4dArray;
END Matrix;
PROCEDURE FrobeniusNorm*( m: Matrix ): NbrRe.Real;
BEGIN
RETURN Array1d.L2Norm( m.data^, 0, LEN( m.data ) );
END FrobeniusNorm;
PROCEDURE Transpose*( u: Matrix ): Matrix;
VAR res: Matrix;
BEGIN
res := u.Copy(); res.Transpose; RETURN res;
END Transpose;
PROCEDURE ":="*( VAR l: Matrix; VAR r: ARRAY OF ARRAY OF Value );
BEGIN
IF l = NIL THEN NEW( l, 0, LEN( r, 1 ), 0, LEN( r, 0 ) ) ELSE l.NewRange( 0, LEN( r, 1 ), 0, LEN( r, 0 ), FALSE ); END;
ArrayXdBytes.CopyMemoryToArrayPart( SYSTEM.ADR( r[0, 0] ), l, LEN( r, 1 ) * LEN( r, 0 ), NIL , NIL );
END ":=";
PROCEDURE ":="*( VAR l: Matrix; r: Vec.Vector );
BEGIN
IF l = NIL THEN NEW( l, 0, 1, 0, r.len[0] ) ELSE l.NewRange( 0, 1, 0, r.len[0], FALSE ); END;
ArrayXdBytes.CopyDataRaw( r, l );
END ":=";
PROCEDURE ":="*( VAR l: Matrix; r: Value );
BEGIN
IF l # NIL THEN ArrayXd.Fill( l, r ); END;
END ":=";
PROCEDURE "+"*( l, r: Matrix ): Matrix;
VAR res: Matrix;
BEGIN
res := l.Alike(); ArrayXd.Add( l, r, res ); RETURN res;
END "+";
PROCEDURE "-"*( l, r: Matrix ): Matrix;
VAR res: Matrix;
BEGIN
res := l.Alike(); ArrayXd.Sub( l, r, res ); RETURN res;
END "-";
PROCEDURE "+"*( l: Matrix; r: Value ): Matrix;
VAR res: Matrix;
BEGIN
res := l.Alike(); ArrayXd.AddAV( l, r, res ); RETURN res;
END "+";
PROCEDURE "+"*( l: Value; r: Matrix ): Matrix;
BEGIN
RETURN r + l
END "+";
PROCEDURE "-"*( l: Matrix; r: Value ): Matrix;
VAR res: Matrix;
BEGIN
res := l.Alike(); ArrayXd.SubAV( l, r, res ); RETURN res;
END "-";
PROCEDURE "-"*( l: Value; r: Matrix ): Matrix;
VAR res: Matrix;
BEGIN
res := r.Alike(); ArrayXd.SubVA( l, r, res ); RETURN res;
END "-";
PROCEDURE "-"*( l: Matrix ): Matrix;
BEGIN
RETURN 0 - l;
END "-";
PROCEDURE "*"*( l: Matrix; r: Value ): Matrix;
VAR res: Matrix;
BEGIN
res := l.Alike(); ArrayXd.MulAV( l, r, res ); RETURN res;
END "*";
PROCEDURE "*"*( l: Value; r: Matrix ): Matrix;
BEGIN
RETURN r * l;
END "*";
PROCEDURE "DIV"*( l: Matrix; r: Value ): Matrix;
VAR res: Matrix;
BEGIN
res := l.Alike(); ArrayXd.DivAV( l, r, res ); RETURN res;
END "DIV";
PROCEDURE "DIV"*( l: Value; r: Matrix ): Matrix;
VAR res: Matrix;
BEGIN
res := r.Alike(); ArrayXd.DivVA( l, r, res ); RETURN res;
END "DIV";
PROCEDURE "MOD"*( l: Matrix; r: Value ): Matrix;
VAR res: Matrix;
BEGIN
res := l.Alike(); ArrayXd.ModAV( l, r, res ); RETURN res;
END "MOD";
PROCEDURE "MOD"*( l: Value; r: Matrix ): Matrix;
VAR res: Matrix;
BEGIN
res := r.Alike(); ArrayXd.ModVA( l, r, res ); RETURN res;
END "MOD";
PROCEDURE "*"*( l: Vec.Vector; r: Matrix ): Vec.Vector;
VAR res: Vec.Vector; rc, rr, lc, i, j: LONGINT; sum: Value;
BEGIN
rc := r.cols; rr := r.rows; lc := l.lenx;
IF lc # rr THEN DataErrors.Error( "The supplied matrix / vector were incompatible." ); HALT( 100 );
ELSE
NEW( res, 0, rc );
FOR i := 0 TO rc - 1 DO
sum := 0;
FOR j := 0 TO lc - 1 DO sum := sum + l.Get( j ) * r.Get( j, i ); END;
res.Set( i, sum );
END;
END;
RETURN res;
END "*";
PROCEDURE "*"*( l: Matrix; r: Vec.Vector ): Vec.Vector;
VAR res: Vec.Vector; lr, lc, rr, i, j: LONGINT; sum: Value;
BEGIN
lr := l.rows; lc := l.cols; rr := r.lenx;
IF lc # rr THEN DataErrors.Error( "The supplied matrix / vector were incompatible." ); HALT( 100 );
ELSE
NEW( res, 0, lr );
FOR i := 0 TO lr - 1 DO
sum := 0;
FOR j := 0 TO lc - 1 DO sum := sum + l.Get( i, j ) * r.Get( j ); END;
res.Set( i, sum );
END;
END;
RETURN res;
END "*";
PROCEDURE "*"*( l, r: Matrix ): Matrix;
VAR res: Matrix; rr, rc, lr, lc, i, j, k: LONGINT; sum: Value;
BEGIN
rc := r.cols;
rr := r.rows;
lr := l.rows; lc := l.cols;
IF lc # rr THEN DataErrors.Error( "The supplied matrices were incompatible." ); HALT( 100 )
ELSE
NEW( res, 0, lr, 0, rc );
FOR i := 0 TO lr - 1 DO
FOR j := 0 TO rc - 1 DO
sum := 0;
FOR k := 0 TO lc - 1 DO sum := sum + l.Get( i, k ) * r.Get( k, j ); END;
res.Set( i, j, sum );
END;
END;
END;
RETURN res;
END "*";
PROCEDURE LoadMatrix( R: DataIO.Reader; VAR obj: OBJECT );
VAR a: Matrix; version: SHORTINT; ver: NbrInt.Integer;
BEGIN
R.RawSInt( version );
IF version = -1 THEN
obj := NIL
ELSE
IF version = VERSION THEN NEW( a, 0, 0, 0, 0 ); a.Read( R ); obj := a
ELSE
ver := version; DataErrors.IntError( ver, "Alien version number encountered." ); HALT( 1000 )
END
END
END LoadMatrix;
PROCEDURE StoreMatrix( W: DataIO.Writer; obj: OBJECT );
VAR a: Matrix;
BEGIN
IF obj = NIL THEN W.RawSInt( -1 ) ELSE W.RawSInt( VERSION ); a := obj( Matrix ); a.Write( W ) END
END StoreMatrix;
PROCEDURE Register;
VAR a: Matrix;
BEGIN
NEW( a, 0, 0, 0, 0 ); DataIO.PlugIn( a, LoadMatrix, StoreMatrix )
END Register;
PROCEDURE Load*( R: DataIO.Reader; VAR obj: Matrix );
VAR ptr: OBJECT;
BEGIN
R.Object( ptr ); obj := ptr( Matrix )
END Load;
PROCEDURE Store*( W: DataIO.Writer; obj: Matrix );
BEGIN
W.Object( obj )
END Store;
BEGIN
Register
END MtxInt.