1:- module(ffimatrix,
    2          [ matrix_write/1,
    3            matrix_new/3,
    4            matrix_new/4,
    5            matrix_mul/3,
    6            matrix_dims/2,
    7            matrix_size/2,
    8            matrix_type/2,
    9            matrix_to_list/2,
   10            matrix_get/3,
   11            matrix_set/3,
   12            matrix_set_all/2,
   13            matrix_add/3,
   14            matrix_inc/2,
   15            matrix_inc/3,
   16            matrix_dec/2,
   17            matrix_dec/3,
   18            matrix_max/2,
   19            matrix_min/2,
   20            matrix_sum/2,
   21            matrix_transpose/2,
   22            matrix_maxarg/2,
   23            matrix_minarg/2,
   24            matrix_agg_lines/2,
   25            matrix_agg_cols/2,
   26            matrix_op/4,
   27            matrix_op_to_all/4,
   28            matrix_select/4,
   29            matrix_createZeroes/3,
   30            matrixSetZeroes/3,
   31            matrixSetOnes/3,
   32            matrix_read/2,
   33            matrix_setAllNumbers/3,
   34            matrix_eye/3,
   35            matrix_equal/2,
   36            matrix_map/3,
   37            matrix_foreach/2,
   38            matrix_foldl/3,
   39            matrix_from_list/2,
   40            cholesky_function/1,
   41            matrix_random/3,
   42            matrix_to_list_of_lists/2
   43          ]).   44
   45:- use_module(library(ffi)).   46:- use_module(library(lambda)).   47
   48:- multifile
   49    ffi:library_path_hook/3,                    % +Spec, -Path, +Options
   50    ffi:cpp_hook/2.   51
   52
   53ffi:cpp_hook(path(clang),[ '-I/usr/local/opt/openblas/include','-E', '-xc',-]).
   54ffi:library_path_hook(liblapack,'/usr/local/opt/openblas/lib/liblapack.dylib',_).
   55ffi:library_path_hook(libblas,'/usr/local/opt/openblas/lib/libblas.dylib',_).
   56
   57:- c_import("#include <lapacke.h>",
   58            [liblapack],
   59            ['LAPACKE_dpotrf'(int, char,int, *double, int)]).
   60
   61cholesky(matrix(doubles, [Order, Order], A)) :-
   62    'LAPACKE_dpotrf'( 101, 'L',Order,A, 3).
 cholesky_function
Given in input M1 make the decomposition of Cholesky into M2
   66cholesky_function(Matrix1) :-
   67    Matrix1=matrix(Type, [Order, Order], _),
   68    matrix_new(Type, [Order, Order], Matrix2),
   69    matrix_sollower(Matrix1, Matrix2),
   70    cholesky(Matrix2),
   71    matrix_write(Matrix2).
   72
   73matrix_sollower(matrix(Type, [NRows, NColumns], Elements1), matrix(Type, [NRows, NColumns], Elements2)) :-
   74    sollowerMatrix(NRows, NColumns, Elements1, Elements2).
   75
   76
   77matrix_copy(matrix(Type, [NRows, NColumns], Elements1), matrix(Type, [NRows, NColumns], Elements2)) :-
   78    NElements is NRows*NColumns,
   79    arrayCopy(NElements, Elements1, Elements2).
   80
   81matrix_to_list_of_lists(M, L) :-
   82    aux0(M, 0, 0, L).
   83
   84aux0(matrix(_, [NRows, _], _), NRows, _, []) :- !.      
   85aux0(M, Row, NColumns, [[]|MoreRows]) :-
   86    M=matrix(_, [_, NColumns], _), !,
   87    Row1 is Row+1,
   88    aux0(M, Row1, 0, MoreRows).
   89aux0(M, Row, Column, [[Element|MoreElements]|MoreRows]) :-
   90    matrix_get(M, [Row, Column], Element),
   91    Column1 is Column+1,
   92    aux0(M, Row, Column1, [MoreElements|MoreRows]).
   93
   94
   95% https://www.dcc.fc.up.pt/~vsc/Yap/documentation.html#matrix
   96:- c_import("#include <cblas.h>",
   97            [libblas],
   98            
   99            [ cblas_dgemm(int,
  100                          int,
  101                          int,
  102                          int,
  103                          int,
  104                          int,
  105                          double,
  106                          *double,
  107                          int,
  108                          *double,
  109                          int,
  110                          double,
  111                          *double,
  112                          int)
  113            ]).
  114
  115:- (multifile user:term_expansion/2).  116:- (dynamic user:term_expansion/2).  117                
  118
  119
  120:-c_import("#include \"../matrixNative.h\"", 
  121    ['../lib/x86_64-darwin/matrixNative'], 
  122    [matrixSetAll(int, int, *double, double),      matrixEye(int, int, *double),      matrixCopy(int, int, *double, *double),      arrayCopy(int, *double, *double),      sollowerMatrix(int, int, *double, *double)]).
 matrix_write(+M1)
Print elements of a matrix
  130matrix_write(matrix(_, [NRows, NColumns], Matrix)) :-
  131    write_matrix(Matrix, NRows, NColumns).
  132
  133write_matrix(Matrix, NRows, NColumns) :-
  134    write_matrix(Matrix, 0, 0, NRows, NColumns).
  135
  136write_matrix(_, NRows, _, NRows, _) :- !.
  137
  138write_matrix(MatrixElements, CurrentRow, NColumns, NRows, NColumns) :- !,
  139    writeln(''),
  140    CurrentRow1 is CurrentRow+1,
  141    write_matrix(MatrixElements, CurrentRow1, 0, NRows, NColumns).
  142
  143write_matrix(MatrixElements, CurrentRow, CurrentColumn, NRows, NColumns) :- !,
  144    Index is CurrentRow*NColumns+CurrentColumn,
  145    c_load(MatrixElements[Index], Element),
  146    write(Element),
  147    write(' '),
  148    CurrentColumn1 is CurrentColumn+1,
  149    write_matrix(MatrixElements, CurrentRow, CurrentColumn1, NRows, NColumns).
  150
  151%!  matrix_new(+Type,+Dimension,-M1)
  152% Create a new matrix. In the matrix it must be specified Type (ints,floats,doubles) and the dimension with Nrows and Ncolumns. The matrix will not be initialised.
  153matrix_new(doubles, [NRows, NColumns], matrix(doubles, [NRows, NColumns], Matrix)) :-
  154    NRows>0,
  155    NColumns>0,
  156    Nelements is NRows*NColumns,
  157    c_alloc(Matrix, double[Nelements]).
  158
  159%!  matrix_new(+Type,+Dimension,+Elements,-M1)
  160% Create a new matrix. In the matrix it must be specified Type (ints,floats,doubles) and the dimension with NRows and NColumns. The matrix will be initialised with elements specified in a list.
  161matrix_new(doubles, [NRows, NColumns], ListElements, matrix(doubles, [NRows, NColumns], Matrix)) :-
  162    NRows>0,
  163    NColumns>0,
  164    NElements is NRows*NColumns,
  165    length(ListElements, NElements),
  166    c_alloc(Matrix, double[]=ListElements).
 matrix_mul(+M1, +M2, -M3)
