ffimatrix.pl

    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:- c_import("#include <atlas/clapack.h>",
   49            [liblapack_atlas],
   50            [clapack_dpotrf(int, int, int, *double, int, [int])]).
   51
   52cholesky(matrix(doubles, [Order, Order], A)) :-
   53    clapack_dpotrf(101, 122, Order, A, 3, _).

cholesky_function: Given in input M1 make the decomposition of Cholesky into M2

   57cholesky_function(Matrix1) :-
   58    Matrix1=matrix(Type, [Order, Order], _),
   59    matrix_new(Type, [Order, Order], Matrix2),
   60    matrix_sollower(Matrix1, Matrix2),
   61    cholesky(Matrix2),
   62    matrix_write(Matrix2).
   63
   64matrix_sollower(matrix(Type, [NRows, NColumns], Elements1), matrix(Type, [NRows, NColumns], Elements2)) :-
   65    sollowerMatrix(NRows, NColumns, Elements1, Elements2).
   66
   67
   68matrix_copy(matrix(Type, [NRows, NColumns], Elements1), matrix(Type, [NRows, NColumns], Elements2)) :-
   69    NElements is NRows*NColumns,
   70    arrayCopy(NElements, Elements1, Elements2).
   71
   72matrix_to_list_of_lists(M, L) :-
   73    aux0(M, 0, 0, L).
   74
   75aux0(matrix(_, [NRows, _], _), NRows, _, []) :- !.      
   76aux0(M, Row, NColumns, [[]|MoreRows]) :-
   77    M=matrix(_, [_, NColumns], _), !,
   78    Row1 is Row+1,
   79    aux0(M, Row1, 0, MoreRows).
   80aux0(M, Row, Column, [[Element|MoreElements]|MoreRows]) :-
   81    matrix_get(M, [Row, Column], Element),
   82    Column1 is Column+1,
   83    aux0(M, Row, Column1, [MoreElements|MoreRows]).
   84
   85
   86
   87
   88% https://www.dcc.fc.up.pt/~vsc/Yap/documentation.html#matrix
   89:- c_import("#include <atlas/cblas.h>",
   90            [libblas],
   91            
   92            [ cblas_dgemm(int,
   93                          int,
   94                          int,
   95                          int,
   96                          int,
   97                          int,
   98                          double,
   99                          *double,
  100                          int,
  101                          *double,
  102                          int,
  103                          double,
  104                          *double,
  105                          int)
  106            ]).
  107
  108:- (multifile user:term_expansion/2).  109:- (dynamic user:term_expansion/2).  110                
  111
  112user:term_expansion((:-import),  (:-Import)) :-
  113    current_prolog_flag(arch, Arch),
  114    atom_concat('../lib/', Arch, A1),
  115    atom_concat(A1, '/matrixNative', A2),
  116    Import=c_import("#include \"../matrixNative.h\"", [A2], [matrixSetAll(int, int, *double, double), matrixEye(int, int, *double), matrixCopy(int, int, *double, *double), arrayCopy(int, *double, *double), sollowerMatrix(int, int, *double, *double)]).
  117
  118
  119:- import.

matrix_write(+M1): Print elements of a matrix

  123matrix_write(matrix(_, [NRows, NColumns], Matrix)) :-
  124    write_matrix(Matrix, NRows, NColumns).
  125
  126write_matrix(Matrix, NRows, NColumns) :-
  127    write_matrix(Matrix, 0, 0, NRows, NColumns).
  128
  129write_matrix(_, NRows, _, NRows, _) :- !.
  130
  131write_matrix(MatrixElements, CurrentRow, NColumns, NRows, NColumns) :- !,
  132    writeln(''),
  133    CurrentRow1 is CurrentRow+1,
  134    write_matrix(MatrixElements, CurrentRow1, 0, NRows, NColumns).
  135
  136write_matrix(MatrixElements, CurrentRow, CurrentColumn, NRows, NColumns) :- !,
  137    Index is CurrentRow*NColumns+CurrentColumn,
  138    c_load(MatrixElements[Index], Element),
  139    write(Element),
  140    write(' '),
  141    CurrentColumn1 is CurrentColumn+1,
  142    write_matrix(MatrixElements, CurrentRow, CurrentColumn1, NRows, NColumns).
  143
  144%!  matrix_new(+Type,+Dimension,-M1)
  145% 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.
  146matrix_new(doubles, [NRows, NColumns], matrix(doubles, [NRows, NColumns], Matrix)) :-
  147    NRows>0,
  148    NColumns>0,
  149    Nelements is NRows*NColumns,
  150    c_alloc(Matrix, double[Nelements]).
  151
  152%!  matrix_new(+Type,+Dimension,+Elements,-M1)
  153% 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.
  154matrix_new(doubles, [NRows, NColumns], ListElements, matrix(doubles, [NRows, NColumns], Matrix)) :-
  155    NRows>0,
  156    NColumns>0,
  157    NElements is NRows*NColumns,
  158    length(ListElements, NElements),
  159    c_alloc(Matrix, double[]=ListElements).

matrix_mul(+M1, +M2, -M3): Multiplication between two matrices using "cblas_dgemm" function from "cblas.h"

  163matrix_mul(matrix(Type, [NRowsMatrix1, NColumnsMatrix1], Matrix1), matrix(Type, [NColumnsMatrix1, NColumnsMatrix2], Matrix2), matrix(Type, [NRowsMatrix1, NColumnsMatrix2], Matrix3)) :-
  164    matrix_new(Type,
  165               [NRowsMatrix1, NColumnsMatrix2],
  166               matrix(Type, [NRowsMatrix1, NColumnsMatrix2], Matrix3)),
  167    cblas_dgemm(102,
  168                111,
  169                111,
  170                NRowsMatrix1,
  171                NColumnsMatrix2,
  172                NColumnsMatrix1,
  173                1,
  174                Matrix1,
  175                NRowsMatrix1,
  176                Matrix2,
  177                NColumnsMatrix1,
  178                0,
  179                Matrix3,
  180                NRowsMatrix1).

matrix_dims(+M1, -Dimension): Calculates the total number of elements in the matrix. The value is contained in Dimension.

  184matrix_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.

  188matrix_size(matrix(_, [NRows, NColumns], _), Nelements) :-
  189    Nelements is NRows*NColumns.

