112 lines
4.3 KiB
Text
112 lines
4.3 KiB
Text
BEGIN # LU decomposition - translation of the Kotlin sample #
|
|
|
|
OP * = ( [,]REAL this, other )[,]REAL:
|
|
IF 1 LWB this /= 2 LWB other OR 1 UPB this /= 2 UPB other
|
|
THEN print( ( "Cannot multiply matrices: " ) );
|
|
print( ( "[", whole( 1 LWB this, 0 ), ":", whole( 1 UPB this, 0 ) ) );
|
|
print( ( ",", whole( 2 LWB this, 0 ), ":", whole( 2 UPB this, 0 ), "]" ) );
|
|
print( ( " X " ) );
|
|
print( ( "[", whole( 1 LWB other, 0 ), ":", whole( 1 UPB other, 0 ) ) );
|
|
print( ( ",", whole( 2 LWB other, 0 ), ":", whole( 2 UPB other, 0 ), "]" ) );
|
|
print( ( newline ) );
|
|
stop
|
|
ELSE [ 1 LWB this : 1 UPB this, 2 LWB other : 2 UPB other ]REAL result;
|
|
FOR i FROM 1 LWB this TO 1 UPB this DO
|
|
FOR j FROM 2 LWB other TO 2 UPB other DO
|
|
result[ i, j ] := 0;
|
|
FOR k FROM 1 LWB other TO 1 UPB other DO
|
|
result[ i, j ] +:= this[ i, k ] * other[ k, j ]
|
|
OD
|
|
OD
|
|
OD;
|
|
result
|
|
FI # * # ;
|
|
|
|
OP PIVOTISE = ( [,]REAL m in )[,]REAL:
|
|
BEGIN
|
|
[,]REAL m = m in[ AT 1, AT 1 ];
|
|
INT n = ( 1 UPB m - 1 LWB m ) + 1;
|
|
[ 1 : n, 1 : n ]REAL im;
|
|
FOR i TO n DO
|
|
FOR j TO n DO im[ i, j ] := 0 OD;
|
|
im[ i, i ] := 1
|
|
OD;
|
|
FOR i TO n DO
|
|
REAL max := ABS m[ i, i ];
|
|
INT row := i;
|
|
FOR j FROM i TO n DO
|
|
IF ABS m[ j, i ] > max THEN
|
|
max := ABS m[ j, i ];
|
|
row := j
|
|
FI
|
|
OD;
|
|
IF i /= row THEN
|
|
[]REAL t = im[ i, : ];
|
|
im[ i, : ] := im[ row, : ];
|
|
im[ row, : ] := t
|
|
FI
|
|
OD;
|
|
im
|
|
END # PIVOTISE # ;
|
|
|
|
MODE DECOMPOSITION = STRUCT( [ 1 : 1, 1 : 1 ]REAL l, u, p );
|
|
|
|
OP LU = ( [,]REAL a in )DECOMPOSITION:
|
|
BEGIN
|
|
[,]REAL a = a in[ AT 1, AT 1 ];
|
|
INT n = ( 1 UPB a in - 1 LWB a in ) + 1;
|
|
[ 1 : n, 1 : n ]REAL l, u;
|
|
[,]REAL p = PIVOTISE a;
|
|
[,]REAL a2 = p * a;
|
|
FOR i TO n DO FOR j TO n DO l[ i, j ] := u[ i, j ] := 0 OD OD;
|
|
FOR j TO n DO
|
|
l[ j, j ] := 1;
|
|
FOR i TO j DO
|
|
REAL sum := 0;
|
|
FOR k TO i DO sum +:= u[ k, j ] * l[ i, k ] OD;
|
|
u[ i, j ] := a2[ i, j ] - sum
|
|
OD;
|
|
FOR i FROM j + 1 TO n DO
|
|
REAL sum := 0;
|
|
FOR k TO j DO sum +:= u[ k, j ] * l[ i, k ] OD;
|
|
l[ i, j ] := ( a2[ i, j ] - sum ) / u[ j, j ]
|
|
OD
|
|
OD;
|
|
DECOMPOSITION( l, u, p )
|
|
END # LU # ;
|
|
|
|
PROC show matrix = ( STRING title, [,]REAL m, INT fw, fdp )VOID:
|
|
BEGIN
|
|
print( ( " ", title, newline ) );
|
|
FOR i FROM 1 LWB m TO 1 UPB m DO
|
|
print( ( " " ) );
|
|
FOR j FROM 2 LWB m TO 2 UPB m DO print( ( fixed( m[ i, j ], -fw, fdp ), " " ) ) OD;
|
|
print( ( newline ) )
|
|
OD
|
|
END # show matrix # ;
|
|
|
|
BEGIN # test the LU decomposition operator #
|
|
[,]REAL a1 = ( ( 1.0, 3.0, 5.0 )
|
|
, ( 2.0, 4.0, 7.0 )
|
|
, ( 1.0, 1.0, 0.0 )
|
|
);
|
|
DECOMPOSITION d1 = LU a1;
|
|
print( ( "EXAMPLE 1:-", newline ) );
|
|
show matrix( "A:", a1, 1, 0 );
|
|
show matrix( "L:", l OF d1, 8, 5 );
|
|
show matrix( "U:", u OF d1, 8, 5 );
|
|
show matrix( "P:", p OF d1, 1, 0 );
|
|
print( ( newline ) );
|
|
[,]REAL a2 = ( ( 11.0, 9.0, 24.0, 2.0 )
|
|
, ( 1.0, 5.0, 2.0, 6.0 )
|
|
, ( 3.0, 17.0, 18.0, 1.0 )
|
|
, ( 2.0, 5.0, 7.0, 1.0 )
|
|
);
|
|
DECOMPOSITION d2 = LU a2;
|
|
print( ( "EXAMPLE 2:-", newline ) );
|
|
show matrix( "A:", a2, 2, 0 );
|
|
show matrix( "L:", l OF d2, 8, 5 );
|
|
show matrix( "U:", u OF d2, 8, 5 );
|
|
show matrix( "P:", p OF d2, 1, 0 )
|
|
END
|
|
END
|