Multiplication between two matrices using "cblas_dgemm" function from "cblas.h"
  170matrix_mul(matrix(Type, [NRowsMatrix1, NColumnsMatrix1], Matrix1), matrix(Type, [NColumnsMatrix1, NColumnsMatrix2], Matrix2), matrix(Type, [NRowsMatrix1, NColumnsMatrix2], Matrix3)) :-
  171    matrix_new(Type,
  172               [NRowsMatrix1, NColumnsMatrix2],
  173               matrix(Type, [NRowsMatrix1, NColumnsMatrix2], Matrix3)),
  174    cblas_dgemm(102,
  175                111,
  176                111,
  177                NRowsMatrix1,
  178                NColumnsMatrix2,
  179                NColumnsMatrix1,
  180                1,
  181                Matrix1,
  182                NRowsMatrix1,
  183                Matrix2,
  184                NColumnsMatrix1,
  185                0,
  186                Matrix3,
  187                NRowsMatrix1).
 matrix_dims(+M1, -Dimension)
Calculates the total number of elements in the matrix. The value is contained in Dimension.
  191matrix_dims(matrix(_, Dimension, _), Dimension).
 matrix_size(+M1, -Nelements)
Calculates the number of elements in the matrix. We have to pass the number of the rows and the number of the column.
  195matrix_size(matrix(_, [NRows, NColumns], _), Nelements) :-
  196    Nelements is NRows*NColumns.
 matrix_type(+M1, -Type)
Function that return the Type of the matrix passed.
  200matrix_type(matrix(Type, _, _), Type).
 matrix_to_list(+M1, -List)
Function that transforms the matrix passed in input in a list.
  204to_list(_, NRows, _, NRows, _, List) :- !,
  205    List=([]).
  206
  207to_list(MatrixElements, CurrentRow, NColumns, Nrows, NColumns, List) :- !,
  208    CurrentRow1 is CurrentRow+1,
  209    to_list(MatrixElements, CurrentRow1, 0, Nrows, NColumns, List).
  210
  211to_list(MatrixElements, CurrentRow, CurrentColumn, NRows, NColumns, [Element|MoreElements]) :-
  212    Index is CurrentRow*NColumns+CurrentColumn,
  213    c_load(MatrixElements[Index], Element),
  214    CurrentColumn1 is CurrentColumn+1,
  215    to_list(MatrixElements, CurrentRow, CurrentColumn1, NRows, NColumns, MoreElements).
  216
  217matrix_to_list(matrix(_, [NRows, NColumns], MatrixElements), List) :-
  218    to_list(MatrixElements, 0, 0, NRows, NColumns, List).
  219
  220%!  matrix_get(+M1,+Position,-Element)
  221% matrix_get(M,[Row,Column],E) unifies E with the element of M at position [Row, Column] 
  222matrix_get(matrix(_, [_, NColumns], MatrixElements), [Row, Column], Element) :-
  223    Index is Row*NColumns+Column,
  224    c_load(MatrixElements[Index], Element).
  225
  226%!  matrix_set(+Type,+Position,-Element)
  227% matrix_set(M,[Row,Column],E) update the element of M at position [Row, Column] with the element specified in Element.
  228matrix_set(matrix(_, [_, NColumns], MatrixElements), [Row, Column], Element) :-
  229    Index is Row*NColumns+Column,
  230    c_store(MatrixElements[Index], Element).
 matrix_map(+Predicate, +M1, -M2)
Function that return true if Predicate can successfully be applied on all elements of M1.
  234matrix_map(Predicate, matrix(Type, [NRows, NColumns], MatrixElements1), matrix(Type, [NRows, NColumns], MatrixElements2)) :-
  235    matrix_new(Type,
  236               [NRows, NColumns],
  237               matrix(_, _, MatrixElements2)),
  238    matrix_map(Predicate, NRows, NColumns, MatrixElements1, MatrixElements2, 0, 0).
  239
  240matrix_map(_, NRows, _, _, _, NRows, _) :- !.
  241
  242matrix_map(Predicate, NRows, NColumns, MatrixElements1, MatrixElements2, CurrentRow, NColumns) :- !,
  243    CurrentRow1 is CurrentRow+1,
  244    matrix_map(Predicate, NRows, NColumns, MatrixElements1, MatrixElements2, CurrentRow1, 0).
  245
  246matrix_map(Predicate, NRows, NColumns, MatrixElements1, MatrixElements2, CurrentRow, CurrentColumn) :-
  247    Index is CurrentRow*NColumns+CurrentColumn,
  248    c_load(MatrixElements1[Index], E1),
  249    call(Predicate, E1, E2),
  250    c_store(MatrixElements2[Index], E2),
  251    CurrentColumn1 is CurrentColumn+1,
  252    matrix_map(Predicate,
  253               NRows,
  254               NColumns,
  255               MatrixElements1,
  256               MatrixElements2,
  257               CurrentRow,
  258               CurrentColumn1).
  259
  260%%%%%%%%%%%%%%%%%%%%%%%%%%%%       TWO FUNCTIONS USED BY MATRIX FOREACH       %%%%%%%%%%%%%%%%%%%%%%%%%%%%
  261
  262% Unify the element with the number specified.
  263myset(_, Number, Number).
  264
  265myset(_, 0).
  266% This function find the max number between two number.
  267findMax(FirstNumber, SecondNumber, FirstNumber) :-
  268    FirstNumber>=SecondNumber, !.
  269
  270findMax(_, SecondNumber, SecondNumber).
  271
  272%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 matrix_foreach(+M1, +Predicate)
