1 #ifdef HAVE_CONFIG_H
2 #include "config.h"
3 #endif
4
5 #include <php.h>
6 #include <cblas.h>
7 #include <lapacke.h>
8 #include "kernel/operators.h"
9
10 /**
11 * Matrix-matrix multiplication i.e. linear transformation of matrices A and B.
12 *
13 * @param return_value
14 * @param a
15 * @param b
16 */
tensor_matmul(zval * return_value,zval * a,zval * b)17 void tensor_matmul(zval * return_value, zval * a, zval * b)
18 {
19 unsigned int i, j;
20 Bucket * row;
21 zval rowC, c;
22
23 zend_array * aa = Z_ARR_P(a);
24 zend_array * ab = Z_ARR_P(b);
25
26 Bucket * ba = aa->arData;
27 Bucket * bb = ab->arData;
28
29 unsigned int m = zend_array_count(aa);
30 unsigned int p = zend_array_count(ab);
31 unsigned int n = zend_array_count(Z_ARR(bb[0].val));
32
33 double * va = emalloc(m * p * sizeof(double));
34 double * vb = emalloc(n * p * sizeof(double));
35 double * vc = emalloc(m * n * sizeof(double));
36
37 for (i = 0; i < m; ++i) {
38 row = Z_ARR(ba[i].val)->arData;
39
40 for (j = 0; j < p; ++j) {
41 va[i * p + j] = zephir_get_doubleval(&row[j].val);
42 }
43 }
44
45 for (i = 0; i < p; ++i) {
46 row = Z_ARR(bb[i].val)->arData;
47
48 for (j = 0; j < n; ++j) {
49 vb[i * n + j] = zephir_get_doubleval(&row[j].val);
50 }
51 }
52
53 cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m, n, p, 1.0, va, p, vb, n, 0.0, vc, n);
54
55 array_init_size(&c, m);
56
57 for (i = 0; i < m; ++i) {
58 array_init_size(&rowC, n);
59
60 for (j = 0; j < n; ++j) {
61 add_next_index_double(&rowC, vc[i * n + j]);
62 }
63
64 add_next_index_zval(&c, &rowC);
65 }
66
67 RETVAL_ARR(Z_ARR(c));
68
69 efree(va);
70 efree(vb);
71 efree(vc);
72 }
73
74 /**
75 * Dot product between vectors A and B.
76 *
77 * @param return_value
78 * @param a
79 * @param b
80 */
tensor_dot(zval * return_value,zval * a,zval * b)81 void tensor_dot(zval * return_value, zval * a, zval * b)
82 {
83 unsigned int i;
84
85 zend_array * aa = Z_ARR_P(a);
86 zend_array * ab = Z_ARR_P(b);
87
88 Bucket * ba = aa->arData;
89 Bucket * bb = ab->arData;
90
91 unsigned int n = zend_array_count(aa);
92
93 double sigma = 0.0;
94
95 for (i = 0; i < n; ++i) {
96 sigma += zephir_get_doubleval(&ba[i].val) * zephir_get_doubleval(&bb[i].val);
97 }
98
99 RETVAL_DOUBLE(sigma);
100 }
101
102 /**
103 * Return the multiplicative inverse of a square matrix A.
104 *
105 * @param return_value
106 * @param a
107 */
tensor_inverse(zval * return_value,zval * a)108 void tensor_inverse(zval * return_value, zval * a)
109 {
110 unsigned int i, j;
111 Bucket * row;
112 zval rowB, b;
113
114 zend_array * aa = Z_ARR_P(a);
115
116 Bucket * ba = aa->arData;
117
118 unsigned int n = zend_array_count(aa);
119
120 double * va = emalloc(n * n * sizeof(double));
121 int * pivots = emalloc(n * sizeof(int));
122
123 for (i = 0; i < n; ++i) {
124 row = Z_ARR(ba[i].val)->arData;
125
126 for (j = 0; j < n; ++j) {
127 va[i * n + j] = zephir_get_doubleval(&row[j].val);
128 }
129 }
130
131 lapack_int status;
132
133 status = LAPACKE_dgetrf(LAPACK_ROW_MAJOR, n, n, va, n, pivots);
134
135 if (status != 0) {
136 RETURN_NULL();
137 }
138
139 status = LAPACKE_dgetri(LAPACK_ROW_MAJOR, n, va, n, pivots);
140
141 if (status != 0) {
142 RETURN_NULL();
143 }
144
145 array_init_size(&b, n);
146
147 for (i = 0; i < n; ++i) {
148 array_init_size(&rowB, n);
149
150 for (j = 0; j < n; ++j) {
151 add_next_index_double(&rowB, va[i * n + j]);
152 }
153
154 add_next_index_zval(&b, &rowB);
155 }
156
157 RETVAL_ARR(Z_ARR(b));
158
159 efree(va);
160 efree(pivots);
161 }
162
163 /**
164 * Return the (Moore-Penrose) pseudoinverse of a general matrix A.
165 *
166 * @param return_value
167 * @param a
168 */
tensor_pseudoinverse(zval * return_value,zval * a)169 void tensor_pseudoinverse(zval * return_value, zval * a)
170 {
171 unsigned int i, j;
172 Bucket * row;
173 zval b, rowB;
174
175 zend_array * aa = Z_ARR_P(a);
176
177 Bucket * ba = aa->arData;
178
179 unsigned int m = zend_array_count(aa);
180 unsigned int n = zend_array_count(Z_ARR(ba[0].val));
181 unsigned int k = MIN(m, n);
182
183 double * va = emalloc(m * n * sizeof(double));
184 double * vu = emalloc(m * m * sizeof(double));
185 double * vs = emalloc(k * sizeof(double));
186 double * vvt = emalloc(n * n * sizeof(double));
187 double * vb = emalloc(n * m * sizeof(double));
188
189 for (i = 0; i < m; ++i) {
190 row = Z_ARR(ba[i].val)->arData;
191
192 for (j = 0; j < n; ++j) {
193 va[i * n + j] = zephir_get_doubleval(&row[j].val);
194 }
195 }
196
197 lapack_int status = LAPACKE_dgesdd(LAPACK_ROW_MAJOR, 'A', m, n, va, n, vs, vu, m, vvt, n);
198
199 if (status != 0) {
200 RETURN_NULL();
201 }
202
203 for (i = 0; i < k; ++i) {
204 cblas_dscal(m, 1.0 / vs[i], &vu[i], m);
205 }
206
207 cblas_dgemm(CblasRowMajor, CblasTrans, CblasTrans, n, m, m, 1.0, vvt, n, vu, m, 0.0, vb, m);
208
209 array_init_size(&b, n);
210
211 for (i = 0; i < n; ++i) {
212 array_init_size(&rowB, m);
213
214 for (j = 0; j < m; ++j) {
215 add_next_index_double(&rowB, vb[i * m + j]);
216 }
217
218 add_next_index_zval(&b, &rowB);
219 }
220
221 RETVAL_ARR(Z_ARR(b));
222
223 efree(va);
224 efree(vu);
225 efree(vs);
226 efree(vvt);
227 efree(vb);
228 }
229
230 /**
231 * Return the row echelon form of matrix A.
232 *
233 * @param return_value
234 * @param a
235 */
tensor_ref(zval * return_value,zval * a)236 void tensor_ref(zval * return_value, zval * a)
237 {
238 unsigned int i, j;
239 Bucket * row;
240 zval rowB, b;
241 zval tuple;
242
243 zend_array * aa = Z_ARR_P(a);
244
245 Bucket * ba = aa->arData;
246
247 unsigned int m = zend_array_count(aa);
248 unsigned int n = zend_array_count(Z_ARR(ba[0].val));
249
250 double * va = emalloc(m * n * sizeof(double));
251 int * pivots = emalloc(MIN(m, n) * sizeof(int));
252
253 for (i = 0; i < m; ++i) {
254 row = Z_ARR(ba[i].val)->arData;
255
256 for (j = 0; j < n; ++j) {
257 va[i * n + j] = zephir_get_doubleval(&row[j].val);
258 }
259 }
260
261 lapack_int status = LAPACKE_dgetrf(LAPACK_ROW_MAJOR, m, n, va, n, pivots);
262
263 if (status != 0) {
264 RETURN_NULL();
265 }
266
267 array_init_size(&b, m);
268
269 long swaps = 0;
270
271 for (i = 0; i < m; ++i) {
272 array_init_size(&rowB, n);
273
274 for (j = 0; j < i; ++j) {
275 add_next_index_double(&rowB, 0.0);
276 }
277
278 for (j = i; j < n; ++j) {
279 add_next_index_double(&rowB, va[i * n + j]);
280 }
281
282 add_next_index_zval(&b, &rowB);
283
284 if (i + 1 != pivots[i]) {
285 ++swaps;
286 }
287 }
288
289 array_init_size(&tuple, 2);
290
291 add_next_index_zval(&tuple, &b);
292 add_next_index_long(&tuple, swaps);
293
294 RETVAL_ARR(Z_ARR(tuple));
295
296 efree(va);
297 efree(pivots);
298 }
299
300 /**
301 * Compute the Cholesky decomposition of matrix A and return the lower triangular matrix.
302 *
303 * @param return_value
304 * @param a
305 */
tensor_cholesky(zval * return_value,zval * a)306 void tensor_cholesky(zval * return_value, zval * a)
307 {
308 unsigned int i, j;
309 Bucket * row;
310 zval rowB, b;
311
312 zend_array * aa = Z_ARR_P(a);
313
314 Bucket * ba = aa->arData;
315
316 unsigned int n = zend_array_count(aa);
317
318 double * va = emalloc(n * n * sizeof(double));
319
320 for (i = 0; i < n; ++i) {
321 row = Z_ARR(ba[i].val)->arData;
322
323 for (j = 0; j < n; ++j) {
324 va[i * n + j] = zephir_get_doubleval(&row[j].val);
325 }
326 }
327
328 lapack_int status = LAPACKE_dpotrf(LAPACK_ROW_MAJOR, 'L', n, va, n);
329
330 if (status != 0) {
331 RETURN_NULL();
332 }
333
334 array_init_size(&b, n);
335
336 for (i = 0; i < n; ++i) {
337 array_init_size(&rowB, n);
338
339 for (j = 0; j <= i; ++j) {
340 add_next_index_double(&rowB, va[i * n + j]);
341 }
342
343 for (j = i + 1; j < n; ++j) {
344 add_next_index_double(&rowB, 0.0);
345 }
346
347 add_next_index_zval(&b, &rowB);
348 }
349
350 RETVAL_ARR(Z_ARR(b));
351
352 efree(va);
353 }
354
355 /**
356 * Compute the LU factorization of matrix A and return a tuple with lower, upper, and permutation matrices.
357 *
358 * @param return_value
359 * @param a
360 */
tensor_lu(zval * return_value,zval * a)361 void tensor_lu(zval * return_value, zval * a)
362 {
363 unsigned int i, j;
364 Bucket * row;
365 zval rowL, l, rowU, u, rowP, p;
366 zval tuple;
367
368 zend_array * aa = Z_ARR_P(a);
369
370 Bucket * ba = aa->arData;
371
372 unsigned int n = zend_array_count(aa);
373
374 double * va = emalloc(n * n * sizeof(double));
375 int * pivots = emalloc(n * sizeof(int));
376
377 for (i = 0; i < n; ++i) {
378 row = Z_ARR(ba[i].val)->arData;
379
380 for (j = 0; j < n; ++j) {
381 va[i * n + j] = zephir_get_doubleval(&row[j].val);
382 }
383 }
384
385 lapack_int status = LAPACKE_dgetrf(LAPACK_ROW_MAJOR, n, n, va, n, pivots);
386
387 if (status != 0) {
388 RETURN_NULL();
389 }
390
391 array_init_size(&l, n);
392 array_init_size(&u, n);
393 array_init_size(&p, n);
394
395 for (i = 0; i < n; ++i) {
396 array_init_size(&rowL, n);
397
398 for (j = 0; j < i; ++j) {
399 add_next_index_double(&rowL, va[i * n + j]);
400 }
401
402 add_next_index_double(&rowL, 1.0);
403
404 for (j = i + 1; j < n; ++j) {
405 add_next_index_double(&rowL, 0.0);
406 }
407
408 add_next_index_zval(&l, &rowL);
409 }
410
411 for (i = 0; i < n; ++i) {
412 array_init_size(&rowU, n);
413
414 for (j = 0; j < i; ++j) {
415 add_next_index_double(&rowU, 0.0);
416 }
417
418 for (j = i; j < n; ++j) {
419 add_next_index_double(&rowU, va[i * n + j]);
420 }
421
422 add_next_index_zval(&u, &rowU);
423 }
424
425 for (i = 0; i < n; ++i) {
426 array_init_size(&rowP, n);
427
428 for (j = 0; j < n; ++j) {
429 if (j == pivots[i] - 1) {
430 add_next_index_long(&rowP, 1);
431 } else {
432 add_next_index_long(&rowP, 0);
433 }
434 }
435
436 add_next_index_zval(&p, &rowP);
437 }
438
439 array_init_size(&tuple, 3);
440
441 add_next_index_zval(&tuple, &l);
442 add_next_index_zval(&tuple, &u);
443 add_next_index_zval(&tuple, &p);
444
445 RETVAL_ARR(Z_ARR(tuple));
446
447 efree(va);
448 efree(pivots);
449 }
450
451 /**
452 * Compute the eigendecomposition of a general matrix A and return the eigenvalues and eigenvectors in a tuple.
453 *
454 * @param return_value
455 * @param a
456 */
tensor_eig(zval * return_value,zval * a)457 void tensor_eig(zval * return_value, zval * a)
458 {
459 unsigned int i, j;
460 Bucket * row;
461 zval eigenvalues;
462 zval eigenvectors;
463 zval eigenvector;
464 zval tuple;
465
466 zend_array * aa = Z_ARR_P(a);
467
468 Bucket * ba = aa->arData;
469
470 unsigned int n = zend_array_count(aa);
471
472 double * va = emalloc(n * n * sizeof(double));
473 double * wr = emalloc(n * sizeof(double));
474 double * wi = emalloc(n * sizeof(double));
475 double * vr = emalloc(n * n * sizeof(double));
476
477 for (i = 0; i < n; ++i) {
478 row = Z_ARR(ba[i].val)->arData;
479
480 for (j = 0; j < n; ++j) {
481 va[i * n + j] = zephir_get_doubleval(&row[j].val);
482 }
483 }
484
485 lapack_int status = LAPACKE_dgeev(LAPACK_ROW_MAJOR, 'N', 'V', n, va, n, wr, wi, NULL, n, vr, n);
486
487 if (status != 0) {
488 RETURN_NULL();
489 }
490
491 array_init_size(&eigenvalues, n);
492 array_init_size(&eigenvectors, n);
493
494 for (i = 0; i < n; ++i) {
495 add_next_index_double(&eigenvalues, wr[i]);
496
497 array_init_size(&eigenvector, n);
498
499 for (j = 0; j < n; ++j) {
500 add_next_index_double(&eigenvector, vr[i * n + j]);
501 }
502
503 add_next_index_zval(&eigenvectors, &eigenvector);
504 }
505
506 array_init_size(&tuple, 2);
507
508 add_next_index_zval(&tuple, &eigenvalues);
509 add_next_index_zval(&tuple, &eigenvectors);
510
511 RETVAL_ARR(Z_ARR(tuple));
512
513 efree(va);
514 efree(wr);
515 efree(wi);
516 efree(vr);
517 }
518
519 /**
520 * Compute the eigendecomposition of a symmetric matrix A and return the eigenvalues and eigenvectors in a tuple.
521 *
522 * @param return_value
523 * @param a
524 */
tensor_eig_symmetric(zval * return_value,zval * a)525 void tensor_eig_symmetric(zval * return_value, zval * a)
526 {
527 unsigned int i, j;
528 Bucket * row;
529 zval eigenvalues;
530 zval eigenvectors;
531 zval eigenvector;
532 zval tuple;
533
534 zend_array * aa = Z_ARR_P(a);
535
536 Bucket * ba = aa->arData;
537
538 unsigned int n = zend_array_count(aa);
539
540 double * va = emalloc(n * n * sizeof(double));
541 double * wr = emalloc(n * sizeof(double));
542
543 for (i = 0; i < n; ++i) {
544 row = Z_ARR(ba[i].val)->arData;
545
546 for (j = 0; j < n; ++j) {
547 va[i * n + j] = zephir_get_doubleval(&row[j].val);
548 }
549 }
550
551 lapack_int status = LAPACKE_dsyev(LAPACK_ROW_MAJOR, 'V', 'U', n, va, n, wr);
552
553 if (status != 0) {
554 RETURN_NULL();
555 }
556
557 array_init_size(&eigenvalues, n);
558 array_init_size(&eigenvectors, n);
559
560 for (i = 0; i < n; ++i) {
561 add_next_index_double(&eigenvalues, wr[i]);
562
563 array_init_size(&eigenvector, n);
564
565 for (j = 0; j < n; ++j) {
566 add_next_index_double(&eigenvector, va[i * n + j]);
567 }
568
569 add_next_index_zval(&eigenvectors, &eigenvector);
570 }
571
572 array_init_size(&tuple, 2);
573
574 add_next_index_zval(&tuple, &eigenvalues);
575 add_next_index_zval(&tuple, &eigenvectors);
576
577 RETVAL_ARR(Z_ARR(tuple));
578
579 efree(va);
580 efree(wr);
581 }
582
583 /**
584 * Compute the singular value decomposition of a matrix A and return the singular values and unitary matrices U and VT in a tuple.
585 *
586 * @param return_value
587 * @param a
588 */
tensor_svd(zval * return_value,zval * a)589 void tensor_svd(zval * return_value, zval * a)
590 {
591 unsigned int i, j;
592 Bucket * row;
593 zval u, rowU;
594 zval s;
595 zval vt, rowVt;
596 zval tuple;
597
598 zend_array * aa = Z_ARR_P(a);
599
600 Bucket * ba = aa->arData;
601
602 unsigned int m = zend_array_count(aa);
603 unsigned int n = zend_array_count(Z_ARR(ba[0].val));
604 unsigned int k = MIN(m, n);
605
606 double * va = emalloc(m * n * sizeof(double));
607 double * vu = emalloc(m * m * sizeof(double));
608 double * vs = emalloc(k * sizeof(double));
609 double * vvt = emalloc(n * n * sizeof(double));
610
611 for (i = 0; i < m; ++i) {
612 row = Z_ARR(ba[i].val)->arData;
613
614 for (j = 0; j < n; ++j) {
615 va[i * n + j] = zephir_get_doubleval(&row[j].val);
616 }
617 }
618
619 lapack_int status = LAPACKE_dgesdd(LAPACK_ROW_MAJOR, 'A', m, n, va, n, vs, vu, m, vvt, n);
620
621 if (status != 0) {
622 RETURN_NULL();
623 }
624
625 array_init_size(&u, m);
626 array_init_size(&s, k);
627 array_init_size(&vt, n);
628
629 for (i = 0; i < m; ++i) {
630 array_init_size(&rowU, m);
631
632 for (j = 0; j < m; ++j) {
633 add_next_index_double(&rowU, vu[i * m + j]);
634 }
635
636 add_next_index_zval(&u, &rowU);
637 }
638
639 for (i = 0; i < k; ++i) {
640 add_next_index_double(&s, vs[i]);
641 }
642
643 for (i = 0; i < n; ++i) {
644 array_init_size(&rowVt, n);
645
646 for (j = 0; j < n; ++j) {
647 add_next_index_double(&rowVt, vvt[i * n + j]);
648 }
649
650 add_next_index_zval(&vt, &rowVt);
651 }
652
653 array_init_size(&tuple, 3);
654
655 add_next_index_zval(&tuple, &u);
656 add_next_index_zval(&tuple, &s);
657 add_next_index_zval(&tuple, &vt);
658
659 RETVAL_ARR(Z_ARR(tuple));
660
661 efree(va);
662 efree(vu);
663 efree(vs);
664 efree(vvt);
665 }
666