matrix_type(+M1, -Type): Function that return the Type of the matrix passed.

  193matrix_type(matrix(Type, _, _), Type).

matrix_to_list(+M1, -List): Function that transforms the matrix passed in input in a list.

  197to_list(_, NRows, _, NRows, _, List) :- !,
  198    List=([]).
  199
  200to_list(MatrixElements, CurrentRow, NColumns, Nrows, NColumns, List) :- !,
  201    CurrentRow1 is CurrentRow+1,
  202    to_list(MatrixElements, CurrentRow1, 0, Nrows, NColumns, List).
  203
  204to_list(MatrixElements, CurrentRow, CurrentColumn, NRows, NColumns, [Element|MoreElements]) :-
  205    Index is CurrentRow*NColumns+CurrentColumn,
  206    c_load(MatrixElements[Index], Element),
  207    CurrentColumn1 is CurrentColumn+1,
  208    to_list(MatrixElements, CurrentRow, CurrentColumn1, NRows, NColumns, MoreElements).
  209
  210matrix_to_list(matrix(_, [NRows, NColumns], MatrixElements), List) :-
  211    to_list(MatrixElements, 0, 0, NRows, NColumns, List).
  212
  213%!  matrix_get(+M1,+Position,-Element)
  214% matrix_get(M,[Row,Column],E) unifies E with the element of M at position [Row, Column] 
  215matrix_get(matrix(_, [_, NColumns], MatrixElements), [Row, Column], Element) :-
  216    Index is Row*NColumns+Column,
  217    c_load(MatrixElements[Index], Element).
  218
  219%!  matrix_set(+Type,+Position,-Element)
  220% matrix_set(M,[Row,Column],E) update the element of M at position [Row, Column] with the element specified in Element.
  221matrix_set(matrix(_, [_, NColumns], MatrixElements), [Row, Column], Element) :-
  222    Index is Row*NColumns+Column,
  223    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.

  227matrix_map(Predicate, matrix(Type, [NRows, NColumns], MatrixElements1), matrix(Type, [NRows, NColumns], MatrixElements2)) :-
  228    matrix_new(Type,
  229               [NRows, NColumns],
  230               matrix(_, _, MatrixElements2)),
  231    matrix_map(Predicate, NRows, NColumns, MatrixElements1, MatrixElements2, 0, 0).
  232
  233matrix_map(_, NRows, _, _, _, NRows, _) :- !.
  234
  235matrix_map(Predicate, NRows, NColumns, MatrixElements1, MatrixElements2, CurrentRow, NColumns) :- !,
  236    CurrentRow1 is CurrentRow+1,
  237    matrix_map(Predicate, NRows, NColumns, MatrixElements1, MatrixElements2, CurrentRow1, 0).
  238
  239matrix_map(Predicate, NRows, NColumns, MatrixElements1, MatrixElements2, CurrentRow, CurrentColumn) :-
  240    Index is CurrentRow*NColumns+CurrentColumn,
  241    c_load(MatrixElements1[Index], E1),
  242    call(Predicate, E1, E2),
  243    c_store(MatrixElements2[Index], E2),
  244    CurrentColumn1 is CurrentColumn+1,
  245    matrix_map(Predicate,
  246               NRows,
  247               NColumns,
  248               MatrixElements1,
  249               MatrixElements2,
  250               CurrentRow,
  251               CurrentColumn1).
  252
  253%%%%%%%%%%%%%%%%%%%%%%%%%%%%       TWO FUNCTIONS USED BY MATRIX FOREACH       %%%%%%%%%%%%%%%%%%%%%%%%%%%%
  254
  255% Unify the element with the number specified.
  256myset(_, Number, Number).
  257
  258myset(_, 0).
  259% This function find the max number between two number.
  260findMax(FirstNumber, SecondNumber, FirstNumber) :-
  261    FirstNumber>=SecondNumber, !.
  262
  263findMax(_, SecondNumber, SecondNumber).
  264
  265%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

matrix_foreach(+M1, +Predicate): For each element E of M1, enumerated by rows, calls Predicate(+E,-F) and replaces E with F.

  269matrix_foreach(matrix(_, [NRows, NColumns], MatrixElements), Predicate) :-
  270    matrix_foreach(Predicate, NRows, NColumns, MatrixElements, 0, 0).
  271
  272matrix_foreach(_, NRows, _, _, NRows, _) :- !.
  273
  274matrix_foreach(Predicate, NRows, NColumns, MatrixElements, CurrentRow, NColumns) :- !,
  275    CurrentRow1 is CurrentRow+1,
  276    matrix_foreach(Predicate, NRows, NColumns, MatrixElements, CurrentRow1, 0).
  277
  278matrix_foreach(Predicate, NRows, NColumns, MatrixElements, CurrentRow, CurrentColumn) :-
  279    Index is CurrentRow*NColumns+CurrentColumn,
  280    c_load(MatrixElements[Index], E),
  281    call(Predicate, E, E1),
  282    c_store(MatrixElements[Index], E1),
  283    CurrentColumn1 is CurrentColumn+1,
  284    matrix_foreach(Predicate,
  285                   NRows,
  286                   NColumns,
  287                   MatrixElements,
  288                   CurrentRow,
  289                   CurrentColumn1).
  290
  291%!  matrix_foldl(+M1,+Predicate,-Max)
  292% Fold a matrix, using arguments of the matrix as left argument.
  293matrix_foldl(matrix(_, [NRows, NColumns], MatrixElements), Predicate, Max) :-
  294    c_load(MatrixElements[0], TmpMax),
  295    matrix_start_foldl(Predicate,
  296                       NRows,
  297                       NColumns,
  298                       MatrixElements,
  299                       TmpMax,
  300                       Max).
  301
  302matrix_start_foldl(_, 1, 1, _, Max, Max) :- !.
  303matrix_start_foldl(Predicate, NRows, NColumns, MatrixElements, Tmp, Max) :-
  304    NColumns>1, !,
  305    matrix_foldl(Predicate,
  306                 NRows,
  307                 NColumns,
  308                 MatrixElements,
  309                 0,
  310                 1,
  311                 Tmp,
  312                 Max).
  313
  314matrix_start_foldl(Predicate, NRows, NColumns, MatrixElements, Tmp, Max) :-
  315    matrix_foldl(Predicate,
  316                 NRows,
  317                 NColumns,
  318                 MatrixElements,
  319                 1,
  320                 0,
  321                 Tmp,
  322                 Max).
  323
  324matrix_foldl(_, NRows, _, _, NRows, _, Max, Max) :- !.
  325
  326matrix_foldl(Predicate, NRows, NColumns, MatrixElements, CurrentRow, NColumns, TmpMax, Max) :- !,
  327    CurrentRow1 is CurrentRow+1,
  328    matrix_foldl(Predicate,
  329                 NRows,
  330                 NColumns,
  331                 MatrixElements,
  332                 CurrentRow1,
  333                 0,
  334                 TmpMax,
  335                 Max).
  336
  337matrix_foldl(Predicate, NRows, NColumns, MatrixElements, CurrentRow, CurrentColumn, TmpMax, Max) :-
  338    Index is CurrentRow*NColumns+CurrentColumn,
  339    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.

  353set_all(_, NRows, _, _, NRows, _) :- !.
  354
  355set_all(Number, NRows, NColumns, MatrixElements, CurrentRow, NColumns) :- !,
  356    CurrentRow1 is CurrentRow+1,
  357    set_all(Number, NRows, NColumns, MatrixElements, CurrentRow1, 0).
  358
  359set_all(Number, NRows, NColumns, MatrixElements, CurrentRow, CurrentColumn) :-
  360    Index is CurrentRow*NColumns+CurrentColumn,
  361    c_store(MatrixElements[Index], Number),
  362    CurrentColumn1 is CurrentColumn+1,
  363    set_all(Number, NRows, NColumns, MatrixElements, CurrentRow, CurrentColumn1).
  364
  365matrix_set_all(number(Type, Number), matrix(Type, [NRows, NColumns], Element)) :-
  366    set_all(Number, NRows, NColumns, Element, 0, 0).
  367
  368%!  matrix_add(+M1,Position,+Number)
  369% matrix_add(M,[Row,Column],Number) add Number to the element of matrix M in position [Row,Column].
  370matrix_add(matrix(_, [_, NColumns], MatrixElements), [Row, Column], Number) :-
  371    Index is Row*NColumns+Column,
  372    c_load(MatrixElements[Index], Tmp),
  373    NewElement is Tmp+Number,
  374    c_store(MatrixElements[Index], NewElement).
  375
  376%!  matrix_inc(+M1,+Position)
  377% matrix_inc(M,[Row,Column]) add 1 (one) to the element of matrix M in position [Row,Column].
  378matrix_inc(matrix(_, [_, NColumns], MatrixElements), [Row, Column]) :-
  379    Index is Row*NColumns+Column,
  380    c_load(MatrixElements[Index], Tmp),
  381    E is Tmp+1,
  382    c_store(MatrixElements[Index], E).
  383
  384%!  matrix_inc(+M1,+Position,-NewElement)
  385% 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.
  386matrix_inc(matrix(_, [_, NColumns], MatrixElements), [Row, Column], NewElement) :-
  387    Index is Row*NColumns+Column,
  388    c_load(MatrixElements[Index], Tmp),
  389    NewElement is Tmp+1,
  390    c_store(MatrixElements[Index], NewElement).
  391
  392%!  matrix_dec(+M1,+Position)
  393% matrix_dec(M,[Row,Column]) decrement 1 (one) to the element of matrix M in position [Row,Column].
  394matrix_dec(matrix(_, [_, NColumns], MatrixElements), [Row, Column]) :-
  395    Index is Row*NColumns+Column,
  396    c_load(MatrixElements[Index], Tmp),
  397    NewElement is Tmp-1,
  398    c_store(MatrixElements[Index], NewElement).
  399
  400%!  matrix_dec(+M1,+Position,-NewElement)
  401% 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.
  402matrix_dec(matrix(_, [_, NColumns], MatrixElements), [Row, Column], NewElement) :-
  403    Index is Row*NColumns+Column,
  404    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.

  410max_matrix(NRows, _, _, NRows, _, Max, Max) :- !.
  411
  412max_matrix(NRows, NColumns, MatrixElements, CurrentRow, NColumns, TmpMax, Max) :- !,
  413    CurrentRow1 is CurrentRow+1,
  414    max_matrix(NRows, NColumns, MatrixElements, CurrentRow1, 0, TmpMax, Max).
  415
  416max_matrix(NRows, NColumns, MatrixElements, CurrentRow, CurrentColumn, TmpMax, Max) :-
  417    Index is CurrentRow*NColumns+CurrentColumn,
  418    c_load(MatrixElements[Index], E),
  419    CurrentColumn1 is CurrentColumn+1,
  420    (   E>TmpMax
  421    ->  TmpMax1 is E
  422    ;   TmpMax1 is TmpMax
  423    ),
  424    max_matrix(NRows,
  425               NColumns,
  426               MatrixElements,
  427               CurrentRow,
  428               CurrentColumn1,
  429               TmpMax1,
  430               Max).
  431
  432matrix_max(matrix(_, [NRows, NColumns], MatrixElements), Max) :-
  433    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.

  438matrix_min(Matrix, Min) :-
  439    matrix_foldl(Matrix,
  440                 \X^Y^Z^(X=<Y->Z=X;Z=Y),
  441                 Min).

matrix_sum(+M1, -Sum): matrix_sum(M,Sum) unify Sum with the sum of the elements in M.

  445matrix_sum(Matrix, Sum) :-
  446    matrix_foldl(Matrix,
  447                 \CurrentSum^Element^NewSum^(NewSum is CurrentSum+Element),
  448                 Sum).