For each element E of M1, enumerated by rows, calls Predicate(+E,-F) and replaces E with F.
  276matrix_foreach(matrix(_, [NRows, NColumns], MatrixElements), Predicate) :-
  277    matrix_foreach(Predicate, NRows, NColumns, MatrixElements, 0, 0).
  278
  279matrix_foreach(_, NRows, _, _, NRows, _) :- !.
  280
  281matrix_foreach(Predicate, NRows, NColumns, MatrixElements, CurrentRow, NColumns) :- !,
  282    CurrentRow1 is CurrentRow+1,
  283    matrix_foreach(Predicate, NRows, NColumns, MatrixElements, CurrentRow1, 0).
  284
  285matrix_foreach(Predicate, NRows, NColumns, MatrixElements, CurrentRow, CurrentColumn) :-
  286    Index is CurrentRow*NColumns+CurrentColumn,
  287    c_load(MatrixElements[Index], E),
  288    call(Predicate, E, E1),
  289    c_store(MatrixElements[Index], E1),
  290    CurrentColumn1 is CurrentColumn+1,
  291    matrix_foreach(Predicate,
  292                   NRows,
  293                   NColumns,
  294                   MatrixElements,
  295                   CurrentRow,
  296                   CurrentColumn1).
  297
  298%!  matrix_foldl(+M1,+Predicate,-Max)
  299% Fold a matrix, using arguments of the matrix as left argument.
  300matrix_foldl(matrix(_, [NRows, NColumns], MatrixElements), Predicate, Max) :-
  301    c_load(MatrixElements[0], TmpMax),
  302    matrix_start_foldl(Predicate,
  303                       NRows,
  304                       NColumns,
  305                       MatrixElements,
  306                       TmpMax,
  307                       Max).
  308
  309matrix_start_foldl(_, 1, 1, _, Max, Max) :- !.
  310matrix_start_foldl(Predicate, NRows, NColumns, MatrixElements, Tmp, Max) :-
  311    NColumns>1, !,
  312    matrix_foldl(Predicate,
  313                 NRows,
  314                 NColumns,
  315                 MatrixElements,
  316                 0,
  317                 1,
  318                 Tmp,
  319                 Max).
  320
  321matrix_start_foldl(Predicate, NRows, NColumns, MatrixElements, Tmp, Max) :-
  322    matrix_foldl(Predicate,
  323                 NRows,
  324                 NColumns,
  325                 MatrixElements,
  326                 1,
  327                 0,
  328                 Tmp,
  329                 Max).
  330
  331matrix_foldl(_, NRows, _, _, NRows, _, Max, Max) :- !.
  332
  333matrix_foldl(Predicate, NRows, NColumns, MatrixElements, CurrentRow, NColumns, TmpMax, Max) :- !,
  334    CurrentRow1 is CurrentRow+1,
  335    matrix_foldl(Predicate,
  336                 NRows,
  337                 NColumns,
  338                 MatrixElements,
  339                 CurrentRow1,
  340                 0,
  341                 TmpMax,
  342                 Max).
  343
  344matrix_foldl(Predicate, NRows, NColumns, MatrixElements, CurrentRow, CurrentColumn, TmpMax, Max) :-
  345    Index is CurrentRow*NColumns+CurrentColumn,
  346    c_load(MatrixElements[Index], E),    call(Predicate, TmpMax, E, TmpMax1),    CurrentColumn1 is CurrentColumn+1,    matrix_foldl(Predicate,                 NRows,                 NColumns,                 MatrixElements,                 CurrentRow,                 CurrentColumn1,                 TmpMax1,                 Max).
 matrix_set_all(+Number, +M1)
matrix_set_all set all element of the matrix to Element.
  360set_all(_, NRows, _, _, NRows, _) :- !.
  361
  362set_all(Number, NRows, NColumns, MatrixElements, CurrentRow, NColumns) :- !,
  363    CurrentRow1 is CurrentRow+1,
  364    set_all(Number, NRows, NColumns, MatrixElements, CurrentRow1, 0).
  365
  366set_all(Number, NRows, NColumns, MatrixElements, CurrentRow, CurrentColumn) :-
  367    Index is CurrentRow*NColumns+CurrentColumn,
  368    c_store(MatrixElements[Index], Number),
  369    CurrentColumn1 is CurrentColumn+1,
  370    set_all(Number, NRows, NColumns, MatrixElements, CurrentRow, CurrentColumn1).
  371
  372matrix_set_all(number(Type, Number), matrix(Type, [NRows, NColumns], Element)) :-
  373    set_all(Number, NRows, NColumns, Element, 0, 0).
  374
  375%!  matrix_add(+M1,Position,+Number)
  376% matrix_add(M,[Row,Column],Number) add Number to the element of matrix M in position [Row,Column].
  377matrix_add(matrix(_, [_, NColumns], MatrixElements), [Row, Column], Number) :-
  378    Index is Row*NColumns+Column,
  379    c_load(MatrixElements[Index], Tmp),
  380    NewElement is Tmp+Number,
  381    c_store(MatrixElements[Index], NewElement).
  382
  383%!  matrix_inc(+M1,+Position)
  384% matrix_inc(M,[Row,Column]) add 1 (one) to the element of matrix M in position [Row,Column].
  385matrix_inc(matrix(_, [_, NColumns], MatrixElements), [Row, Column]) :-
  386    Index is Row*NColumns+Column,
  387    c_load(MatrixElements[Index], Tmp),
  388    E is Tmp+1,
  389    c_store(MatrixElements[Index], E).
  390
  391%!  matrix_inc(+M1,+Position,-NewElement)
  392% matrix_inc(M,[Row,Column],NewElement) add 1 (one) to the element of matrix M in position [Row,Column] and unify the new element with NewElement.
  393matrix_inc(matrix(_, [_, NColumns], MatrixElements), [Row, Column], NewElement) :-
  394    Index is Row*NColumns+Column,
  395    c_load(MatrixElements[Index], Tmp),
  396    NewElement is Tmp+1,
  397    c_store(MatrixElements[Index], NewElement).
  398
  399%!  matrix_dec(+M1,+Position)
  400% matrix_dec(M,[Row,Column]) decrement 1 (one) to the element of matrix M in position [Row,Column].
  401matrix_dec(matrix(_, [_, NColumns], MatrixElements), [Row, Column]) :-
  402    Index is Row*NColumns+Column,
  403    c_load(MatrixElements[Index], Tmp),
  404    NewElement is Tmp-1,
  405    c_store(MatrixElements[Index], NewElement).
  406
  407%!  matrix_dec(+M1,+Position,-NewElement)
  408% matrix_dec(M,[Row,Column],NewElement) decrement 1 (one) to the element of matrix M in position [Row,Column] and unify the new element with NewElement.
  409matrix_dec(matrix(_, [_, NColumns], MatrixElements), [Row, Column], NewElement) :-
  410    Index is Row*NColumns+Column,
  411    c_load(MatrixElements[Index], Tmp),    NewElement is Tmp-1,    c_store(MatrixElements[Index], NewElement).
 matrix_max(+M1, -Max)
matrix_max(M,Max) unify Max with the maximum element of the matrix M.
  417max_matrix(NRows, _, _, NRows, _, Max, Max) :- !.
  418
  419max_matrix(NRows, NColumns, MatrixElements, CurrentRow, NColumns, TmpMax, Max) :- !,
  420    CurrentRow1 is CurrentRow+1,
  421    max_matrix(NRows, NColumns, MatrixElements, CurrentRow1, 0, TmpMax, Max).
  422
  423max_matrix(NRows, NColumns, MatrixElements, CurrentRow, CurrentColumn, TmpMax, Max) :-
  424    Index is CurrentRow*NColumns+CurrentColumn,
  425    c_load(MatrixElements[Index], E),
  426    CurrentColumn1 is CurrentColumn+1,
  427    (   E>TmpMax
  428    ->  TmpMax1 is E
  429    ;   TmpMax1 is TmpMax
  430    ),
  431    max_matrix(NRows,
  432               NColumns,
  433               MatrixElements,
  434               CurrentRow,
  435               CurrentColumn1,
  436               TmpMax1,
  437               Max).
  438
  439matrix_max(matrix(_, [NRows, NColumns], MatrixElements), Max) :-
  440    c_load(MatrixElements[0], TmpMax),    max_matrix(NRows, NColumns, MatrixElements, 0, 0, TmpMax, Max).
 matrix_min(+M1, -Min)
