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, 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).
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)]).
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).
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).
191matrix_dims(matrix(_, Dimension, _), Dimension).
195matrix_size(matrix(_, [NRows, NColumns], _), Nelements) :-
196 Nelements is NRows*NColumns.
200matrix_type(matrix(Type, _, _), Type).
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).
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
261
263myset(_, Number, Number).
264
265myset(_, 0).
267findMax(FirstNumber, SecondNumber, FirstNumber) :-
268 FirstNumber>=SecondNumber, !.
269
270findMax(_, SecondNumber, SecondNumber).
271
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).
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).
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).
445matrix_min(Matrix, Min) :-
446 matrix_foldl(Matrix,
447 \X^Y^Z^(X=<Y->Z=X;Z=Y),
448 Min).
452matrix_sum(Matrix, Sum) :-
453 matrix_foldl(Matrix,
454 \CurrentSum^Element^NewSum^(NewSum is CurrentSum+Element),
455 Sum).
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).
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]).
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]).
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).
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).
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).
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
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
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
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).
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).
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).
817matrix_read([NRows, NColumns], M1) :-
818 matrix_new(doubles, [NRows, NColumns], M1),
819 matrix_foreach(M1, \_^Y^read(Y)).
823matrixSetZeroes(NRows, NColumns, Elements) :-
824 matrixSetAll(NRows, NColumns, Elements, 0).
828matrixSetOnes(NRows, NColumns, Elements) :-
829 matrixSetAll(NRows, NColumns, Elements, 1).
833matrix_setAllNumbers(NRows, NColumns, Elements) :-
834 read(Number),
835 matrixSetAll(NRows, NColumns, Elements, Number).
839matrix_eye(NRows, NColumns, Elements) :-
840 matrixEye(NRows, NColumns, Elements).
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).
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
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)