matrix_transpose(+M1, -M2): matrix_transpose(M1,M2) transpose matrix M1 to M2.

  452transpose_matrix(NRows, _, _, _, NRows, _) :- !.
  453
  454transpose_matrix(NRows, NColumns, MatrixElements1, MatrixElements2, CurrentRow, NColumns) :- !,
  455    CurrentRow1 is CurrentRow+1,
  456    transpose_matrix(NRows, NColumns, MatrixElements1, MatrixElements2, CurrentRow1, 0).
  457
  458transpose_matrix(NRows, NColumns, MatrixElements1, MatrixElements2, CurrentRow, CurrentColumn) :-
  459    Index is CurrentRow*NColumns+CurrentColumn,
  460    IndexMatrix2 is CurrentColumn*NRows+CurrentRow,
  461    c_load(MatrixElements1[Index], E),
  462    c_store(MatrixElements2[IndexMatrix2], E),
  463    CurrentColumn1 is CurrentColumn+1,
  464    transpose_matrix(NRows,
  465                     NColumns,
  466                     MatrixElements1,
  467                     MatrixElements2,
  468                     CurrentRow,
  469                     CurrentColumn1).
  470
  471
  472matrix_transpose(matrix(Type, [NRows, NColumns], MatrixElements1), Matrix2) :-
  473    matrix_new(Type, [NColumns, NRows], Matrix2),
  474    Matrix2=matrix(_, _, MatrixElements2),
  475    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.

  479maxarg_matrix(NRows, _, _, NRows, _, _, [NRowsMax, NColumnsMax], [NRowsMax, NColumnsMax]) :- !.
  480
  481maxarg_matrix(NRows, NColumns, MatrixElements, CurrentRow, NColumns, TmpMax, [NRowTmp, NColumnTmp], [NRowMax, NColomnMax]) :- !,
  482    CurrentRow1 is CurrentRow+1,
  483    maxarg_matrix(NRows,
  484                  NColumns,
  485                  MatrixElements,
  486                  CurrentRow1,
  487                  0,
  488                  TmpMax,
  489                  [NRowTmp, NColumnTmp],
  490                  [NRowMax, NColomnMax]).
  491
  492                maxarg_matrix(NRows, NColumns, MatrixElements, CurrentRow, CurrentColumn, TmpMax, [NRowTmp, NColumnTmp], [NRowMax, NColomnMax]) :-
  493    Index is CurrentRow*NColumns+CurrentColumn,
  494    c_load(MatrixElements[Index], Element),
  495    CurrentColumn1 is CurrentColumn+1,
  496    (   Element>=TmpMax
  497    ->  ( TmpMax1 is Element,
  498          NRowTmp1 is CurrentRow,
  499          NColumnTmp1 is CurrentColumn
  500        ),
  501        maxarg_matrix(NRows,
  502                      NColumns,
  503                      MatrixElements,
  504                      CurrentRow,
  505                      CurrentColumn1,
  506                      TmpMax1,
  507                      [NRowTmp1, NColumnTmp1],
  508                      [NRowMax, NColomnMax])
  509    ;   TmpMax1 is TmpMax,
  510        NRowTmp1 is NRowTmp,
  511        NColumnTmp1 is NColumnTmp,
  512        maxarg_matrix(NRows,
  513                      NColumns,
  514                      MatrixElements,
  515                      CurrentRow,
  516                      CurrentColumn1,
  517                      TmpMax1,
  518                      [NRowTmp1, NColumnTmp1],
  519                      [NRowMax, NColomnMax])
  520    ).
  521
  522matrix_maxarg(matrix(_, [NRows, NColumns], MatrixElements), [NRowMax, NColomnMax]) :-
  523    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.

  535minarg_matrix(NRows, _, _, NRows, _, _, [NRowTmp, NColumnTmp], [NRowTmp, NColumnTmp]) :- !.
  536
  537minarg_matrix(NRows, NColumns, MatrixElements, CurrentRow, NColumns, TmpMax, [NRowTmp, NColumnTmp], [NRowMax, NColumnMax]) :- !,
  538    CurrentRow1 is CurrentRow+1,
  539    minarg_matrix(NRows,
  540                  NColumns,
  541                  MatrixElements,
  542                  CurrentRow1,
  543                  0,
  544                  TmpMax,
  545                  [NRowTmp, NColumnTmp],
  546                  [NRowMax, NColumnMax]).
  547
  548                minarg_matrix(NRows, NColumns, MatrixElements, CurrentRow, CurrentColumn, TmpMax, [NRowTmp, NColumnTmp], [NRowMax, NColumnMax]) :-
  549    Index is CurrentRow*NColumns+CurrentColumn,
  550    c_load(MatrixElements[Index], Element),
  551    CurrentColumn1 is CurrentColumn+1,
  552    (   Element=<TmpMax
  553    ->  TmpMax1 is Element,
  554        NRowTmp1 is CurrentRow,
  555        NColumnTmp1 is CurrentColumn
  556    ;   TmpMax1 is TmpMax,
  557        NRowTmp1 is NRowTmp,
  558        NColumnTmp1 is NColumnTmp
  559    ),
  560    minarg_matrix(NRows,
  561                  NColumns,
  562                  MatrixElements,
  563                  CurrentRow,
  564                  CurrentColumn1,
  565                  TmpMax1,
  566                  [NRowTmp1, NColumnTmp1],
  567                  [NRowMax, NColumnMax]).
  568
  569matrix_minarg(matrix(_, [NRows, NColumns], MatrixElements), [NRowMax, NColumnMax]) :-
  570    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.

  582agg_lines_matrix(NRows, _, _, _, NRows, _, _) :- !.
  583
  584agg_lines_matrix(NRows, NColumns, MatrixElements1, MatrixElements2, CurrentRow, NColumns, SumRow) :- !,
  585    c_store(MatrixElements2[CurrentRow], SumRow),
  586    SumRow1 is 0,
  587    CurrentRow1 is CurrentRow+1,
  588    agg_lines_matrix(NRows,
  589                     NColumns,
  590                     MatrixElements1,
  591                     MatrixElements2,
  592                     CurrentRow1,
  593                     0,
  594                     SumRow1).
  595
  596agg_lines_matrix(NRows, NColumns, MatrixElements1, MatrixElements2, CurrentRow, CurrentColumn, SumRow) :-
  597    Index is CurrentRow*NColumns+CurrentColumn,
  598    c_load(MatrixElements1[Index], Element),
  599    SumRow1 is SumRow+Element,
  600    CurrentColumn1 is CurrentColumn+1,
  601    agg_lines_matrix(NRows,
  602                     NColumns,
  603                     MatrixElements1,
  604                     MatrixElements2,
  605                     CurrentRow,
  606                     CurrentColumn1,
  607                     SumRow1).
  608
  609
  610matrix_agg_lines(matrix(Type, [NRows, NColumns], MatrixElements1), AggregationMatrix) :-
  611    matrix_new(Type, [NRows, 1], AggregationMatrix),
  612    AggregationMatrix=matrix(_, _, MatrixElements2),
  613    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.

  617agg_cols_matrix(_, NColumns, _, _, _, NColumns, _) :- !.
  618
  619agg_cols_matrix(NRows, NColumns, MatrixElements1, MatrixElements2, NRows, CurrentColumn, SumRow) :- !,
  620    c_store(MatrixElements2[CurrentColumn], SumRow),
  621    SumRow1 is 0,
  622    CurrentColumn1 is CurrentColumn+1,
  623    agg_cols_matrix(NRows,
  624                    NColumns,
  625                    MatrixElements1,
  626                    MatrixElements2,
  627                    0,
  628                    CurrentColumn1,
  629                    SumRow1).
  630
  631agg_cols_matrix(NRows, NColumns, MatrixElements1, MatrixElements2, CurrentRow, CurrentColumn, SumRow) :-
  632    Index is CurrentRow*NColumns+CurrentColumn,
  633    c_load(MatrixElements1[Index], Element),
  634    SumRow1 is SumRow+Element,
  635    CurrentRow1 is CurrentRow+1,
  636    agg_cols_matrix(NRows,
  637                    NColumns,
  638                    MatrixElements1,
  639                    MatrixElements2,
  640                    CurrentRow1,
  641                    CurrentColumn,
  642                    SumRow1).
  643
  644
  645matrix_agg_cols(matrix(Type, [NRows, NColumns], MatrixElements), AggregationMatrix) :-
  646    matrix_new(Type, [1, NColumns], AggregationMatrix),
  647    AggregationMatrix=matrix(_, _, MatrixElements2),
  648    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.

  652matrix_map(Predicate, matrix(Type, [NRows, NColumns], MatrixElements1), matrix(Type, [NRows, NColumns], MatrixElements2), matrix(Type, [NRows, NColumns], MatrixElements3)) :-
  653    matrix_new(Type,
  654               [NRows, NColumns],
  655               matrix(_, _, MatrixElements3)),
  656    matrix_map(Predicate,
  657               NRows,
  658               NColumns,
  659               MatrixElements1,
  660               MatrixElements2,
  661               MatrixElements3,
  662               0,
  663               0).
  664
  665matrix_map(_, NRows, _, _, _, _, NRows, _) :- !.
  666
  667matrix_map(Predicate, NRows, NColumns, MatrixElements1, MatrixElements2, MatrixElements3, CurrentRow, NColumns) :- !,
  668    CurrentRow1 is CurrentRow+1,
  669    matrix_map(Predicate,
  670               NRows,
  671               NColumns,
  672               MatrixElements1,
  673               MatrixElements2,
  674               MatrixElements3,
  675               CurrentRow1,
  676               0).
  677
  678matrix_map(Predicate, NRows, NColumns, MatrixElements1, MatrixElements2, MatrixElements3, CurrentRow, CurrentColumn) :-
  679    Index is CurrentRow*NColumns+CurrentColumn,
  680    c_load(MatrixElements1[Index], Element1),
  681    c_load(MatrixElements2[Index], Element2),
  682    call(Predicate, Element1, Element2, Element3),
  683    c_store(MatrixElements3[Index], Element3),
  684    CurrentColumn1 is CurrentColumn+1,
  685    matrix_map(Predicate,
  686               NRows,
  687               NColumns,
  688               MatrixElements1,
  689               MatrixElements2,
  690               MatrixElements3,
  691               CurrentRow,
  692               CurrentColumn1).
  693
  694matrix_op(matrix(Type, [NRows, NColumns], Element1), matrix(Type, [NRows, NColumns], Element2), Operator, Matrix3) :-
  695    LambdaExpression= \E1^E2^E3^(Predicate=..[Operator, E1, E2], E3 is Predicate),
  696    matrix_map(LambdaExpression,
  697               matrix(Type, [NRows, NColumns], Element1),
  698               matrix(Type, [NRows, NColumns], Element2),
  699               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.

  703matrix_op_to_all(matrix(Type, [NRows, NColumns], MatrixElements), Operator, Operand, Matrix2) :-
  704    LambdaExpression= \E1^E2^(Predicate=..[Operator, E1, Operand], E2 is Predicate),
  705    matrix_map(LambdaExpression,
  706               matrix(Type, [NRows, NColumns], MatrixElements),
  707               Matrix2).
  708
  709% function that verifies that the inserted rows or columns are valid
  710function_call(_, []) :- !.
  711function_call(P, [H|T]) :-
  712    call(P, H),
  713    function_call(P, T).
  714
  715function_call_lesser(L, M) :-
  716    P= \H^(H<M),
  717    function_call(P, L).
  718
  719% matrix_select_rows(M1,RowList,M2) select a list of row (RowList) from M1 and this rows will be copied in M2.
  720matrix_select_rows(matrix(Type, [NRows, NColumns], MatrixElements1), RowList, matrix(Type, [NRows2, NColumns], MatrixElements2)) :-
  721    function_call_lesser(RowList, NRows),
  722    length(RowList, NRows2),
  723    matrix_new(Type,
  724               [NRows2, NColumns],
  725               matrix(_, _, MatrixElements2)),
  726    matrix_select_rows(NColumns, MatrixElements1, MatrixElements2, RowList, 0, 0).
  727
  728matrix_select_rows(_, _, _, [], _, _) :- !.
  729matrix_select_rows(NColumns, MatrixElements1, MatrixElements2, [_|MoreRows], RowIndex, NColumns) :- !,
  730    RowIndex1 is RowIndex+1,
  731    matrix_select_rows(NColumns, MatrixElements1, MatrixElements2, MoreRows, RowIndex1, 0).
  732
  733matrix_select_rows(NColumns, MatrixElements1, MatrixElements2, [Row|MoreRows], RowIndex, ColumnIndex) :-
  734    Index1 is Row*NColumns+ColumnIndex,
  735    c_load(MatrixElements1[Index1], E),
  736    Index2 is RowIndex*NColumns+ColumnIndex,
  737    c_store(MatrixElements2[Index2], E),
  738    ColumnIndex1 is ColumnIndex+1,
  739    matrix_select_rows(NColumns,
  740                       MatrixElements1,
  741                       MatrixElements2,
  742                       [Row|MoreRows],
  743                       RowIndex,
  744                       ColumnIndex1).
  745
  746% matrix_select_cols(M1,ColList,M2) select a list of columns (ColList) from M1 and this columns will be copied in M2.
  747matrix_select_cols(matrix(Type, [NRows, NColumnsMatrix1], MatrixElements1), List, matrix(Type, [NRows, NColumnsMatrix2], MatrixElements2)) :-
  748    function_call_lesser(List, NColumnsMatrix1),
  749    length(List, NColumnsMatrix2),
  750    matrix_new(Type,
  751               [NRows, NColumnsMatrix2],
  752               matrix(_, _, MatrixElements2)),
  753    matrix_select_cols(NRows,
  754                       NColumnsMatrix1,
  755                       NColumnsMatrix2,
  756                       MatrixElements1,
  757                       MatrixElements2,
  758                       List,
  759                       0,
  760                       0).
  761
  762matrix_select_cols(_, _, _, _, _, [], _, _) :- !.
  763matrix_select_cols(NRows, NColumnsMatrix1, NColumnsMatrix2, MatrixElements1, MatrixElements2, [_|MoreCols], NRows, ColumnIndex) :- !,
  764    NewNColumns is ColumnIndex+1,
  765    matrix_select_cols(NRows,
  766                       NColumnsMatrix1,
  767                       NColumnsMatrix2,
  768                       MatrixElements1,
  769                       MatrixElements2,
  770                       MoreCols,
  771                       0,
  772                       NewNColumns).
  773
  774matrix_select_cols(NRows, NColumnsMatrix1, NColumnsMatrix2, MatrixElements1, MatrixElements2, [Cols|MoreCols], RowIndex, ColumnIndex) :-
  775    Index1 is RowIndex*NColumnsMatrix1+Cols,
  776    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.

  791matrix_select(matrix(Type, [NRows, NColumns], MatrixElements), 1, List, Result) :-
  792    matrix_select_rows(matrix(Type, [NRows, NColumns], MatrixElements),
  793                       List,
  794                       Result).
  795
  796matrix_select(matrix(Type, [NRows, NColumns], MatrixElements), 2, List, Result) :-
  797    matrix_select_cols(matrix(Type, [NRows, NColumns], MatrixElements),
  798                       List,
  799                       Result).

matrix_createZeroes(+Type, +Dimension, -M1): matrix_createZeroes(M1) create a new matrix with all elements set to 0.

  803matrix_createZeroes(Type, [NRows, NColumns], M1) :-
  804    matrix_new(Type, [NRows, NColumns], M1),
  805    M1=matrix(Type, [NRows, NColumns], Elements),
  806    matrixSetZeroes(NRows, NColumns, Elements).

matrix_read(+Dimension, -M1): Read from keyboard values and these will be insert in M1.

  810matrix_read([NRows, NColumns], M1) :-
  811    matrix_new(doubles, [NRows, NColumns], M1),
  812    matrix_foreach(M1, \_^Y^read(Y)).

matrixSetZeroes(+NRows, +NColumns, +Elements): matrixSetZeroes sets to 0 all its argument matrix's elements.

  816matrixSetZeroes(NRows, NColumns, Elements) :-
  817    matrixSetAll(NRows, NColumns, Elements, 0).

matrixSetOnes(+NRows, +NColumns, +Elements): matrixSetOnes(M1) create a new matrix with all elements set to 1.

  821matrixSetOnes(NRows, NColumns, Elements) :-
  822    matrixSetAll(NRows, NColumns, Elements, 1).

matrix_setAllNumbers(+NRows, +NColumns, +Elements): matrix_setAll sets to Number given in input all the matrix's elements.

  826matrix_setAllNumbers(NRows, NColumns, Elements) :-
  827    read(Number),
  828    matrixSetAll(NRows, NColumns, Elements, Number).

matrix_eye(+NRows, +NColumns, +Elements): matrix_eye create a identity matrix

  832matrix_eye(NRows, NColumns, Elements) :-
  833    matrixEye(NRows, NColumns, Elements).

matrix_equal(+M1, +M2): matrix_equal(M1,M2) compares two matrix (M1 and M2).

  837equal_matrix(NRows, _, _, _, NRows, _) :- !.
  838
  839equal_matrix(NRows, NColumns, MatrixElements1, MatrixElements2, CurrentRow, NColumns) :- !,
  840    CurrentRow1 is CurrentRow+1,
  841    equal_matrix(NRows, NColumns, MatrixElements1, MatrixElements2, CurrentRow1, 0).
  842
  843equal_matrix(NRows, NColumns, MatrixElements1, MatrixElements2, CurrentRow, CurrentColumn) :-
  844    Index is CurrentRow*NColumns+CurrentColumn,
  845    c_load(MatrixElements1[Index], Element1),
  846    c_load(MatrixElements2[Index], Element2),
  847    CurrentColumn1 is CurrentColumn+1,
  848    Element1==Element2,
  849    equal_matrix(NRows, NColumns, MatrixElements1, MatrixElements2, CurrentRow, CurrentColumn1).
  850
  851matrix_equal(matrix(Type, [NRows, NColumns], MatrixElements1), matrix(Type, [NRows, NColumns], MatrixElements2)) :-
  852    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.

  856from_list(List, Elements) :-
  857    from_list(List, 0, Elements).
  858
  859from_list([], _, _) :- !.
  860from_list([Element|MoreElements], Index, Elements) :-
  861    c_store(Elements[Index], Element),
  862    Index1 is Index+1,
  863    from_list(MoreElements, Index1, Elements).
  864
  865matrix_from_list(List, matrix(_, [_, _], Elements)) :-
  866    from_list(List, Elements).
  867
  868example(Elements) :-
  869    matrix_new(doubles, [3, 3], Matrix1),
  870    Matrix1=matrix(_, [3, 3], Elements),
  871    matrix_eye(3, 3, Elements).
  872
  873
  874matrix_random(NRows, NColumns, Matrix) :-
  875    matrix_new(doubles, [NRows, NColumns], Matrix),
  876    matrix_foreach(Matrix,
  877                   \X^Y^(random(E), myset(X, Y, E))).
  878
  879aux1 :-
  880    matrix_random(2, 2, Matrix1),
  881    matrix_transpose(Matrix1, Matrix2),
  882    matrix_mul(Matrix1, Matrix2, Matrix3),
  883    matrix_write(Matrix3).
  884
  885%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%       START TESTS       %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  886:- begin_tests(matrix).  887
  888    test(matrix_mul1) :-
  889    matrix_new(doubles, [2, 2], [1, 2, 3, 4], Matrix1),
  890    matrix_new(doubles, [2, 2], [4, 3, 2, 1], Matrix2),
  891    matrix_new(doubles, [2, 2], [13.0, 20.0, 5.0, 8.0], Expected),
  892    matrix_mul(Matrix1, Matrix2, Matrix3),
  893    matrix_equal(Matrix3, Expected).
  894
  895    test(matrix_dims1) :-
  896    matrix_new(doubles, [3, 3], [1, 2, 3, 4, 5, 6, 7, 8, 9], Matrix1),
  897    matrix_dims(Matrix1, [3, 3]).
  898
  899    test(matrix_size1) :-
  900    matrix_new(doubles, [3, 3], [1, 2, 3, 4, 5, 6, 7, 8, 9], Matrix1),
  901    matrix_size(Matrix1, 9).
  902    
  903    test(matrix_type1) :-
  904    matrix_new(doubles, [3, 3], [1, 2, 3, 4, 5, 6, 7, 8, 9], Matrix1),
  905    matrix_type(Matrix1, doubles).
  906
  907    test(matrix_to_list1) :-
  908    matrix_new(doubles, [3, 3], [1, 2, 3, 4, 5, 6, 7, 8, 9], Matrix1),
  909    matrix_to_list(Matrix1, [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0]).
  910
  911    test(matrix_get1) :-
  912    matrix_new(doubles, [3, 3], [1, 2, 3, 4, 5, 6, 7, 8, 9], Matrix1),
  913    matrix_get(Matrix1, [0, 1], 2.0).
  914
  915    test(matrix_set1) :-
  916    matrix_new(doubles, [3, 3], [1, 2, 3, 4, 5, 6, 7, 8, 9], Matrix1),
  917    matrix_set(Matrix1, [1, 2], 5),
  918    matrix_get(Matrix1, [1, 2], 5.0).
  919
  920    test(matrix_set_all1) :-
  921    matrix_new(doubles, [2, 2], Matrix1),
  922    matrix_new(doubles, [2, 2], [7.0, 7.0, 7.0, 7.0], Expected),
  923    matrix_set_all(number(doubles, 7), Matrix1),
  924    matrix_equal(Matrix1, Expected).
  925
  926    test(matrix_add1) :-
  927    matrix_new(doubles, [3, 3], [1, 2, 3, 4, 5, 6, 7, 8, 9], Matrix1),
  928    matrix_add(Matrix1, [1, 0], 7),
  929    matrix_get(Matrix1, [1, 0], 11.0).
  930
  931    test(matrix_inc1) :-
  932    matrix_new(doubles, [3, 3], [1, 2, 3, 4, 5, 6, 7, 8, 9], Matrix1),
  933    matrix_inc(Matrix1, [0, 1], 3.0).
  934
  935    test(matrix_dec1) :-
  936    matrix_new(doubles, [3, 3], [1, 2, 3, 4, 5, 6, 7, 8, 9], Matrix1),
  937    matrix_dec(Matrix1, [0, 1], 1.0).
  938
  939    test(matrix_max1) :-
  940    matrix_new(doubles, [3, 3], [1, 2, 3, 11, 5, 6, 7, 8, 9], Matrix1),
  941    matrix_max(Matrix1, 11.0).
  942
  943    test(matrix_min1) :-
  944    matrix_new(doubles, [3, 3], [3, 2, 1, 11, 5, 6, 7, 8, 9], Matrix1),
  945    matrix_min(Matrix1, 1.0).
  946
  947    test(matrix_sum1) :-
  948    matrix_new(doubles, [3, 3], [1, 2, 3, 11, 5, 6, 7, 8, 9], Matrix1),
  949    matrix_sum(Matrix1, 52.0).
  950
  951    test(matrix_transpose1) :-
  952    matrix_new(doubles, [2, 3], [1, 2, 3, 11, 5, 6], Matrix1),
  953    matrix_new(doubles, [3, 2], [1, 11, 2, 5, 3, 6], Expected),
  954    matrix_transpose(Matrix1, Matrix2),
  955    matrix_equal(Matrix2, Expected).
  956
  957    test(matrix_maxarg1) :-
  958    matrix_new(doubles, [2, 3], [1, 2, 3, 11, 5, 6], Matrix1),
  959    matrix_maxarg(Matrix1, [1, 0]).
  960
  961    test(matrix_minarg1) :-
  962    matrix_new(doubles, [2, 3], [13, 2, 3, 11, 1, 6], Matrix1),
  963    matrix_minarg(Matrix1, [1, 1]).
  964
  965    test(matrix_agg_lines1) :-
  966    matrix_new(doubles, [2, 3], [13, 2, 3, 10, 1, 6], Matrix1),
  967    matrix_new(doubles, [2, 1], [18, 17], Expected),
  968    matrix_agg_lines(Matrix1, AggregationMatrix),
  969    matrix_equal(AggregationMatrix, Expected).
  970
  971    test(matrix_agg_cols1) :-
  972    matrix_new(doubles, [2, 3], [13, 2, 3, 10, 1, 6], Matrix1),
  973    matrix_new(doubles, [1, 3], [23, 3, 9], Expected),
  974    matrix_agg_cols(Matrix1, AggregationMatrix),
  975    matrix_equal(AggregationMatrix, Expected).
  976
  977    test(matrix_op1) :-
  978    matrix_new(doubles, [2, 2], [13, 2, 3, 10], Matrix1),
  979    matrix_new(doubles, [2, 2], [2, 8, 13, 9], Matrix2),
  980    matrix_write(Matrix1),
  981    matrix_op(Matrix1, Matrix2, +, Matrix3),
  982    matrix_new(doubles, [2, 2], [15, 10, 16, 19], Expected),
  983    matrix_equal(Matrix3, Expected).
  984
  985    test(matrix_op_to_all1) :-
  986    matrix_new(doubles, [2, 3], [13, 2, 3, 10, 1, 6], Matrix1),
  987    matrix_new(doubles, [2, 3], [16, 5, 6, 13, 4, 9], Expected),
  988    matrix_op_to_all(Matrix1, +, 3, Matrix2),
  989    matrix_equal(Matrix2, Expected).
  990
  991    test(matrix_select1) :-
  992    matrix_new(doubles, [3, 2], [13, 2, 3, 10, 1, 6], Matrix1),
  993    matrix_new(doubles, [1, 2], [13, 2], Expected),
  994    matrix_select(Matrix1, 1, [0], Matrix2),
  995    matrix_equal(Matrix2, Expected).
  996
  997    test(matrix_select2) :-
  998    matrix_new(doubles, [3, 2], [13, 2, 3, 10, 1, 6], Matrix1),
  999    matrix_new(doubles, [2, 2], [13, 2, 3, 10], Expected),
 1000    matrix_select(Matrix1, 1, [0, 1], Matrix2),
 1001    matrix_equal(Matrix2, Expected).
 1002
 1003    test(matrix_select3) :-
 1004    matrix_new(doubles, [2, 3], [13, 2, 3, 10, 1, 6], Matrix1),
 1005    matrix_new(doubles, [2, 2], [13, 3, 10, 6], Expected),
 1006    matrix_select(Matrix1, 2, [0, 2], Matrix2),
 1007    matrix_equal(Matrix2, Expected).
 1008
 1009    test(matrixCreateZeroes1) :-
 1010    matrix_createZeroes(doubles, [3, 3], Matrix1),
 1011    matrix_new(doubles, [3, 3], [0, 0, 0, 0, 0, 0, 0, 0, 0], Expected),
 1012    matrix_equal(Matrix1, Expected).
 1013
 1014    test(matrixSetZeroes1) :-
 1015    matrix_new(doubles, [3, 3], Matrix1),
 1016    Matrix1=matrix(_, [3, 3], Elements),
 1017    matrix_new(doubles, [3, 3], [0, 0, 0, 0, 0, 0, 0, 0, 0], Expected),
 1018    matrixSetZeroes(3, 3, Elements),
 1019    matrix_equal(Matrix1, Expected).
 1020
 1021    test(matrixSetOnes1) :-
 1022    matrix_new(doubles, [3, 3], Matrix1),
 1023    Matrix1=matrix(_, [3, 3], Elements),
 1024    matrix_new(doubles, [3, 3], [1, 1, 1, 1, 1, 1, 1, 1, 1], Expected),
 1025    matrixSetOnes(3, 3, Elements),
 1026    matrix_equal(Matrix1, Expected).
 1027
 1028    test(matrix_setAllNumbers1) :-
 1029    matrix_new(doubles, [3, 3], Matrix1),
 1030    Matrix1=matrix(_, [3, 3], Elements),
 1031    Number is 5.0,
 1032    matrixSetAll(3, 3, Elements, Number),
 1033    matrix_new(doubles, [3, 3], [5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0], Expected),
 1034    matrix_equal(Matrix1, Expected).
 1035
 1036    test(matrix_eye1) :-
 1037    matrix_new(doubles, [3, 3], Matrix1),
 1038    Matrix1=matrix(_, [3, 3], Elements),
 1039    matrix_eye(3, 3, Elements),
 1040    matrix_new(doubles, [3, 3], [1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0], Expected),
 1041    matrix_equal(Matrix1, Expected).
 1042
 1043    test(matrix_equal1) :-
 1044    matrix_new(doubles, [3, 3], [1, 2, 3, 4, 5, 6, 7, 8, 9], Matrix1),
 1045    matrix_new(doubles, [3, 2], [1, 3, 4, 6, 7, 9], Expected),
 1046    matrix_select(Matrix1, 2, [0, 2], Matrix2),
 1047    matrix_equal(Matrix2, Expected).
 1048
 1049    test(matrix_map1) :-
 1050    matrix_new(doubles, [2, 3], [1, 2, 3, 11, 5, 6], Matrix1),
 1051    matrix_map(\E1^E2^(E2 is E1*2),
 1052               Matrix1,
 1053               Matrix2),
 1054    matrix_new(doubles, [2, 3], [2, 4, 6, 22, 10, 12], Expected),
 1055    matrix_equal(Matrix2, Expected).
 1056
 1057    test(matrix_foreach1) :-
 1058    matrix_new(doubles, [3, 3], [1, 2, 3, 4, 5, 6, 7, 8, 9], Matrix1),
 1059    matrix_new(doubles, [3, 3], [3, 3, 3, 3, 3, 3, 3, 3, 3], Expected),
 1060    matrix_foreach(Matrix1,
 1061                   \X^Y^myset(X, Y, 3)),
 1062    matrix_equal(Matrix1, Expected).
 1063
 1064    test(matrix_foldl1) :-
 1065    matrix_new(doubles, [3, 3], [1, 2, 3, 4, 5, 6, 7, 8, 9], Matrix1),
 1066    matrix_foldl(Matrix1, findMax, 9.0).
 1067
 1068    test(matrix_from_list1) :-
 1069    matrix_new(doubles, [3, 3], [1, 2, 3, 4, 5, 6, 7, 8, 9], Matrix1),
 1070    matrix_from_list([11, 12, 13, 14, 15, 16, 17, 18, 19], Matrix1),
 1071    matrix_new(doubles,
 1072               [3, 3],
 1073               [11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 17.0, 18.0, 19.0],
 1074               Expected),
 1075    matrix_equal(Matrix1, Expected).
 1076
 1077    test(cholesky1) :-
 1078    matrix_new(doubles, [3, 3], [25, 15, -5, 15, 18, 0, -5, 0, 11], Matrix1),
 1079    matrix_new(doubles, [3, 3], Matrix2),
 1080    matrix_new(doubles, [3, 3], [5.0, 0.0, 0.0, 3.0, 3.0, 0.0, -1.0, 1.0, 3.0], Expected),
 1081    matrix_sollower(Matrix1, Matrix2),
 1082    cholesky(Matrix2),
 1083    matrix_equal(Matrix2, Expected).
 1084
 1085end_tests(matrix)