Commit e7e7538f authored by Mathieu Hemery's avatar Mathieu Hemery
Browse files

Reflow and comment

parent 53fcbbf0
:- module(matrix,
[
constant_list/3,
constant_matrix/4,
replace/4,
get_element/3,
replace_element/3,
init_matrix/1,
write_matrix/1,
ludcmp_matrix/2,
ludcmp_matrix/0,
matrix_sum/4,
matrix_prod/3,
scalar_prod/3,
matrix_apply/3,
matrix_sum/3,
matrix_sub/3,
create_diag_matrix/3,
sub_matrix/1,
row_sum/3,
transpose_matrix/2,
lusubst/2
[
constant_list/3,
constant_matrix/4,
replace/4,
get_element/3,
replace_element/3,
init_matrix/1,
write_matrix/1,
ludcmp_matrix/2,
ludcmp_matrix/0,
matrix_sum/4,
matrix_prod/3,
scalar_prod/3,
matrix_apply/3,
matrix_sum/3,
matrix_sub/3,
create_diag_matrix/3,
sub_matrix/1,
row_sum/3,
transpose_matrix/2,
lusubst/2
]
).
/* Module Matrix
*
* Several operations are done on a "global variable" of SWI-prolog (using nb_set/getval)
* but this is due to historical version of an efficient implementation in GNU-prolog and
* may not be necessary nor useful now.
*/
/* Tools to create a list or a matrix with every elements/coefficients be the same */
/* constant_list(+Length,+Element,-List) */
%! constant_list(+Length, +Element, -List)
%
% create a list with every element set up to Element
constant_list(0,_,[]) :-
!.
!.
constant_list(N1,X,[X|L]) :-
N1 > 0, N is N1 - 1,
constant_list(N,X,L).
%! constant_matrix(+N_row, +N_column, +Element, -List)
%
% create a matrix with every element set up to Element
constant_matrix(RowSize,ColumnSize,Value,Matrix) :-
constant_list(ColumnSize,Value,Row), % The row size is the number of column
constant_list(RowSize,Row,Matrix).
constant_list(ColumnSize,Value,Row), % The row size is the number of column
constant_list(RowSize,Row,Matrix).
/* Allows to replace the i-th element of the list, where the 1st element is index by 1 */
/* replace(+List,+Index,+NewElement,-NewList) */
%! replace(+List, +Index, +NewElement, -NewList)
%
% Allows to replace the i-th element of the list, where the 1st element is index by 1
replace([_|Tail], 1, NewElement, [NewElement|Tail]) :-
!.
!.
replace([Element|Tail], I, NewElement, [Element|Tail2]):- I > 1, I1 is I-1, replace(Tail, I1, NewElement, Tail2).
replace([Element|Tail], I, NewElement, [Element|Tail2]):-
I > 1, I1 is I-1,
replace(Tail, I1, NewElement, Tail2).
%! init_matrix(+Matrix)
%
% Creates the global variable which will stores the matrix on which we are doing operations
init_matrix([Head|Tail]):-
length([Head|Tail],RowSize),
length(Head,ColumnSize),
nb_setval(stored_matrix,[Head|Tail]),
nb_setval(formal_m,RowSize),
nb_setval(formal_n,ColumnSize).
/* Get the (i,j)-th element of the stored matrix */
%! get_element(+RowIndex, +ColumnIndex, -Value)
%
% get the (i,j)-th element of the stored_matrix
get_element(RowIndex,ColumnIndex,Value) :-
nb_getval(stored_matrix,Matrix),
nth1(RowIndex,Matrix,Row),
nth1(ColumnIndex,Row,Value).
nb_getval(stored_matrix,Matrix),
nth1(RowIndex,Matrix,Row),
nth1(ColumnIndex,Row,Value).
/* Replace the the (i,j)-th element of the stored matrix by Value, and stored the result in the global variable */
replace_element(RowIndex,ColumnIndex,Value) :-
nb_getval(stored_matrix,Matrix),
nth1(RowIndex,Matrix,Row),
replace(Row,ColumnIndex,Value,NewRow),
replace(Matrix,RowIndex,NewRow,NewMatrix),
nb_setval(stored_matrix,NewMatrix).
%! replace_element(+RowIndex, +ColumnIndex, +Value) :-
%
% set the (i,j)-th element of the stored matrix to Value
/* Creates the global variable which will stores the matrix on which we are doing operations */
/* init_matrix(+Matrix) */
replace_element(RowIndex,ColumnIndex,Value) :-
nb_getval(stored_matrix,Matrix),
nth1(RowIndex,Matrix,Row),
replace(Row,ColumnIndex,Value,NewRow),
replace(Matrix,RowIndex,NewRow,NewMatrix),
nb_setval(stored_matrix,NewMatrix).
init_matrix([Head|Tail]):-
length([Head|Tail],RowSize),
length(Head,ColumnSize),
nb_setval(stored_matrix,[Head|Tail]),
nb_setval(formal_m,RowSize),
nb_setval(formal_n,ColumnSize).
/* Some tools to write matrix in a pretty way, for debugging for examples*/
%! write_matrix(+Matrix)
write_matrix([]).
write_matrix([H|T]):-
write_tab(H),
write_matrix(T).
write_tab(H),
write_matrix(T).
%! write_tab(+List)
write_tab([]) :-
nl.
write_tab([]):-nl.
write_tab([H|T]):-
write(H),write('\t'),
write_tab(T).
write(H),write('\t'),
write_tab(T).
/* Decompose in place the matrix into two matrices,lower diagonal and upper diagonal (the diagonal of the lower one is set to unity)
The transformed matrix is directly given by the second argument */
/* ludcmp_matrix(+Matrix,-TransformedMatrix) */
%! ludcmp_matrix(+Matrix, -TransformedMatrix)
%
% Perform an LU decomposition, that is find two matrices L and U such that Matrix = L.U,
% L and U begin respectively lower and upper diagonal.
% By convention, L has only 1 on its diagonal so that we gather them in a single matrix M
% such that M[i,j] = L[i,j] if i<j and M[i,j] = U[i,j] if i>=j. M is the transformed
% matrix and is stored in the global variable at the end of the call.
ludcmp_matrix(A,B):-
init_matrix(A),
ludcmp_matrix,
nb_getval(stored_matrix,B).
init_matrix(A),
ludcmp_matrix,
nb_getval(stored_matrix,B).
/* LU-decomposition directly on the stored matrix, and the global variable is updated to this new matrix */
ludcmp_matrix:-
ludcmp_matrix :-
% Value if pivot element is zero (singular)
TINY=1.0E-100,
nb_getval(formal_m,M),
......@@ -217,62 +252,86 @@ ludcmp_matrix:-
)
).
/* For one element (i,j) of the stored matrix, gives back the sum for K between 0 and N-1 of M(I,K)M(K,J) */
%! matrix_sum(+N, +I, +J, -Sum)
%
% if M denotes the stored matrix, return the sum for K between 1 and N of M(I,K)M(K,J)
matrix_sum(SumSize,RowIndex,ColumnIndex,Sum):-
get_element(RowIndex,ColumnIndex,A),
matrix_sum_rec(1,SumSize,RowIndex,ColumnIndex,A,Sum).
get_element(RowIndex,ColumnIndex,A),
matrix_sum_rec(1,SumSize,RowIndex,ColumnIndex,A,Sum).
matrix_sum_rec(SumSize,SumSize,_,_,Sum,Sum):-
!.
!.
matrix_sum_rec(SumIndex,SumSize,RowIndex,ColumnIndex,PartialSum,Sum):-
get_element(RowIndex,SumIndex,IK),
get_element(SumIndex,ColumnIndex,KJ),
PartialSum2 is PartialSum - IK*KJ,
NewSumIndex is SumIndex+1,
matrix_sum_rec(NewSumIndex,SumSize,RowIndex,ColumnIndex,PartialSum2,Sum).
get_element(RowIndex,SumIndex,IK),
get_element(SumIndex,ColumnIndex,KJ),
PartialSum2 is PartialSum - IK*KJ,
NewSumIndex is SumIndex+1,
matrix_sum_rec(NewSumIndex,SumSize,RowIndex,ColumnIndex,PartialSum2,Sum).
/* Product of two matrixes */
%! matrix_prod(+A, +B, -P)
%
% P = A.B
matrix_prod(MatrixA,MatrixB,Result):-
transpose_matrix(MatrixB,TransposedMatrixB),
matrix_prod_rec(MatrixA,TransposedMatrixB,Result).
transpose_matrix(MatrixB,TransposedMatrixB),
matrix_prod_rec(MatrixA,TransposedMatrixB,Result).
matrix_prod_rec([],_,[]).
matrix_prod_rec([RowA|MatrixA],MatrixB,[RowResult|MatrixResult]):-
matrix_prod_sub(RowA,MatrixB,RowResult),
matrix_prod_rec(MatrixA,MatrixB,MatrixResult).
matrix_prod_sub(RowA,MatrixB,RowResult),
matrix_prod_rec(MatrixA,MatrixB,MatrixResult).
matrix_prod_sub(_,[],[]).
matrix_prod_sub(Row,[RowB|MatrixB],[ElementResult|RowResult]):-
scalar_prod(Row,RowB,ElementResult),
matrix_prod_sub(Row,MatrixB,RowResult).
scalar_prod(Row,RowB,ElementResult),
matrix_prod_sub(Row,MatrixB,RowResult).
%! scalar_prod(+X, +Y, -Product)
%
% Comput the scalar product of X and Y
scalar_prod(VecX, VecY, SP) :-
scalar_prod(VecX, VecY, 0, SP).
scalar_prod([],[],0).
scalar_prod([ElementX|VectorX],[ElementY|VectorY],Result):-
scalar_prod(VectorX,VectorY,PartialResult),
Result is PartialResult + ElementX*ElementY.
scalar_prod([], [], SP, SP).
/* Apply matrix A to vector B */
scalar_prod([X|TailX], [Y|TailY], Res_tempo, Acc):-
Res is Res_tempo + X*Y.
scalar_prod(TailX, TailY, Res, Acc).
%! matrix_apply(+M, +V, -Res)
%
% compute the dot product of M.V
matrix_apply([],_,[]).
matrix_apply([Row|Matrix],Vector,[ResultElement|ResultVector]):-
scalar_prod(Row,Vector,ResultElement),
matrix_apply(Matrix,Vector,ResultVector).
scalar_prod(Row,Vector,ResultElement),
matrix_apply(Matrix,Vector,ResultVector).
/* Sum of two matrixes */
%! matrix_sum(+A, +B, -Sum)
%
% Compute the element sum of two matrices
matrix_sum([],[],[]).
matrix_sum([RowA|MatrixA],[RowB|MatrixB],[SumRow|SumMatrix]):-
row_sum(RowA,RowB,SumRow),
matrix_sum(MatrixA,MatrixB,SumMatrix).
row_sum(RowA,RowB,SumRow),
matrix_sum(MatrixA,MatrixB,SumMatrix).
row_sum([],[],[]).
row_sum([ElementA|RowA],[ElementB|Row2],[SumElement|SumRow]):-
SumElement is ElementA+ElementB,
row_sum(RowA,Row2,SumRow).
SumElement is ElementA+ElementB,
row_sum(RowA,Row2,SumRow).
/* Sub of two matrixes */
......@@ -286,67 +345,78 @@ row_sub([ElementA|RowA],[ElementB|Row2],[SubElement|SubRow]):-
SubElement is ElementA-ElementB,
row_sub(RowA,Row2,SubRow).
/* Create a diagonal matrix and stores it in the global variable */
%! create_diag_matrix(+Size, +Value, -Diag_matrix)
%
% Store a diagonal matrix in the global variable and return it
create_diag_matrix(Size,Value,DiagMatrix):-
constant_matrix(Size,Size,0,Matrix),
nb_setval(stored_matrix,Matrix),
nb_setval(formal_n,Size),
nb_setval(formal_m,Size),
forall(
between(1,Size,RowIndex),
replace_element(RowIndex,RowIndex,Value)
),
nb_getval(stored_matrix,DiagMatrix).
constant_matrix(Size,Size,0,Matrix),
nb_setval(stored_matrix,Matrix),
nb_setval(formal_n,Size),
nb_setval(formal_m,Size),
forall(
between(1,Size,RowIndex),
replace_element(RowIndex,RowIndex,Value)
),
nb_getval(stored_matrix,DiagMatrix).
/* Substract one matrix from the stored matrix */
% sub_matrix(+Matrix)
%
% Substract Matrix from the global matrix
sub_matrix([Row|Matrix]):-
sub_matrix([Row|Matrix],1).
sub_matrix([Row|Matrix],1).
sub_matrix([],_).
sub_matrix([Row|Matrix],RowIndex):-
sub_matrix_rec(Row,RowIndex,1),
NewRowIndex is RowIndex+1,
sub_matrix(Matrix,NewRowIndex).
sub_matrix_rec(Row,RowIndex,1),
NewRowIndex is RowIndex+1,
sub_matrix(Matrix,NewRowIndex).
sub_matrix_rec([],_,_).
sub_matrix_rec([Element|Row],RowIndex,ColumnIndex):-
get_element(RowIndex,ColumnIndex,IJ),
NewElement is IJ-Element,
replace_element(RowIndex,ColumnIndex,NewElement),!,
NewColumnIndex is ColumnIndex+1,
sub_matrix_rec(Row,RowIndex,NewColumnIndex).
get_element(RowIndex,ColumnIndex,IJ),
NewElement is IJ-Element,
replace_element(RowIndex,ColumnIndex,NewElement),!,
NewColumnIndex is ColumnIndex+1,
sub_matrix_rec(Row,RowIndex,NewColumnIndex).
/* Transpose matrix and stores it in the global variable */
%! transpose_matrix(+Matrix, -Transposed)
%
% Transpose Matrix and store it in the global variable (which is now the transposed
% matrix)
transpose_matrix(Matrix,TransposedMatrix):-
transpose_matrix(Matrix),
nb_getval(stored_matrix,TransposedMatrix).
transpose_matrix(Matrix),
nb_getval(stored_matrix,TransposedMatrix).
transpose_matrix([Head|Tail]):-
length([Head|Tail],ColumnSize), %These are column and row size of the TRANSPOSED matrix
length(Head,RowSize),
constant_matrix(RowSize,ColumnSize,0,Matrix),!,
nb_setval(stored_matrix,Matrix),
nb_setval(formal_m,RowSize),
nb_setval(formal_n,ColumnSize),
fill_matrix_transposed([Head|Tail],1).
length([Head|Tail],ColumnSize), %These are column and row size of the TRANSPOSED matrix
length(Head,RowSize),
constant_matrix(RowSize,ColumnSize,0,Matrix),!,
nb_setval(stored_matrix,Matrix),
nb_setval(formal_m,RowSize),
nb_setval(formal_n,ColumnSize),
fill_matrix_transposed([Head|Tail],1).
fill_matrix_transposed([],_).
fill_matrix_transposed([Head|Tail],RowIndex):-
fill_matrix_transposed_rec(Head,RowIndex,1),
NewRowIndex is RowIndex+1,
fill_matrix_transposed(Tail,NewRowIndex).
fill_matrix_transposed_rec(Head,RowIndex,1),
NewRowIndex is RowIndex+1,
fill_matrix_transposed(Tail,NewRowIndex).
fill_matrix_transposed_rec([],_,_).
fill_matrix_transposed_rec([Head|Tail],RowIndex,ColumnIndex):-
replace_element(ColumnIndex,RowIndex,Head),!,
NewColumnIndex is ColumnIndex+1,
fill_matrix_transposed_rec(Tail,RowIndex,NewColumnIndex).
replace_element(ColumnIndex,RowIndex,Head),!,
NewColumnIndex is ColumnIndex+1,
fill_matrix_transposed_rec(Tail,RowIndex,NewColumnIndex).
/* Use LU decomposition to solve Ax=b by fwd/backwd substitution
Supposing A has already been decomposed using the ludcmp predicate */
%! lusubst(+B, -X)
%
% Use LU decomposition to solve Ax=b by fwd/backwd substitution
% Supposing A has already been decomposed using the ludcmp predicate
lusubst(B,X):-
nb_getval(formal_n,N),
......@@ -398,4 +468,4 @@ lusubst(B,X):-
replace(Vector,II,YY,NewVector),nb_setval(formal_vector,NewVector)
)
),
nb_getval(formal_vector,X).
\ No newline at end of file
nb_getval(formal_vector,X).
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment