1:- module(matrix,
    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          ]).   41
   42:- use_module(library(ffi)).   43:- use_module(library(lambda)).   44
   45
   46
   47% https://www.dcc.fc.up.pt/~vsc/Yap/documentation.html#matrix
   48:- c_import("#include <atlas/cblas.h>",
   49            [libblas],
   50            
   51            [ cblas_dgemm(int,
   52                          int,
   53                          int,
   54                          int,
   55                          int,
   56                          int,
   57                          double,
   58                          *double,
   59                          int,
   60                          *double,
   61                          int,
   62                          double,
   63                          *double,
   64                          int)
   65            ]).
   66
   67:- (multifile user:term_expansion/2).   68:- (dynamic user:term_expansion/2).   69                
   70
   71user:term_expansion((:-import),  (:-Import)) :-
   72    current_prolog_flag(arch, Arch),
   73    atom_concat('../lib/', Arch, A1),
   74    atom_concat(A1, '/matrixNative', A2),
   75    Import=c_import("#include \"../matrixNative.h\"", [A2], [matrixSetAll(int, int, *double, double), matrixEye(int, int, *double)]).
   76
   77
   78:- import.
 matrix_write(+M1)
Print elements of a matrix
   82matrix_write(matrix(_, [NRows, NColumns], Matrix)) :-
   83    write_matrix(Matrix, NRows, NColumns).
   84
   85write_matrix(Matrix, NRows, NColumns) :-
   86    write_matrix(Matrix, 0, 0, NRows, NColumns).
   87
   88write_matrix(_, NRows, _, NRows, _) :-
   89    !.
   90
   91write_matrix(MatrixElements, CurrentRow, NColumns, NRows, NColumns) :-
   92    !,
   93    writeln(''),
   94    CurrentRow1 is CurrentRow+1,
   95    write_matrix(MatrixElements, CurrentRow1, 0, NRows, NColumns).
   96
   97write_matrix(MatrixElements, CurrentRow, CurrentColumn, NRows, NColumns) :-
   98    !,
   99    Index is CurrentRow*NColumns+CurrentColumn,
  100    c_load(MatrixElements[Index], Element),
  101    write(Element),
  102    write(' '),
  103    CurrentColumn1 is CurrentColumn+1,
  104    write_matrix(MatrixElements, CurrentRow, CurrentColumn1, NRows, NColumns).
  105
  106%!  matrix_new(+Type,+Dimension,-M1)
  107% 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.
  108matrix_new(doubles, [NRows, NColumns], matrix(doubles, [NRows, NColumns], Matrix)) :-
  109    NRows>0,
  110    NColumns>0,
  111    Nelements is NRows*NColumns,
  112    c_alloc(Matrix, double[Nelements]).
  113
  114%!  matrix_new(+Type,+Dimension,+Elements,-M1)
  115% 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.
  116matrix_new(doubles, [NRows, NColumns], ListElements, matrix(doubles, [NRows, NColumns], Matrix)) :-
  117    NRows>0,
  118    NColumns>0,
  119    NElements is NRows*NColumns,
  120    length(ListElements, NElements),
  121    c_alloc(Matrix, double[]=ListElements).
 matrix_mul(+M1, +M2, -M3)
Multiplication between two matrices using "cblas_dgemm" function from "cblas.h"
  125matrix_mul(matrix(Type, [NRowsMatrix1, NColumnsMatrix1], Matrix1), matrix(Type, [NColumnsMatrix1, NColumnsMatrix2], Matrix2), matrix(Type, [NRowsMatrix1, NColumnsMatrix2], Matrix3)) :-
  126    matrix_new(Type,
  127               [NRowsMatrix1, NColumnsMatrix2],
  128               matrix(Type, [NRowsMatrix1, NColumnsMatrix2], Matrix3)),
  129    cblas_dgemm(102,
  130                111,
  131                111,
  132                NRowsMatrix1,
  133                NColumnsMatrix2,
  134                NColumnsMatrix1,
  135                1,
  136                Matrix1,
  137                NRowsMatrix1,
  138                Matrix2,
  139                NColumnsMatrix1,
  140                0,
  141                Matrix3,
  142                NRowsMatrix1).
 matrix_dims(+M1, -Dimension)
Calculates the total number of elements in the matrix. The value is contained in Dimension.
  146matrix_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.
  150matrix_size(matrix(_, [NRows, NColumns], _), Nelements) :-
  151    Nelements is NRows*NColumns.
 matrix_type(+M1, -Type)
Function that return the Type of the matrix passed.
  155matrix_type(matrix(Type, _, _), Type).
 matrix_to_list(+M1, -List)
Function that transforms the matrix passed in input in a list.
  159to_list(_, NRows, _, NRows, _, List) :-
  160    !,
  161    List=([]).
  162
  163to_list(MatrixElements, CurrentRow, NColumns, Nrows, NColumns, List) :-
  164    !,
  165    CurrentRow1 is CurrentRow+1,
  166    to_list(MatrixElements, CurrentRow1, 0, Nrows, NColumns, List).
  167
  168to_list(MatrixElements, CurrentRow, CurrentColumn, NRows, NColumns, [Element|MoreElements]) :-
  169    Index is CurrentRow*NColumns+CurrentColumn,
  170    c_load(MatrixElements[Index], Element),
  171    CurrentColumn1 is CurrentColumn+1,
  172    to_list(MatrixElements, CurrentRow, CurrentColumn1, NRows, NColumns, MoreElements).
  173
  174matrix_to_list(matrix(_, [NRows, NColumns], MatrixElements), List) :-
  175    to_list(MatrixElements, 0, 0, NRows, NColumns, List).
  176
  177%!  matrix_get(+M1,+Position,-Element)
  178% matrix_get(M,[Row,Column],E) unifies E with the element of M at position [Row, Column] 
  179matrix_get(matrix(_, [_, NColumns], MatrixElements), [Row, Column], Element) :-
  180    Index is Row*NColumns+Column,
  181    c_load(MatrixElements[Index], Element).
  182
  183%!  matrix_set(+Type,+Position,-Element)
  184% matrix_set(M,[Row,Column],E) update the element of M at position [Row, Column] with the element specified in Element.
  185matrix_set(matrix(_, [_, NColumns], MatrixElements), [Row, Column], Element) :-
  186    Index is Row*NColumns+Column,
  187    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.
  191matrix_map(Predicate, matrix(Type, [NRows, NColumns], MatrixElements1), matrix(Type, [NRows, NColumns], MatrixElements2)) :-
  192    matrix_new(Type,
  193               [NRows, NColumns],
  194               matrix(_, _, MatrixElements2)),
  195    matrix_map(Predicate, NRows, NColumns, MatrixElements1, MatrixElements2, 0, 0).
  196
  197matrix_map(_, NRows, _, _, _, NRows, _) :-
  198    !.
  199
  200matrix_map(Predicate, NRows, NColumns, MatrixElements1, MatrixElements2, CurrentRow, NColumns) :-
  201    !,
  202    CurrentRow1 is CurrentRow+1,
  203    matrix_map(Predicate, NRows, NColumns, MatrixElements1, MatrixElements2, CurrentRow1, 0).
  204
  205matrix_map(Predicate, NRows, NColumns, MatrixElements1, MatrixElements2, CurrentRow, CurrentColumn) :-
  206    Index is CurrentRow*NColumns+CurrentColumn,
  207    c_load(MatrixElements1[Index], E1),
  208    call(Predicate, E1, E2),
  209    c_store(MatrixElements2[Index], E2),
  210    CurrentColumn1 is CurrentColumn+1,
  211    matrix_map(Predicate,
  212               NRows,
  213               NColumns,
  214               MatrixElements1,
  215               MatrixElements2,
  216               CurrentRow,
  217               CurrentColumn1).
  218
  219%%%%%%%%%%%%%%%%%%%%%%%%%%%%       TWO FUNCTIONS USED BY MATRIX FOREACH       %%%%%%%%%%%%%%%%%%%%%%%%%%%%
  220
  221% Unify the element with the number specified.
  222myset(_, Number, Number).
  223
  224myset(_, 0).
  225% This function find the max number between two number.
  226findMax(FirstNumber, SecondNumber, FirstNumber) :-
  227    FirstNumber>=SecondNumber,
  228    !.
  229
  230findMax(_, SecondNumber, SecondNumber).
  231
  232%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 matrix_foreach(+M1, +Predicate)