matrix_min(M,Min) unify Min with the minimum element of the matrix M.
  445matrix_min(Matrix, Min) :-
  446    matrix_foldl(Matrix,
  447                 \X^Y^Z^(X=<Y->Z=X;Z=Y),
  448                 Min).
 matrix_sum(+M1, -Sum)
matrix_sum(M,Sum) unify Sum with the sum of the elements in M.
  452matrix_sum(Matrix, Sum) :-
  453    matrix_foldl(Matrix,
  454                 \CurrentSum^Element^NewSum^(NewSum is CurrentSum+Element),
  455                 Sum).
 matrix_transpose(+M1, -M2)
matrix_transpose(M1,M2) transpose matrix M1 to M2.
  459transpose_matrix(NRows, _, _, _, NRows, _) :- !.
  460
  461transpose_matrix(NRows, NColumns, MatrixElements1, MatrixElements2, CurrentRow, NColumns) :- !,
  462    CurrentRow1 is CurrentRow+1,
  463    transpose_matrix(NRows, NColumns, MatrixElements1, MatrixElements2, CurrentRow1, 0).
  464
  465transpose_matrix(NRows, NColumns, MatrixElements1, MatrixElements2, CurrentRow, CurrentColumn) :-
  466    Index is CurrentRow*NColumns+CurrentColumn,
  467    IndexMatrix2 is CurrentColumn*NRows+CurrentRow,
  468    c_load(MatrixElements1[Index], E),
  469    c_store(MatrixElements2[IndexMatrix2], E),
  470    CurrentColumn1 is CurrentColumn+1,
  471    transpose_matrix(NRows,
  472                     NColumns,
  473                     MatrixElements1,
  474                     MatrixElements2,
  475                     CurrentRow,
  476                     CurrentColumn1).
  477
  478
  479matrix_transpose(matrix(Type, [NRows, NColumns], MatrixElements1), Matrix2) :-
  480    matrix_new(Type, [NColumns, NRows], Matrix2),
  481    Matrix2=matrix(_, _, MatrixElements2),
  482    transpose_matrix(NRows, NColumns, MatrixElements1, MatrixElements2, 0, 0).
 matrix_maxarg(+M1, -MaxPosition)
matrix_maxarg(M,[NRowMax, NColumnMax]) unify [NRowsMax, NColumnMax] with the position of the maximum in matrix M.
  486maxarg_matrix(NRows, _, _, NRows, _, _, [NRowsMax, NColumnsMax], [NRowsMax, NColumnsMax]) :- !.
  487
  488maxarg_matrix(NRows, NColumns, MatrixElements, CurrentRow, NColumns, TmpMax, [NRowTmp, NColumnTmp], [NRowMax, NColomnMax]) :- !,
  489    CurrentRow1 is CurrentRow+1,
  490    maxarg_matrix(NRows,
  491                  NColumns,
  492                  MatrixElements,
  493                  CurrentRow1,
  494                  0,
  495                  TmpMax,
  496                  [NRowTmp, NColumnTmp],
  497                  [NRowMax, NColomnMax]).
  498
  499                maxarg_matrix(NRows, NColumns, MatrixElements, CurrentRow, CurrentColumn, TmpMax, [NRowTmp, NColumnTmp], [NRowMax, NColomnMax]) :-
  500    Index is CurrentRow*NColumns+CurrentColumn,
  501    c_load(MatrixElements[Index], Element),
  502    CurrentColumn1 is CurrentColumn+1,
  503    (   Element>=TmpMax
  504    ->  ( TmpMax1 is Element,
  505          NRowTmp1 is CurrentRow,
  506          NColumnTmp1 is CurrentColumn
  507        ),
  508        maxarg_matrix(NRows,
  509                      NColumns,
  510                      MatrixElements,
  511                      CurrentRow,
  512                      CurrentColumn1,
  513                      TmpMax1,
  514                      [NRowTmp1, NColumnTmp1],
  515                      [NRowMax, NColomnMax])
  516    ;   TmpMax1 is TmpMax,
  517        NRowTmp1 is NRowTmp,
  518        NColumnTmp1 is NColumnTmp,
  519        maxarg_matrix(NRows,
  520                      NColumns,
  521                      MatrixElements,
  522                      CurrentRow,
  523                      CurrentColumn1,
  524                      TmpMax1,
  525                      [NRowTmp1, NColumnTmp1],
  526                      [NRowMax, NColomnMax])
  527    ).
  528
  529matrix_maxarg(matrix(_, [NRows, NColumns], MatrixElements), [NRowMax, NColomnMax]) :-
  530    c_load(MatrixElements[0], TmpMax),    maxarg_matrix(NRows,                  NColumns,                  MatrixElements,                  0,                  0,                  TmpMax,                  [0, 0],                  [NRowMax, NColomnMax]).
 matrix_minarg(+M1, -MinPosition)
matrix_minarg(M,[NRowMax, NColumnMax]) unify [NRowsMax, NColumnMax] with the position of the minimum in matrix M.
  542minarg_matrix(NRows, _, _, NRows, _, _, [NRowTmp, NColumnTmp], [NRowTmp, NColumnTmp]) :- !.
  543
  544minarg_matrix(NRows, NColumns, MatrixElements, CurrentRow, NColumns, TmpMax, [NRowTmp, NColumnTmp], [NRowMax, NColumnMax]) :- !,
  545    CurrentRow1 is CurrentRow+1,
  546    minarg_matrix(NRows,
  547                  NColumns,
  548                  MatrixElements,
  549                  CurrentRow1,
  550                  0,
  551                  TmpMax,
  552                  [NRowTmp, NColumnTmp],
  553                  [NRowMax, NColumnMax]).
  554
  555                minarg_matrix(NRows, NColumns, MatrixElements, CurrentRow, CurrentColumn, TmpMax, [NRowTmp, NColumnTmp], [NRowMax, NColumnMax]) :-
  556    Index is CurrentRow*NColumns+CurrentColumn,
  557    c_load(MatrixElements[Index], Element),
  558    CurrentColumn1 is CurrentColumn+1,
  559    (   Element=<TmpMax
  560    ->  TmpMax1 is Element,
  561        NRowTmp1 is CurrentRow,
  562        NColumnTmp1 is CurrentColumn
  563    ;   TmpMax1 is TmpMax,
  564        NRowTmp1 is NRowTmp,
  565        NColumnTmp1 is NColumnTmp
  566    ),
  567    minarg_matrix(NRows,
  568                  NColumns,
  569                  MatrixElements,
  570                  CurrentRow,
  571                  CurrentColumn1,
  572                  TmpMax1,
  573                  [NRowTmp1, NColumnTmp1],
  574                  [NRowMax, NColumnMax]).
  575
  576matrix_minarg(matrix(_, [NRows, NColumns], MatrixElements), [NRowMax, NColumnMax]) :-
  577    c_load(MatrixElements[0], TmpMax),    minarg_matrix(NRows,                  NColumns,                  MatrixElements,                  0,                  0,                  TmpMax,                  [0, 0],                  [NRowMax, NColumnMax]).
 matrix_agg_lines(+M1, +AggregationMatrix)
matrix_agg_lines (M,AggregationMatrix) sum all the elements of every row of M and the new value will be insert in the AggregationMatrix.
  589agg_lines_matrix(NRows, _, _, _, NRows, _, _) :- !.
  590
  591agg_lines_matrix(NRows, NColumns, MatrixElements1, MatrixElements2, CurrentRow, NColumns, SumRow) :- !,
  592    c_store(MatrixElements2[CurrentRow], SumRow),
  593    SumRow1 is 0,
  594    CurrentRow1 is CurrentRow+1,
  595    agg_lines_matrix(NRows,
  596                     NColumns,
  597                     MatrixElements1,
  598                     MatrixElements2,
  599                     CurrentRow1,
  600                     0,
  601                     SumRow1).
  602
  603agg_lines_matrix(NRows, NColumns, MatrixElements1, MatrixElements2, CurrentRow, CurrentColumn, SumRow) :-
  604    Index is CurrentRow*NColumns+CurrentColumn,
  605    c_load(MatrixElements1[Index], Element),
  606    SumRow1 is SumRow+Element,
  607    CurrentColumn1 is CurrentColumn+1,
  608    agg_lines_matrix(NRows,
  609                     NColumns,
  610                     MatrixElements1,
  611                     MatrixElements2,
  612                     CurrentRow,
  613                     CurrentColumn1,
  614                     SumRow1).
  615
  616
  617matrix_agg_lines(matrix(Type, [NRows, NColumns], MatrixElements1), AggregationMatrix) :-
  618    matrix_new(Type, [NRows, 1], AggregationMatrix),
  619    AggregationMatrix=matrix(_, _, MatrixElements2),
  620    agg_lines_matrix(NRows, NColumns, MatrixElements1, MatrixElements2, 0, 0, 0).
 matrix_agg_cols(+M1, +AggregationMatrix)
matrix_agg_cols (M,AggregationMatrix) sum all the elements of every column of M and the new value will be insert in the AggregationMatrix.
  624agg_cols_matrix(_, NColumns, _, _, _, NColumns, _) :- !.
  625
  626agg_cols_matrix(NRows, NColumns, MatrixElements1, MatrixElements2, NRows, CurrentColumn, SumRow) :- !,
  627    c_store(MatrixElements2[CurrentColumn], SumRow),
  628    SumRow1 is 0,
  629    CurrentColumn1 is CurrentColumn+1,
  630    agg_cols_matrix(NRows,
  631                    NColumns,
  632                    MatrixElements1,
  633                    MatrixElements2,
  634                    0,
  635                    CurrentColumn1,
  636                    SumRow1).
  637
  638agg_cols_matrix(NRows, NColumns, MatrixElements1, MatrixElements2, CurrentRow, CurrentColumn, SumRow) :-
  639    Index is CurrentRow*NColumns+CurrentColumn,
  640    c_load(MatrixElements1[Index], Element),
  641    SumRow1 is SumRow+Element,
  642    CurrentRow1 is CurrentRow+1,
  643    agg_cols_matrix(NRows,
  644                    NColumns,
  645                    MatrixElements1,
  646                    MatrixElements2,
  647                    CurrentRow1,
  648                    CurrentColumn,
  649                    SumRow1).
  650
  651
  652matrix_agg_cols(matrix(Type, [NRows, NColumns], MatrixElements), AggregationMatrix) :-
  653    matrix_new(Type, [1, NColumns], AggregationMatrix),
  654    AggregationMatrix=matrix(_, _, MatrixElements2),
  655    agg_cols_matrix(NRows, NColumns, MatrixElements, MatrixElements2, 0, 0, 0).
 matrix_op(+M1, +M2, +Operator, -M3)
matrix_op(M1,M2,Op,M3) apply Op an operator entered in input (+,-,*,/) element by element of matrix M1 and M2 and the result will be insert in M3.
  659matrix_map(Predicate, matrix(Type, [NRows, NColumns], MatrixElements1), matrix(Type, [NRows, NColumns], MatrixElements2), matrix(Type, [NRows, NColumns], MatrixElements3)) :-
  660    matrix_new(Type,
  661               [NRows, NColumns],
  662               matrix(_, _, MatrixElements3)),
  663    matrix_map(Predicate,
  664               NRows,
  665               NColumns,
  666               MatrixElements1,
  667               MatrixElements2,
  668               MatrixElements3,
  669               0,
  670               0).
  671
  672matrix_map(_, NRows, _, _, _, _, NRows, _) :- !.
  673
  674matrix_map(Predicate, NRows, NColumns, MatrixElements1, MatrixElements2, MatrixElements3, CurrentRow, NColumns) :- !,
  675    CurrentRow1 is CurrentRow+1,
  676    matrix_map(Predicate,
  677               NRows,
  678               NColumns,
  679               MatrixElements1,
  680               MatrixElements2,
  681               MatrixElements3,
  682               CurrentRow1,
  683               0).
  684
  685matrix_map(Predicate, NRows, NColumns, MatrixElements1, MatrixElements2, MatrixElements3, CurrentRow, CurrentColumn) :-
  686    Index is CurrentRow*NColumns+CurrentColumn,
  687    c_load(MatrixElements1[Index], Element1),
  688    c_load(MatrixElements2[Index], Element2),
  689    call(Predicate, Element1, Element2, Element3),
  690    c_store(MatrixElements3[Index], Element3),
  691    CurrentColumn1 is CurrentColumn+1,
  692    matrix_map(Predicate,
  693               NRows,
  694               NColumns,
  695               MatrixElements1,
  696               MatrixElements2,
  697               MatrixElements3,
  698               CurrentRow,
  699               CurrentColumn1).
  700
  701matrix_op(matrix(Type, [NRows, NColumns], Element1), matrix(Type, [NRows, NColumns], Element2), Operator, Matrix3) :-
  702    LambdaExpression= \E1^E2^E3^(Predicate=..[Operator, E1, E2], E3 is Predicate),
  703    matrix_map(LambdaExpression,
  704               matrix(Type, [NRows, NColumns], Element1),
  705               matrix(Type, [NRows, NColumns], Element2),
  706               Matrix3).
 matrix_op_to_all(+M1, +Operator, +Operand, -M2)
matrix_op_to_all(M1,Operator,Operand,M2) apply Operator (+,-,*,/) and Operand (an integer, a float...) element by element at matrix M1 and the result will be insert in M2. Example: input (+3) sum 3 to each elements of M1.
  710matrix_op_to_all(matrix(Type, [NRows, NColumns], MatrixElements), Operator, Operand, Matrix2) :-
  711    LambdaExpression= \E1^E2^(Predicate=..[Operator, E1, Operand], E2 is Predicate),
  712    matrix_map(LambdaExpression,
  713               matrix(Type, [NRows, NColumns], MatrixElements),
  714               Matrix2).
  715
  716% function that verifies that the inserted rows or columns are valid
  717function_call(_, []) :- !.
  718function_call(P, [H|T]) :-
  719    call(P, H),
  720    function_call(P, T).
  721
  722function_call_lesser(L, M) :-
  723    P= \H^(H<M),
  724    function_call(P, L).
  725
  726% matrix_select_rows(M1,RowList,M2) select a list of row (RowList) from M1 and this rows will be copied in M2.
  727matrix_select_rows(matrix(Type, [NRows, NColumns], MatrixElements1), RowList, matrix(Type, [NRows2, NColumns], MatrixElements2)) :-
  728    function_call_lesser(RowList, NRows),
  729    length(RowList, NRows2),
  730    matrix_new(Type,
  731               [NRows2, NColumns],
  732               matrix(_, _, MatrixElements2)),
  733    matrix_select_rows(NColumns, MatrixElements1, MatrixElements2, RowList, 0, 0).
  734
  735matrix_select_rows(_, _, _, [], _, _) :- !.
  736matrix_select_rows(NColumns, MatrixElements1, MatrixElements2, [_|MoreRows], RowIndex, NColumns) :- !,
  737    RowIndex1 is RowIndex+1,
  738    matrix_select_rows(NColumns, MatrixElements1, MatrixElements2, MoreRows, RowIndex1, 0).
  739
  740matrix_select_rows(NColumns, MatrixElements1, MatrixElements2, [Row|MoreRows], RowIndex, ColumnIndex) :-
  741    Index1 is Row*NColumns+ColumnIndex,
  742    c_load(MatrixElements1[Index1], E),
  743    Index2 is RowIndex*NColumns+ColumnIndex,
  744    c_store(MatrixElements2[Index2], E),
  745    ColumnIndex1 is ColumnIndex+1,
  746    matrix_select_rows(NColumns,
  747                       MatrixElements1,
  748                       MatrixElements2,
  749                       [Row|MoreRows],
  750                       RowIndex,
  751                       ColumnIndex1).
  752
  753% matrix_select_cols(M1,ColList,M2) select a list of columns (ColList) from M1 and this columns will be copied in M2.
  754matrix_select_cols(matrix(Type, [NRows, NColumnsMatrix1], MatrixElements1), List, matrix(Type, [NRows, NColumnsMatrix2], MatrixElements2)) :-
  755    function_call_lesser(List, NColumnsMatrix1),
  756    length(List, NColumnsMatrix2),
  757    matrix_new(Type,
  758               [NRows, NColumnsMatrix2],
  759               matrix(_, _, MatrixElements2)),
  760    matrix_select_cols(NRows,
  761                       NColumnsMatrix1,
  762                       NColumnsMatrix2,
  763                       MatrixElements1,
  764                       MatrixElements2,
  765                       List,
  766                       0,
  767                       0).
  768
  769matrix_select_cols(_, _, _, _, _, [], _, _) :- !.
  770matrix_select_cols(NRows, NColumnsMatrix1, NColumnsMatrix2, MatrixElements1, MatrixElements2, [_|MoreCols], NRows, ColumnIndex) :- !,
  771    NewNColumns is ColumnIndex+1,
  772    matrix_select_cols(NRows,
  773                       NColumnsMatrix1,
  774                       NColumnsMatrix2,
  775                       MatrixElements1,
  776                       MatrixElements2,
  777                       MoreCols,
  778                       0,
  779                       NewNColumns).
  780
  781matrix_select_cols(NRows, NColumnsMatrix1, NColumnsMatrix2, MatrixElements1, MatrixElements2, [Cols|MoreCols], RowIndex, ColumnIndex) :-
  782    Index1 is RowIndex*NColumnsMatrix1+Cols,
  783    c_load(MatrixElements1[Index1], E),    Index2 is RowIndex*NColumnsMatrix2+ColumnIndex,    c_store(MatrixElements2[Index2], E),    RowIndex1 is RowIndex+1,    matrix_select_cols(NRows,                       NColumnsMatrix1,                       NColumnsMatrix2,                       MatrixElements1,                       MatrixElements2,                       [Cols|MoreCols],                       RowIndex1,                       ColumnIndex).
 matrix_select(+M1, +Number, +List, +Result)
matrix_select(M1,Number,List,Result) select from M1 a list of rows or columns (to select rows Number will be 1, to select columns Number will be 2) and these will be copied in Result.
  798matrix_select(matrix(Type, [NRows, NColumns], MatrixElements), 1, List, Result) :-
  799    matrix_select_rows(matrix(Type, [NRows, NColumns], MatrixElements),
  800                       List,
  801                       Result).
  802
  803matrix_select(matrix(Type, [NRows, NColumns], MatrixElements), 2, List, Result) :-
  804    matrix_select_cols(matrix(Type, [NRows, NColumns], MatrixElements),
  805                       List,
  806                       Result).
 matrix_createZeroes(+Type, +Dimension, -M1)
matrix_createZeroes(M1) create a new matrix with all elements set to 0.
  810matrix_createZeroes(Type, [NRows, NColumns], M1) :-
  811    matrix_new(Type, [NRows, NColumns], M1),
  812    M1=matrix(Type, [NRows, NColumns], Elements),
  813    matrixSetZeroes(NRows, NColumns, Elements).
 matrix_read(+Dimension, -M1)
Read from keyboard values and these will be insert in M1.
  817matrix_read([NRows, NColumns], M1) :-
  818    matrix_new(doubles, [NRows, NColumns], M1),
  819    matrix_foreach(M1, \_^Y^read(Y)).
 matrixSetZeroes(+NRows, +NColumns, +Elements)
matrixSetZeroes sets to 0 all its argument matrix's elements.
  823matrixSetZeroes(NRows, NColumns, Elements) :-
  824    matrixSetAll(NRows, NColumns, Elements, 0).
 matrixSetOnes(+NRows, +NColumns, +Elements)
matrixSetOnes(M1) create a new matrix with all elements set to 1.
  828matrixSetOnes(NRows, NColumns, Elements) :-
  829    matrixSetAll(NRows, NColumns, Elements, 1).
 matrix_setAllNumbers(+NRows, +NColumns, +Elements)
matrix_setAll sets to Number given in input all the matrix's elements.
  833matrix_setAllNumbers(NRows, NColumns, Elements) :-
  834    read(Number),
  835    matrixSetAll(NRows, NColumns, Elements, Number).
 matrix_eye(+NRows, +NColumns, +Elements)
matrix_eye create a identity matrix
  839matrix_eye(NRows, NColumns, Elements) :-
  840    matrixEye(NRows, NColumns, Elements).
 matrix_equal(+M1, +M2)