For each element E of M1, enumerated by rows, calls Predicate(+E,-F) and replaces E with F.
  236matrix_foreach(matrix(_, [NRows, NColumns], MatrixElements), Predicate) :-
  237    matrix_foreach(Predicate, NRows, NColumns, MatrixElements, 0, 0).
  238
  239matrix_foreach(_, NRows, _, _, NRows, _) :-
  240    !.
  241
  242matrix_foreach(Predicate, NRows, NColumns, MatrixElements, CurrentRow, NColumns) :-
  243    !,
  244    CurrentRow1 is CurrentRow+1,
  245    matrix_foreach(Predicate, NRows, NColumns, MatrixElements, CurrentRow1, 0).
  246
  247matrix_foreach(Predicate, NRows, NColumns, MatrixElements, CurrentRow, CurrentColumn) :-
  248    Index is CurrentRow*NColumns+CurrentColumn,
  249    c_load(MatrixElements[Index], E),
  250    call(Predicate, E, E1),
  251    c_store(MatrixElements[Index], E1),
  252    CurrentColumn1 is CurrentColumn+1,
  253    matrix_foreach(Predicate,
  254                   NRows,
  255                   NColumns,
  256                   MatrixElements,
  257                   CurrentRow,
  258                   CurrentColumn1).
  259
  260%!  matrix_foldl(+M1,+Predicate,-Max)
  261% Fold a matrix, using arguments of the matrix as left argument.
  262matrix_foldl(matrix(_, [NRows, NColumns], MatrixElements), Predicate, Max) :-
  263    c_load(MatrixElements[0], TmpMax),
  264    matrix_start_foldl(Predicate,
  265                       NRows,
  266                       NColumns,
  267                       MatrixElements,
  268                       TmpMax,
  269                       Max).
  270
  271matrix_start_foldl(_, 1, 1, _, Max, Max) :-
  272    !.
  273matrix_start_foldl(Predicate, NRows, NColumns, MatrixElements, Tmp, Max) :-
  274    NColumns>1,
  275    !,
  276    matrix_foldl(Predicate,
  277                 NRows,
  278                 NColumns,
  279                 MatrixElements,
  280                 0,
  281                 1,
  282                 Tmp,
  283                 Max).
  284
  285matrix_start_foldl(Predicate, NRows, NColumns, MatrixElements, Tmp, Max) :-
  286    matrix_foldl(Predicate,
  287                 NRows,
  288                 NColumns,
  289                 MatrixElements,
  290                 1,
  291                 0,
  292                 Tmp,
  293                 Max).
  294
  295matrix_foldl(_, NRows, _, _, NRows, _, Max, Max) :-
  296    !.
  297
  298matrix_foldl(Predicate, NRows, NColumns, MatrixElements, CurrentRow, NColumns, TmpMax, Max) :-
  299    !,
  300    CurrentRow1 is CurrentRow+1,
  301    matrix_foldl(Predicate,
  302                 NRows,
  303                 NColumns,
  304                 MatrixElements,
  305                 CurrentRow1,
  306                 0,
  307                 TmpMax,
  308                 Max).
  309
  310matrix_foldl(Predicate, NRows, NColumns, MatrixElements, CurrentRow, CurrentColumn, TmpMax, Max) :-
  311    Index is CurrentRow*NColumns+CurrentColumn,
  312    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.
  326set_all(_, NRows, _, _, NRows, _) :-
  327    !.
  328
  329set_all(Number, NRows, NColumns, MatrixElements, CurrentRow, NColumns) :-
  330    !,
  331    CurrentRow1 is CurrentRow+1,
  332    set_all(Number, NRows, NColumns, MatrixElements, CurrentRow1, 0).
  333
  334set_all(Number, NRows, NColumns, MatrixElements, CurrentRow, CurrentColumn) :-
  335    Index is CurrentRow*NColumns+CurrentColumn,
  336    c_store(MatrixElements[Index], Number),
  337    CurrentColumn1 is CurrentColumn+1,
  338    set_all(Number, NRows, NColumns, MatrixElements, CurrentRow, CurrentColumn1).
  339
  340matrix_set_all(number(Type, Number), matrix(Type, [NRows, NColumns], Element)) :-
  341    set_all(Number, NRows, NColumns, Element, 0, 0).
  342
  343%!  matrix_add(+M1,Position,+Number)
  344% matrix_add(M,[Row,Column],Number) add Number to the element of matrix M in position [Row,Column].
  345matrix_add(matrix(_, [_, NColumns], MatrixElements), [Row, Column], Number) :-
  346    Index is Row*NColumns+Column,
  347    c_load(MatrixElements[Index], Tmp),
  348    NewElement is Tmp+Number,
  349    c_store(MatrixElements[Index], NewElement).
  350
  351%!  matrix_inc(+M1,+Position)
  352% matrix_inc(M,[Row,Column]) add 1 (one) to the element of matrix M in position [Row,Column].
  353matrix_inc(matrix(_, [_, NColumns], MatrixElements), [Row, Column]) :-
  354    Index is Row*NColumns+Column,
  355    c_load(MatrixElements[Index], Tmp),
  356    E is Tmp+1,
  357    c_store(MatrixElements[Index], E).
  358
  359%!  matrix_inc(+M1,+Position,-NewElement)
  360% 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.
  361matrix_inc(matrix(_, [_, NColumns], MatrixElements), [Row, Column], NewElement) :-
  362    Index is Row*NColumns+Column,
  363    c_load(MatrixElements[Index], Tmp),
  364    NewElement is Tmp+1,
  365    c_store(MatrixElements[Index], NewElement).
  366
  367%!  matrix_dec(+M1,+Position)
  368% matrix_dec(M,[Row,Column]) decrement 1 (one) to the element of matrix M in position [Row,Column].
  369matrix_dec(matrix(_, [_, NColumns], MatrixElements), [Row, Column]) :-
  370    Index is Row*NColumns+Column,
  371    c_load(MatrixElements[Index], Tmp),
  372    NewElement is Tmp-1,
  373    c_store(MatrixElements[Index], NewElement).
  374
  375%!  matrix_dec(+M1,+Position,-NewElement)
  376% 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.
  377matrix_dec(matrix(_, [_, NColumns], MatrixElements), [Row, Column], NewElement) :-
  378    Index is Row*NColumns+Column,
  379    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.
  385max_matrix(NRows, _, _, NRows, _, Max, Max) :-
  386    !.
  387
  388max_matrix(NRows, NColumns, MatrixElements, CurrentRow, NColumns, TmpMax, Max) :-
  389    !,
  390    CurrentRow1 is CurrentRow+1,
  391    max_matrix(NRows, NColumns, MatrixElements, CurrentRow1, 0, TmpMax, Max).
  392
  393max_matrix(NRows, NColumns, MatrixElements, CurrentRow, CurrentColumn, TmpMax, Max) :-
  394    Index is CurrentRow*NColumns+CurrentColumn,
  395    c_load(MatrixElements[Index], E),
  396    CurrentColumn1 is CurrentColumn+1,
  397    (   E>TmpMax
  398    ->  TmpMax1 is E
  399    ;   TmpMax1 is TmpMax
  400    ),
  401    max_matrix(NRows,
  402               NColumns,
  403               MatrixElements,
  404               CurrentRow,
  405               CurrentColumn1,
  406               TmpMax1,
  407               Max).
  408
  409matrix_max(matrix(_, [NRows, NColumns], MatrixElements), Max) :-
  410    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.
  415matrix_min(Matrix, Min) :-
  416    matrix_foldl(Matrix,
  417                 \X^Y^Z^(X=<Y->Z=X;Z=Y),
  418                 Min).
 matrix_sum(+M1, -Sum)
matrix_sum(M,Sum) unify Sum with the sum of the elements in M.
  422matrix_sum(Matrix, Sum) :-
  423    matrix_foldl(Matrix,
  424                 \CurrentSum^Element^NewSum^(NewSum is CurrentSum+Element),
  425                 Sum).
 matrix_transpose(+M1, -M2)