matrix_equal(M1,M2) compares two matrix (M1 and M2).
  844equal_matrix(NRows, _, _, _, NRows, _) :- !.
  845
  846equal_matrix(NRows, NColumns, MatrixElements1, MatrixElements2, CurrentRow, NColumns) :- !,
  847    CurrentRow1 is CurrentRow+1,
  848    equal_matrix(NRows, NColumns, MatrixElements1, MatrixElements2, CurrentRow1, 0).
  849
  850equal_matrix(NRows, NColumns, MatrixElements1, MatrixElements2, CurrentRow, CurrentColumn) :-
  851    Index is CurrentRow*NColumns+CurrentColumn,
  852    c_load(MatrixElements1[Index], Element1),
  853    c_load(MatrixElements2[Index], Element2),
  854    CurrentColumn1 is CurrentColumn+1,
  855    Element1==Element2,
  856    equal_matrix(NRows, NColumns, MatrixElements1, MatrixElements2, CurrentRow, CurrentColumn1).
  857
  858matrix_equal(matrix(Type, [NRows, NColumns], MatrixElements1), matrix(Type, [NRows, NColumns], MatrixElements2)) :-
  859    equal_matrix(NRows, NColumns, MatrixElements1, MatrixElements2, 0, 0).
 matrix_from_list(+List, +M1)