matrix_transpose(M1,M2) transpose matrix M1 to M2.
  429transpose_matrix(NRows, _, _, _, NRows, _) :-
  430    !.
  431
  432transpose_matrix(NRows, NColumns, MatrixElements1, MatrixElements2, CurrentRow, NColumns) :-
  433    !,
  434    CurrentRow1 is CurrentRow+1,
  435    transpose_matrix(NRows, NColumns, MatrixElements1, MatrixElements2, CurrentRow1, 0).
  436
  437transpose_matrix(NRows, NColumns, MatrixElements1, MatrixElements2, CurrentRow, CurrentColumn) :-
  438    Index is CurrentRow*NColumns+CurrentColumn,
  439    IndexMatrix2 is CurrentColumn*NRows+CurrentRow,
  440    c_load(MatrixElements1[Index], E),
  441    c_store(MatrixElements2[IndexMatrix2], E),
  442    CurrentColumn1 is CurrentColumn+1,
  443    transpose_matrix(NRows,
  444                     NColumns,
  445                     MatrixElements1,
  446                     MatrixElements2,
  447                     CurrentRow,
  448                     CurrentColumn1).
  449
  450
  451matrix_transpose(matrix(Type, [NRows, NColumns], MatrixElements1), Matrix2) :-
  452    matrix_new(Type, [NColumns, NRows], Matrix2),
  453    Matrix2=matrix(_, _, MatrixElements2),
  454    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.
  458maxarg_matrix(NRows, _, _, NRows, _, _, [NRowsMax, NColumnsMax], [NRowsMax, NColumnsMax]) :-
  459    !.
  460
  461maxarg_matrix(NRows, NColumns, MatrixElements, CurrentRow, NColumns, TmpMax, [NRowTmp, NColumnTmp], [NRowMax, NColomnMax]) :-
  462    !,
  463    CurrentRow1 is CurrentRow+1,
  464    maxarg_matrix(NRows,
  465                  NColumns,
  466                  MatrixElements,
  467                  CurrentRow1,
  468                  0,
  469                  TmpMax,
  470                  [NRowTmp, NColumnTmp],
  471                  [NRowMax, NColomnMax]).
  472
  473                maxarg_matrix(NRows, NColumns, MatrixElements, CurrentRow, CurrentColumn, TmpMax, [NRowTmp, NColumnTmp], [NRowMax, NColomnMax]) :-
  474    Index is CurrentRow*NColumns+CurrentColumn,
  475    c_load(MatrixElements[Index], Element),
  476    CurrentColumn1 is CurrentColumn+1,
  477    (   Element>=TmpMax
  478    ->  ( TmpMax1 is Element,
  479          NRowTmp1 is CurrentRow,
  480          NColumnTmp1 is CurrentColumn
  481        ),
  482        maxarg_matrix(NRows,
  483                      NColumns,
  484                      MatrixElements,
  485                      CurrentRow,
  486                      CurrentColumn1,
  487                      TmpMax1,
  488                      [NRowTmp1, NColumnTmp1],
  489                      [NRowMax, NColomnMax])
  490    ;   TmpMax1 is TmpMax,
  491        NRowTmp1 is NRowTmp,
  492        NColumnTmp1 is NColumnTmp,
  493        maxarg_matrix(NRows,
  494                      NColumns,
  495                      MatrixElements,
  496                      CurrentRow,
  497                      CurrentColumn1,
  498                      TmpMax1,
  499                      [NRowTmp1, NColumnTmp1],
  500                      [NRowMax, NColomnMax])
  501    ).
  502
  503matrix_maxarg(matrix(_, [NRows, NColumns], MatrixElements), [NRowMax, NColomnMax]) :-
  504    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.
  516minarg_matrix(NRows, _, _, NRows, _, _, [NRowTmp, NColumnTmp], [NRowTmp, NColumnTmp]) :-
  517    !.
  518
  519minarg_matrix(NRows, NColumns, MatrixElements, CurrentRow, NColumns, TmpMax, [NRowTmp, NColumnTmp], [NRowMax, NColumnMax]) :-
  520    !,
  521    CurrentRow1 is CurrentRow+1,
  522    minarg_matrix(NRows,
  523                  NColumns,
  524                  MatrixElements,
  525                  CurrentRow1,
  526                  0,
  527                  TmpMax,
  528                  [NRowTmp, NColumnTmp],
  529                  [NRowMax, NColumnMax]).
  530
  531                minarg_matrix(NRows, NColumns, MatrixElements, CurrentRow, CurrentColumn, TmpMax, [NRowTmp, NColumnTmp], [NRowMax, NColumnMax]) :-
  532    Index is CurrentRow*NColumns+CurrentColumn,
  533    c_load(MatrixElements[Index], Element),
  534    CurrentColumn1 is CurrentColumn+1,
  535    (   Element=<TmpMax
  536    ->  TmpMax1 is Element,
  537        NRowTmp1 is CurrentRow,
  538        NColumnTmp1 is CurrentColumn
  539    ;   TmpMax1 is TmpMax,
  540        NRowTmp1 is NRowTmp,
  541        NColumnTmp1 is NColumnTmp
  542    ),
  543    minarg_matrix(NRows,
  544                  NColumns,
  545                  MatrixElements,
  546                  CurrentRow,
  547                  CurrentColumn1,
  548                  TmpMax1,
  549                  [NRowTmp1, NColumnTmp1],
  550                  [NRowMax, NColumnMax]).
  551
  552matrix_minarg(matrix(_, [NRows, NColumns], MatrixElements), [NRowMax, NColumnMax]) :-
  553    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.
  565agg_lines_matrix(NRows, _, _, _, NRows, _, _) :-
  566    !.
  567
  568agg_lines_matrix(NRows, NColumns, MatrixElements1, MatrixElements2, CurrentRow, NColumns, SumRow) :-
  569    !,
  570    c_store(MatrixElements2[CurrentRow], SumRow),
  571    SumRow1 is 0,
  572    CurrentRow1 is CurrentRow+1,
  573    agg_lines_matrix(NRows,
  574                     NColumns,
  575                     MatrixElements1,
  576                     MatrixElements2,
  577                     CurrentRow1,
  578                     0,
  579                     SumRow1).
  580
  581agg_lines_matrix(NRows, NColumns, MatrixElements1, MatrixElements2, CurrentRow, CurrentColumn, SumRow) :-
  582    Index is CurrentRow*NColumns+CurrentColumn,
  583    c_load(MatrixElements1[Index], Element),
  584    SumRow1 is SumRow+Element,
  585    CurrentColumn1 is CurrentColumn+1,
  586    agg_lines_matrix(NRows,
  587                     NColumns,
  588                     MatrixElements1,
  589                     MatrixElements2,
  590                     CurrentRow,
  591                     CurrentColumn1,
  592                     SumRow1).
  593
  594
  595matrix_agg_lines(matrix(Type, [NRows, NColumns], MatrixElements1), AggregationMatrix) :-
  596    matrix_new(Type, [NRows, 1], AggregationMatrix),
  597    AggregationMatrix=matrix(_, _, MatrixElements2),
  598    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.
  602agg_cols_matrix(_, NColumns, _, _, _, NColumns, _) :-
  603    !.
  604
  605agg_cols_matrix(NRows, NColumns, MatrixElements1, MatrixElements2, NRows, CurrentColumn, SumRow) :-
  606    !,
  607    c_store(MatrixElements2[CurrentColumn], SumRow),
  608    SumRow1 is 0,
  609    CurrentColumn1 is CurrentColumn+1,
  610    agg_cols_matrix(NRows,
  611                    NColumns,
  612                    MatrixElements1,
  613                    MatrixElements2,
  614                    0,
  615                    CurrentColumn1,
  616                    SumRow1).
  617
  618agg_cols_matrix(NRows, NColumns, MatrixElements1, MatrixElements2, CurrentRow, CurrentColumn, SumRow) :-
  619    Index is CurrentRow*NColumns+CurrentColumn,
  620    c_load(MatrixElements1[Index], Element),
  621    SumRow1 is SumRow+Element,
  622    CurrentRow1 is CurrentRow+1,
  623    agg_cols_matrix(NRows,
  624                    NColumns,
  625                    MatrixElements1,
  626                    MatrixElements2,
  627                    CurrentRow1,
  628                    CurrentColumn,
  629                    SumRow1).
  630
  631
  632matrix_agg_cols(matrix(Type, [NRows, NColumns], MatrixElements), AggregationMatrix) :-
  633    matrix_new(Type, [1, NColumns], AggregationMatrix),
  634    AggregationMatrix=matrix(_, _, MatrixElements2),
  635    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.
  639matrix_map(Predicate, matrix(Type, [NRows, NColumns], MatrixElements1), matrix(Type, [NRows, NColumns], MatrixElements2), matrix(Type, [NRows, NColumns], MatrixElements3)) :-
  640    matrix_new(Type,
  641               [NRows, NColumns],
  642               matrix(_, _, MatrixElements3)),
  643    matrix_map(Predicate,
  644               NRows,
  645               NColumns,
  646               MatrixElements1,
  647               MatrixElements2,
  648               MatrixElements3,
  649               0,
  650               0).
  651
  652matrix_map(_, NRows, _, _, _, _, NRows, _) :-
  653    !.
  654
  655matrix_map(Predicate, NRows, NColumns, MatrixElements1, MatrixElements2, MatrixElements3, CurrentRow, NColumns) :-
  656    !,
  657    CurrentRow1 is CurrentRow+1,
  658    matrix_map(Predicate,
  659               NRows,
  660               NColumns,
  661               MatrixElements1,
  662               MatrixElements2,
  663               MatrixElements3,
  664               CurrentRow1,
  665               0).
  666
  667matrix_map(Predicate, NRows, NColumns, MatrixElements1, MatrixElements2, MatrixElements3, CurrentRow, CurrentColumn) :-
  668    Index is CurrentRow*NColumns+CurrentColumn,
  669    c_load(MatrixElements1[Index], Element1),
  670    c_load(MatrixElements2[Index], Element2),
  671    call(Predicate, Element1, Element2, Element3),
  672    c_store(MatrixElements3[Index], Element3),
  673    CurrentColumn1 is CurrentColumn+1,
  674    matrix_map(Predicate,
  675               NRows,
  676               NColumns,
  677               Element1,
  678               Element2,
  679               Element3,
  680               CurrentRow,
  681               CurrentColumn1).
  682
  683matrix_op(matrix(Type, [NRows, NColumns], Element1), matrix(Type, [NRows, NColumns], Element2), Operator, Matrix3) :-
  684    LambdaExpression= \E1^E2^E3^(Predicate=..[Operator, E1, E2], E3 is Predicate),
  685    matrix_map(LambdaExpression,
  686               matrix(Type, [NRows, NColumns], Element1),
  687               matrix(Type, [NRows, NColumns], Element2),
  688               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.
  692matrix_op_to_all(matrix(Type, [NRows, NColumns], MatrixElements), Operator, Operand, Matrix2) :-
  693    LambdaExpression= \E1^E2^(Predicate=..[Operator, E1, Operand], E2 is Predicate),
  694    matrix_map(LambdaExpression,
  695               matrix(Type, [NRows, NColumns], MatrixElements),
  696               Matrix2).
  697
  698% function that verifies that the inserted rows or columns are valid
  699function_call(_, []) :-
  700    !.
  701function_call(P, [H|T]) :-
  702    call(P, H),
  703    function_call(P, T).
  704
  705function_call_lesser(L, M) :-
  706    P= \H^(H<M),
  707    function_call(P, L).
  708
  709% matrix_select_rows(M1,RowList,M2) select a list of row (RowList) from M1 and this rows will be copied in M2.
  710matrix_select_rows(matrix(Type, [NRows, NColumns], MatrixElements1), RowList, matrix(Type, [NRows2, NColumns], MatrixElements2)) :-
  711    function_call_lesser(RowList, NRows),
  712    length(RowList, NRows2),
  713    matrix_new(Type,
  714               [NRows2, NColumns],
  715               matrix(_, _, MatrixElements2)),
  716    matrix_select_rows(NColumns, MatrixElements1, MatrixElements2, RowList, 0, 0).
  717
  718matrix_select_rows(_, _, _, [], _, _) :-
  719    !.
  720matrix_select_rows(NColumns, MatrixElements1, MatrixElements2, [_|MoreRows], RowIndex, NColumns) :-
  721    !,
  722    RowIndex1 is RowIndex+1,
  723    matrix_select_rows(NColumns, MatrixElements1, MatrixElements2, MoreRows, RowIndex1, 0).
  724
  725matrix_select_rows(NColumns, MatrixElements1, MatrixElements2, [Row|MoreRows], RowIndex, ColumnIndex) :-
  726    Index1 is Row*NColumns+ColumnIndex,
  727    c_load(MatrixElements1[Index1], E),
  728    Index2 is RowIndex*NColumns+ColumnIndex,
  729    c_store(MatrixElements2[Index2], E),
  730    ColumnIndex1 is ColumnIndex+1,
  731    matrix_select_rows(NColumns,
  732                       MatrixElements1,
  733                       MatrixElements2,
  734                       [Row|MoreRows],
  735                       RowIndex,
  736                       ColumnIndex1).
  737
  738% matrix_select_cols(M1,ColList,M2) select a list of columns (ColList) from M1 and this columns will be copied in M2.
  739matrix_select_cols(matrix(Type, [NRows, NColumnsMatrix1], MatrixElements1), List, matrix(Type, [NRows, NColumnsMatrix2], MatrixElements2)) :-
  740    function_call_lesser(List, NColumnsMatrix1),
  741    length(List, NColumnsMatrix2),
  742    matrix_new(Type,
  743               [NRows, NColumnsMatrix2],
  744               matrix(_, _, MatrixElements2)),
  745    matrix_select_cols(NRows,
  746                       NColumnsMatrix1,
  747                       NColumnsMatrix2,
  748                       MatrixElements1,
  749                       MatrixElements2,
  750                       List,
  751                       0,
  752                       0).
  753
  754matrix_select_cols(_, _, _, _, _, [], _, _) :-
  755    !.
  756matrix_select_cols(NRows, NColumnsMatrix1, NColumnsMatrix2, MatrixElements1, MatrixElements2, [_|MoreCols], NRows, ColumnIndex) :-
  757    !,
  758    NewNColumns is ColumnIndex+1,
  759    matrix_select_cols(NRows,
  760                       NColumnsMatrix1,
  761                       NColumnsMatrix2,
  762                       MatrixElements1,
  763                       MatrixElements2,
  764                       MoreCols,
  765                       0,
  766                       NewNColumns).
  767
  768matrix_select_cols(NRows, NColumnsMatrix1, NColumnsMatrix2, MatrixElements1, MatrixElements2, [Cols|MoreCols], RowIndex, ColumnIndex) :-
  769    Index1 is RowIndex*NColumnsMatrix1+Cols,
  770    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.
  785matrix_select(matrix(Type, [NRows, NColumns], MatrixElements), 1, List, Result) :-
  786    matrix_select_rows(matrix(Type, [NRows, NColumns], MatrixElements),
  787                       List,
  788                       Result).
  789
  790matrix_select(matrix(Type, [NRows, NColumns], MatrixElements), 2, List, Result) :-
  791    matrix_select_cols(matrix(Type, [NRows, NColumns], MatrixElements),
  792                       List,
  793                       Result).
 matrix_createZeroes(+Type, +Dimension, -M1)
matrix_createZeroes(M1) create a new matrix with all elements set to 0.
  797matrix_createZeroes(Type, [NRows, NColumns], M1) :-
  798    matrix_new(Type, [NRows, NColumns], M1),
  799    M1=matrix(Type, [NRows, NColumns], Elements),
  800    matrixSetZeroes(NRows, NColumns, Elements).
 matrix_read(+Dimension, -M1)
Read from keyboard values and these will be insert in M1.
  804matrix_read([NRows, NColumns], M1) :-
  805    matrix_new(doubles, [NRows, NColumns], M1),
  806    matrix_foreach(M1, \_^Y^read(Y)).
 matrixSetZeroes(+NRows, +NColumns, +Elements)
matrixSetZeroes sets to 0 all its argument matrix's elements.
  810matrixSetZeroes(NRows, NColumns, Elements) :-
  811    matrixSetAll(NRows, NColumns, Elements, 0).
 matrixSetOnes(+NRows, +NColumns, +Elements)
matrixSetOnes(M1) create a new matrix with all elements set to 1.
  815matrixSetOnes(NRows, NColumns, Elements) :-
  816    matrixSetAll(NRows, NColumns, Elements, 1).
 matrix_setAllNumbers(+NRows, +NColumns, +Elements)
matrix_setAll sets to Number given in input all the matrix's elements.
  820matrix_setAllNumbers(NRows, NColumns, Elements) :-
  821    read(Number),
  822    matrixSetAll(NRows, NColumns, Elements, Number).
 matrix_eye(+NRows, +NColumns, +Elements)
matrix_eye create a identity matrix
  826matrix_eye(NRows, NColumns, Elements) :-
  827    matrixEye(NRows, NColumns, Elements).
 matrix_equal(+M1, +M2)
matrix_equal(M1,M2) compares two matrix (M1 and M2).
  831equal_matrix(NRows, _, _, _, NRows, _) :-
  832    !.
  833
  834equal_matrix(NRows, NColumns, MatrixElements1, MatrixElements2, CurrentRow, NColumns) :-
  835    !,
  836    CurrentRow1 is CurrentRow+1,
  837    equal_matrix(NRows, NColumns, MatrixElements1, MatrixElements2, CurrentRow1, 0).
  838
  839equal_matrix(NRows, NColumns, MatrixElements1, MatrixElements2, CurrentRow, CurrentColumn) :-
  840    Index is CurrentRow*NColumns+CurrentColumn,
  841    c_load(MatrixElements1[Index], Element1),
  842    c_load(MatrixElements2[Index], Element2),
  843    CurrentColumn1 is CurrentColumn+1,
  844    Element1==Element2,
  845    equal_matrix(NRows, NColumns, MatrixElements1, MatrixElements2, CurrentRow, CurrentColumn1).
  846
  847matrix_equal(matrix(Type, [NRows, NColumns], MatrixElements1), matrix(Type, [NRows, NColumns], MatrixElements2)) :-
  848    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.
  852from_list(List, Elements) :-
  853    from_list(List, 0, Elements).
  854
  855from_list([], _, _) :-
  856    !.
  857from_list([Element|MoreElements], Index, Elements) :-
  858    c_store(Elements[Index], Element),
  859    Index1 is Index+1,
  860    from_list(MoreElements, Index1, Elements).
  861
  862matrix_from_list(List, matrix(_, [_, _], Elements)) :-
  863    from_list(List, Elements).
  864
  865%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%       START TESTS       %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  866:- begin_tests(matrix).  867
  868    test(matrix_mul1) :-
  869    matrix_new(doubles, [2, 2], [1, 2, 3, 4], Matrix1),
  870    matrix_new(doubles, [2, 2], [4, 3, 2, 1], Matrix2),
  871    matrix_new(doubles, [2, 2], [13.0, 20.0, 5.0, 8.0], Expected),
  872    matrix_mul(Matrix1, Matrix2, Matrix3),
  873    matrix_equal(Matrix3, Expected).
  874
  875    test(matrix_dims1) :-
  876    matrix_new(doubles, [3, 3], [1, 2, 3, 4, 5, 6, 7, 8, 9], Matrix1),
  877    matrix_dims(Matrix1, [3, 3]).
  878
  879    test(matrix_size1) :-
  880    matrix_new(doubles, [3, 3], [1, 2, 3, 4, 5, 6, 7, 8, 9], Matrix1),
  881    matrix_size(Matrix1, 9).
  882    
  883    test(matrix_type1) :-
  884    matrix_new(doubles, [3, 3], [1, 2, 3, 4, 5, 6, 7, 8, 9], Matrix1),
  885    matrix_type(Matrix1, doubles).
  886
  887    test(matrix_to_list1) :-
  888    matrix_new(doubles, [3, 3], [1, 2, 3, 4, 5, 6, 7, 8, 9], Matrix1),
  889    matrix_to_list(Matrix1, [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0]).
  890
  891    test(matrix_get1) :-
  892    matrix_new(doubles, [3, 3], [1, 2, 3, 4, 5, 6, 7, 8, 9], Matrix1),
  893    matrix_get(Matrix1, [0, 1], 2.0).
  894
  895    test(matrix_set1) :-
  896    matrix_new(doubles, [3, 3], [1, 2, 3, 4, 5, 6, 7, 8, 9], Matrix1),
  897    matrix_set(Matrix1, [1, 2], 5),
  898    matrix_get(Matrix1, [1, 2], 5.0).
  899
  900    test(matrix_set_all1) :-
  901    matrix_new(doubles, [2, 2], Matrix1),
  902    matrix_new(doubles, [2, 2], [7.0, 7.0, 7.0, 7.0], Expected),
  903    matrix_set_all(number(doubles, 7), Matrix1),
  904    matrix_equal(Matrix1, Expected).
  905
  906    test(matrix_add1) :-
  907    matrix_new(doubles, [3, 3], [1, 2, 3, 4, 5, 6, 7, 8, 9], Matrix1),
  908    matrix_add(Matrix1, [1, 0], 7),
  909    matrix_get(Matrix1, [1, 0], 11.0).
  910
  911    test(matrix_inc1) :-
  912    matrix_new(doubles, [3, 3], [1, 2, 3, 4, 5, 6, 7, 8, 9], Matrix1),
  913    matrix_inc(Matrix1, [0, 1], 3.0).
  914
  915    test(matrix_dec1) :-
  916    matrix_new(doubles, [3, 3], [1, 2, 3, 4, 5, 6, 7, 8, 9], Matrix1),
  917    matrix_dec(Matrix1, [0, 1], 1.0).
  918
  919    test(matrix_max1) :-
  920    matrix_new(doubles, [3, 3], [1, 2, 3, 11, 5, 6, 7, 8, 9], Matrix1),
  921    matrix_max(Matrix1, 11.0).
  922
  923    test(matrix_min1) :-
  924    matrix_new(doubles, [3, 3], [3, 2, 1, 11, 5, 6, 7, 8, 9], Matrix1),
  925    matrix_min(Matrix1, 1.0).
  926
  927    test(matrix_sum1) :-
  928    matrix_new(doubles, [3, 3], [1, 2, 3, 11, 5, 6, 7, 8, 9], Matrix1),
  929    matrix_sum(Matrix1, 52.0).
  930
  931    test(matrix_transpose1) :-
  932    matrix_new(doubles, [2, 3], [1, 2, 3, 11, 5, 6], Matrix1),
  933    matrix_new(doubles, [3, 2], [1, 11, 2, 5, 3, 6], Expected),
  934    matrix_transpose(Matrix1, Matrix2),
  935    matrix_equal(Matrix2, Expected).
  936
  937    test(matrix_maxarg1) :-
  938    matrix_new(doubles, [2, 3], [1, 2, 3, 11, 5, 6], Matrix1),
  939    matrix_maxarg(Matrix1, [1, 0]).
  940
  941    test(matrix_minarg1) :-
  942    matrix_new(doubles, [2, 3], [13, 2, 3, 11, 1, 6], Matrix1),
  943    matrix_minarg(Matrix1, [1, 1]).
  944
  945    test(matrix_agg_lines1) :-
  946    matrix_new(doubles, [2, 3], [13, 2, 3, 10, 1, 6], Matrix1),
  947    matrix_new(doubles, [2, 1], [18, 17], Expected),
  948    matrix_agg_lines(Matrix1, AggregationMatrix),
  949    matrix_equal(AggregationMatrix, Expected).
  950
  951    test(matrix_agg_cols1) :-
  952    matrix_new(doubles, [2, 3], [13, 2, 3, 10, 1, 6], Matrix1),
  953    matrix_new(doubles, [1, 3], [23, 3, 9], Expected),
  954    matrix_agg_cols(Matrix1, AggregationMatrix),
  955    matrix_equal(AggregationMatrix, Expected).
  956
  957    test(matrix_op1) :-
  958    matrix_new(doubles, [2, 2], [13, 2, 3, 10], Matrix1),
  959    matrix_new(doubles, [2, 2], [2, 8, 13, 9], Matrix2),
  960    matrix_write(Matrix1),
  961    matrix_op(Matrix1, Matrix2, +, Matrix3),
  962    matrix_new(doubles, [2, 2], [15, 10, 16, 19], Expected),
  963    matrix_equal(Matrix3, Expected).
  964
  965    test(matrix_op_to_all1) :-
  966    matrix_new(doubles, [2, 3], [13, 2, 3, 10, 1, 6], Matrix1),
  967    matrix_new(doubles, [2, 3], [16, 5, 6, 13, 4, 9], Expected),
  968    matrix_op_to_all(Matrix1, +, 3, Matrix2),
  969    matrix_equal(Matrix2, Expected).
  970
  971    test(matrix_select1) :-
  972    matrix_new(doubles, [3, 2], [13, 2, 3, 10, 1, 6], Matrix1),
  973    matrix_new(doubles, [1, 2], [13, 2], Expected),
  974    matrix_select(Matrix1, 1, [0], Matrix2),
  975    matrix_equal(Matrix2, Expected).
  976
  977    test(matrix_select2) :-
  978    matrix_new(doubles, [3, 2], [13, 2, 3, 10, 1, 6], Matrix1),
  979    matrix_new(doubles, [2, 2], [13, 2, 3, 10], Expected),
  980    matrix_select(Matrix1, 1, [0, 1], Matrix2),
  981    matrix_equal(Matrix2, Expected).
  982
  983    test(matrix_select3) :-
  984    matrix_new(doubles, [2, 3], [13, 2, 3, 10, 1, 6], Matrix1),
  985    matrix_new(doubles, [2, 2], [13, 3, 10, 6], Expected),
  986    matrix_select(Matrix1, 2, [0, 2], Matrix2),
  987    matrix_equal(Matrix2, Expected).
  988
  989    test(matrixCreateZeroes1) :-
  990    matrix_createZeroes(doubles, [3, 3], Matrix1),
  991    matrix_new(doubles, [3, 3], [0, 0, 0, 0, 0, 0, 0, 0, 0], Expected),
  992    matrix_equal(Matrix1, Expected).
  993
  994    test(matrixSetZeroes1) :-
  995    matrix_new(doubles, [3, 3], Matrix1),
  996    Matrix1=matrix(_, [3, 3], Elements),
  997    matrix_new(doubles, [3, 3], [0, 0, 0, 0, 0, 0, 0, 0, 0], Expected),
  998    matrixSetZeroes(3, 3, Elements),
  999    matrix_equal(Matrix1, Expected).
 1000
 1001    test(matrixSetOnes1) :-
 1002    matrix_new(doubles, [3, 3], Matrix1),
 1003    Matrix1=matrix(_, [3, 3], Elements),
 1004    matrix_new(doubles, [3, 3], [1, 1, 1, 1, 1, 1, 1, 1, 1], Expected),
 1005    matrixSetOnes(3, 3, Elements),
 1006    matrix_equal(Matrix1, Expected).
 1007
 1008    test(matrix_setAllNumbers1) :-
 1009    matrix_new(doubles, [3, 3], Matrix1),
 1010    Matrix1=matrix(_, [3, 3], Elements),
 1011    Number is 5.0,
 1012    matrixSetAll(3, 3, Elements, Number),
 1013    matrix_new(doubles, [3, 3], [5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0], Expected),
 1014    matrix_equal(Matrix1, Expected).
 1015
 1016    test(matrix_eye1) :-
 1017    matrix_new(doubles, [3, 3], Matrix1),
 1018    Matrix1=matrix(_, [3, 3], Elements),
 1019    matrix_eye(3, 3, Elements),
 1020    matrix_new(doubles, [3, 3], [1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0], Expected),
 1021    matrix_equal(Matrix1, Expected).
 1022
 1023    test(matrix_equal1) :-
 1024    matrix_new(doubles, [3, 3], [1, 2, 3, 4, 5, 6, 7, 8, 9], Matrix1),
 1025    matrix_new(doubles, [3, 2], [1, 3, 4, 6, 7, 9], Expected),
 1026    matrix_select(Matrix1, 2, [0, 2], Matrix2),
 1027    matrix_equal(Matrix2, Expected).
 1028
 1029    test(matrix_map1) :-
 1030    matrix_new(doubles, [2, 3], [1, 2, 3, 11, 5, 6], Matrix1),
 1031    matrix_map(\E1^E2^(E2 is E1*2),
 1032               Matrix1,
 1033               Matrix2),
 1034    matrix_new(doubles, [2, 3], [2, 4, 6, 22, 10, 12], Expected),
 1035    matrix_equal(Matrix2, Expected).
 1036
 1037    test(matrix_foreach1) :-
 1038    matrix_new(doubles, [3, 3], [1, 2, 3, 4, 5, 6, 7, 8, 9], Matrix1),
 1039    matrix_new(doubles, [3, 3], [3, 3, 3, 3, 3, 3, 3, 3, 3], Expected),
 1040    matrix_foreach(Matrix1,
 1041                   \X^Y^myset(X, Y, 3)),
 1042    matrix_equal(Matrix1, Expected).
 1043
 1044    test(matrix_foldl1) :-
 1045    matrix_new(doubles, [3, 3], [1, 2, 3, 4, 5, 6, 7, 8, 9], Matrix1),
 1046    matrix_foldl(Matrix1, findMax, 9.0).
 1047
 1048    test(matrix_from_list1) :-
 1049    matrix_new(doubles, [3, 3], [1, 2, 3, 4, 5, 6, 7, 8, 9], M1),
 1050    matrix_from_list([11, 12, 13, 14, 15, 16, 17, 18, 19], M1),
 1051    matrix_new(doubles,
 1052               [3, 3],
 1053               [11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 17.0, 18.0, 19.0],
 1054               Expected),
 1055    matrix_equal(M1, Expected).
 1056
 1057end_tests(matrix)