matrix_from_list(List,M1) sect the elements of the matrix M1 with new element passed in a List.
  863from_list(List, Elements) :-
  864    from_list(List, 0, Elements).
  865
  866from_list([], _, _) :- !.
  867from_list([Element|MoreElements], Index, Elements) :-
  868    c_store(Elements[Index], Element),
  869    Index1 is Index+1,
  870    from_list(MoreElements, Index1, Elements).
  871
  872matrix_from_list(List, matrix(_, [_, _], Elements)) :-
  873    from_list(List, Elements).
  874
  875example(Elements) :-
  876    matrix_new(doubles, [3, 3], Matrix1),
  877    Matrix1=matrix(_, [3, 3], Elements),
  878    matrix_eye(3, 3, Elements).
  879
  880
  881matrix_random(NRows, NColumns, Matrix) :-
  882    matrix_new(doubles, [NRows, NColumns], Matrix),
  883    matrix_foreach(Matrix,
  884                   \X^Y^(random(E), myset(X, Y, E))).
  885
  886aux1 :-
  887    matrix_random(2, 2, Matrix1),
  888    matrix_transpose(Matrix1, Matrix2),
  889    matrix_mul(Matrix1, Matrix2, Matrix3),
  890    matrix_write(Matrix3).
  891
  892%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%       START TESTS       %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  893:- begin_tests(matrix).  894
  895    test(matrix_mul1) :-
  896    matrix_new(doubles, [2, 2], [1, 2, 3, 4], Matrix1),
  897    matrix_new(doubles, [2, 2], [4, 3, 2, 1], Matrix2),
  898    matrix_new(doubles, [2, 2], [13.0, 20.0, 5.0, 8.0], Expected),
  899    matrix_mul(Matrix1, Matrix2, Matrix3),
  900    matrix_equal(Matrix3, Expected).
  901
  902    test(matrix_dims1) :-
  903    matrix_new(doubles, [3, 3], [1, 2, 3, 4, 5, 6, 7, 8, 9], Matrix1),
  904    matrix_dims(Matrix1, [3, 3]).
  905
  906    test(matrix_size1) :-
  907    matrix_new(doubles, [3, 3], [1, 2, 3, 4, 5, 6, 7, 8, 9], Matrix1),
  908    matrix_size(Matrix1, 9).
  909    
  910    test(matrix_type1) :-
  911    matrix_new(doubles, [3, 3], [1, 2, 3, 4, 5, 6, 7, 8, 9], Matrix1),
  912    matrix_type(Matrix1, doubles).
  913
  914    test(matrix_to_list1) :-
  915    matrix_new(doubles, [3, 3], [1, 2, 3, 4, 5, 6, 7, 8, 9], Matrix1),
  916    matrix_to_list(Matrix1, [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0]).
  917
  918    test(matrix_get1) :-
  919    matrix_new(doubles, [3, 3], [1, 2, 3, 4, 5, 6, 7, 8, 9], Matrix1),
  920    matrix_get(Matrix1, [0, 1], 2.0).
  921
  922    test(matrix_set1) :-
  923    matrix_new(doubles, [3, 3], [1, 2, 3, 4, 5, 6, 7, 8, 9], Matrix1),
  924    matrix_set(Matrix1, [1, 2], 5),
  925    matrix_get(Matrix1, [1, 2], 5.0).
  926
  927    test(matrix_set_all1) :-
  928    matrix_new(doubles, [2, 2], Matrix1),
  929    matrix_new(doubles, [2, 2], [7.0, 7.0, 7.0, 7.0], Expected),
  930    matrix_set_all(number(doubles, 7), Matrix1),
  931    matrix_equal(Matrix1, Expected).
  932
  933    test(matrix_add1) :-
  934    matrix_new(doubles, [3, 3], [1, 2, 3, 4, 5, 6, 7, 8, 9], Matrix1),
  935    matrix_add(Matrix1, [1, 0], 7),
  936    matrix_get(Matrix1, [1, 0], 11.0).
  937
  938    test(matrix_inc1) :-
  939    matrix_new(doubles, [3, 3], [1, 2, 3, 4, 5, 6, 7, 8, 9], Matrix1),
  940    matrix_inc(Matrix1, [0, 1], 3.0).
  941
  942    test(matrix_dec1) :-
  943    matrix_new(doubles, [3, 3], [1, 2, 3, 4, 5, 6, 7, 8, 9], Matrix1),
  944    matrix_dec(Matrix1, [0, 1], 1.0).
  945
  946    test(matrix_max1) :-
  947    matrix_new(doubles, [3, 3], [1, 2, 3, 11, 5, 6, 7, 8, 9], Matrix1),
  948    matrix_max(Matrix1, 11.0).
  949
  950    test(matrix_min1) :-
  951    matrix_new(doubles, [3, 3], [3, 2, 1, 11, 5, 6, 7, 8, 9], Matrix1),
  952    matrix_min(Matrix1, 1.0).
  953
  954    test(matrix_sum1) :-
  955    matrix_new(doubles, [3, 3], [1, 2, 3, 11, 5, 6, 7, 8, 9], Matrix1),
  956    matrix_sum(Matrix1, 52.0).
  957
  958    test(matrix_transpose1) :-
  959    matrix_new(doubles, [2, 3], [1, 2, 3, 11, 5, 6], Matrix1),
  960    matrix_new(doubles, [3, 2], [1, 11, 2, 5, 3, 6], Expected),
  961    matrix_transpose(Matrix1, Matrix2),
  962    matrix_equal(Matrix2, Expected).
  963
  964    test(matrix_maxarg1) :-
  965    matrix_new(doubles, [2, 3], [1, 2, 3, 11, 5, 6], Matrix1),
  966    matrix_maxarg(Matrix1, [1, 0]).
  967
  968    test(matrix_minarg1) :-
  969    matrix_new(doubles, [2, 3], [13, 2, 3, 11, 1, 6], Matrix1),
  970    matrix_minarg(Matrix1, [1, 1]).
  971
  972    test(matrix_agg_lines1) :-
  973    matrix_new(doubles, [2, 3], [13, 2, 3, 10, 1, 6], Matrix1),
  974    matrix_new(doubles, [2, 1], [18, 17], Expected),
  975    matrix_agg_lines(Matrix1, AggregationMatrix),
  976    matrix_equal(AggregationMatrix, Expected).
  977
  978    test(matrix_agg_cols1) :-
  979    matrix_new(doubles, [2, 3], [13, 2, 3, 10, 1, 6], Matrix1),
  980    matrix_new(doubles, [1, 3], [23, 3, 9], Expected),
  981    matrix_agg_cols(Matrix1, AggregationMatrix),
  982    matrix_equal(AggregationMatrix, Expected).
  983
  984    test(matrix_op1) :-
  985    matrix_new(doubles, [2, 2], [13, 2, 3, 10], Matrix1),
  986    matrix_new(doubles, [2, 2], [2, 8, 13, 9], Matrix2),
  987    matrix_write(Matrix1),
  988    matrix_op(Matrix1, Matrix2, +, Matrix3),
  989    matrix_new(doubles, [2, 2], [15, 10, 16, 19], Expected),
  990    matrix_equal(Matrix3, Expected).
  991
  992    test(matrix_op_to_all1) :-
  993    matrix_new(doubles, [2, 3], [13, 2, 3, 10, 1, 6], Matrix1),
  994    matrix_new(doubles, [2, 3], [16, 5, 6, 13, 4, 9], Expected),
  995    matrix_op_to_all(Matrix1, +, 3, Matrix2),
  996    matrix_equal(Matrix2, Expected).
  997
  998    test(matrix_select1) :-
  999    matrix_new(doubles, [3, 2], [13, 2, 3, 10, 1, 6], Matrix1),
 1000    matrix_new(doubles, [1, 2], [13, 2], Expected),
 1001    matrix_select(Matrix1, 1, [0], Matrix2),
 1002    matrix_equal(Matrix2, Expected).
 1003
 1004    test(matrix_select2) :-
 1005    matrix_new(doubles, [3, 2], [13, 2, 3, 10, 1, 6], Matrix1),
 1006    matrix_new(doubles, [2, 2], [13, 2, 3, 10], Expected),
 1007    matrix_select(Matrix1, 1, [0, 1], Matrix2),
 1008    matrix_equal(Matrix2, Expected).
 1009
 1010    test(matrix_select3) :-
 1011    matrix_new(doubles, [2, 3], [13, 2, 3, 10, 1, 6], Matrix1),
 1012    matrix_new(doubles, [2, 2], [13, 3, 10, 6], Expected),
 1013    matrix_select(Matrix1, 2, [0, 2], Matrix2),
 1014    matrix_equal(Matrix2, Expected).
 1015
 1016    test(matrixCreateZeroes1) :-
 1017    matrix_createZeroes(doubles, [3, 3], Matrix1),
 1018    matrix_new(doubles, [3, 3], [0, 0, 0, 0, 0, 0, 0, 0, 0], Expected),
 1019    matrix_equal(Matrix1, Expected).
 1020
 1021    test(matrixSetZeroes1) :-
 1022    matrix_new(doubles, [3, 3], Matrix1),
 1023    Matrix1=matrix(_, [3, 3], Elements),
 1024    matrix_new(doubles, [3, 3], [0, 0, 0, 0, 0, 0, 0, 0, 0], Expected),
 1025    matrixSetZeroes(3, 3, Elements),
 1026    matrix_equal(Matrix1, Expected).
 1027
 1028    test(matrixSetOnes1) :-
 1029    matrix_new(doubles, [3, 3], Matrix1),
 1030    Matrix1=matrix(_, [3, 3], Elements),
 1031    matrix_new(doubles, [3, 3], [1, 1, 1, 1, 1, 1, 1, 1, 1], Expected),
 1032    matrixSetOnes(3, 3, Elements),
 1033    matrix_equal(Matrix1, Expected).
 1034
 1035    test(matrix_setAllNumbers1) :-
 1036    matrix_new(doubles, [3, 3], Matrix1),
 1037    Matrix1=matrix(_, [3, 3], Elements),
 1038    Number is 5.0,
 1039    matrixSetAll(3, 3, Elements, Number),
 1040    matrix_new(doubles, [3, 3], [5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0], Expected),
 1041    matrix_equal(Matrix1, Expected).
 1042
 1043    test(matrix_eye1) :-
 1044    matrix_new(doubles, [3, 3], Matrix1),
 1045    Matrix1=matrix(_, [3, 3], Elements),
 1046    matrix_eye(3, 3, Elements),
 1047    matrix_new(doubles, [3, 3], [1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0], Expected),
 1048    matrix_equal(Matrix1, Expected).
 1049
 1050    test(matrix_equal1) :-
 1051    matrix_new(doubles, [3, 3], [1, 2, 3, 4, 5, 6, 7, 8, 9], Matrix1),
 1052    matrix_new(doubles, [3, 2], [1, 3, 4, 6, 7, 9], Expected),
 1053    matrix_select(Matrix1, 2, [0, 2], Matrix2),
 1054    matrix_equal(Matrix2, Expected).
 1055
 1056    test(matrix_map1) :-
 1057    matrix_new(doubles, [2, 3], [1, 2, 3, 11, 5, 6], Matrix1),
 1058    matrix_map(\E1^E2^(E2 is E1*2),
 1059               Matrix1,
 1060               Matrix2),
 1061    matrix_new(doubles, [2, 3], [2, 4, 6, 22, 10, 12], Expected),
 1062    matrix_equal(Matrix2, Expected).
 1063
 1064    test(matrix_foreach1) :-
 1065    matrix_new(doubles, [3, 3], [1, 2, 3, 4, 5, 6, 7, 8, 9], Matrix1),
 1066    matrix_new(doubles, [3, 3], [3, 3, 3, 3, 3, 3, 3, 3, 3], Expected),
 1067    matrix_foreach(Matrix1,
 1068                   \X^Y^myset(X, Y, 3)),
 1069    matrix_equal(Matrix1, Expected).
 1070
 1071    test(matrix_foldl1) :-
 1072    matrix_new(doubles, [3, 3], [1, 2, 3, 4, 5, 6, 7, 8, 9], Matrix1),
 1073    matrix_foldl(Matrix1, findMax, 9.0).
 1074
 1075    test(matrix_from_list1) :-
 1076    matrix_new(doubles, [3, 3], [1, 2, 3, 4, 5, 6, 7, 8, 9], Matrix1),
 1077    matrix_from_list([11, 12, 13, 14, 15, 16, 17, 18, 19], Matrix1),
 1078    matrix_new(doubles,
 1079               [3, 3],
 1080               [11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 17.0, 18.0, 19.0],
 1081               Expected),
 1082    matrix_equal(Matrix1, Expected).
 1083
 1084    test(cholesky1) :-
 1085    matrix_new(doubles, [3, 3], [25, 15, -5, 15, 18, 0, -5, 0, 11], Matrix1),
 1086    matrix_new(doubles, [3, 3], Matrix2),
 1087    matrix_new(doubles, [3, 3], [5.0, 0.0, 0.0, 3.0, 3.0, 0.0, -1.0, 1.0, 3.0], Expected),
 1088    matrix_sollower(Matrix1, Matrix2),
 1089    cholesky(Matrix2),
 1090    matrix_equal(Matrix2, Expected).
 1091
 1092end_tests(matrix)