Skip to content
This repository was archived by the owner on Feb 18, 2020. It is now read-only.

Commit 2de42bd

Browse files
henrique-borbahenrique-borba
committed
CArray::svd
1 parent 34ef83b commit 2de42bd

File tree

4 files changed

+122
-4
lines changed

4 files changed

+122
-4
lines changed

‎kernel/iterators.c

Lines changed: 0 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,3 @@
1-
//
2-
// Created by Henrique Borba on 19/11/2018.
3-
//
41
#include "php.h"
52
#include "iterators.h"
63
#include "carray.h"

‎kernel/linalg.c

Lines changed: 69 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -608,4 +608,73 @@ CArray_Vdot(CArray * a, CArray * b, MemoryPointer * out)
608608
//Py_XDECREF(ap2);
609609
//Py_XDECREF(ret);
610610
return NULL;
611+
}
612+
613+
614+
CArray **
615+
CArray_Svd(CArray * a, int full_matrices, int compute_uv, MemoryPointer * out)
616+
{
617+
int m, n;
618+
int lda, ldu, ldvt, info, lwork;
619+
int * iwork;
620+
double * s, * u, * vt;
621+
CArray * u_ca, * s_ca, *vh_ca, ** rtn;
622+
623+
if (CArray_NDIM(a) != 2) {
624+
throw_valueerror_exception("Expected 2D array");
625+
return NULL;
626+
}
627+
628+
m = CArray_DIMS(a)[0];
629+
n = CArray_DIMS(a)[1];
630+
631+
lda = m;
632+
ldu = m;
633+
ldvt = n;
634+
lwork = -1;
635+
636+
s = emalloc(sizeof(double) * n);
637+
u = emalloc(sizeof(double) * (ldu * m));
638+
vt = emalloc(sizeof(double) * (ldvt * n));
639+
640+
info = LAPACKE_dgesdd(LAPACK_ROW_MAJOR, 'S', m, n, DDATA(a), lda, s, u, ldu, vt, ldvt);
641+
642+
if( info > 0 ) {
643+
throw_valueerror_exception( "The algorithm computing SVD failed to converge." );
644+
return NULL;
645+
}
646+
647+
u_ca = emalloc(sizeof(CArray));
648+
u_ca = CArray_NewFromDescr_int(u_ca, CArray_DESCR(a), CArray_NDIM(a), CArray_DIMS(a), CArray_STRIDES(a),
649+
NULL, 0, NULL, 1, 0);
650+
efree(u_ca->data);
651+
u_ca->data = (char *)u;
652+
CArrayDescriptor_INCREF(CArray_DESCR(a));
653+
654+
s_ca = emalloc(sizeof(CArray));
655+
s_ca = CArray_NewFromDescr_int(s_ca, CArray_DESCR(a), 1, &n, NULL,
656+
NULL, 0, NULL, 1, 0);
657+
efree(s_ca->data);
658+
s_ca->data = (char *)s;
659+
CArrayDescriptor_INCREF(CArray_DESCR(a));
660+
661+
vh_ca = emalloc(sizeof(CArray));
662+
vh_ca = CArray_NewFromDescr_int(vh_ca, CArray_DESCR(a), CArray_NDIM(a), CArray_DIMS(a), CArray_STRIDES(a),
663+
NULL, 0, NULL, 1, 0);
664+
efree(vh_ca->data);
665+
vh_ca->data = (char *)vt;
666+
CArrayDescriptor_INCREF(CArray_DESCR(a));
667+
668+
if (out != NULL) {
669+
add_to_buffer(&(out[0]), u_ca, sizeof(CArray));
670+
add_to_buffer(&(out[1]), s_ca, sizeof(CArray));
671+
add_to_buffer(&(out[2]), vh_ca, sizeof(CArray));
672+
}
673+
674+
rtn = emalloc(sizeof(CArray *) * 3);
675+
rtn[0] = u_ca;
676+
rtn[1] = s_ca;
677+
rtn[2] = vh_ca;
678+
679+
return rtn;
611680
}

‎kernel/linalg.h

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -12,4 +12,5 @@ CArray * CArray_Inv(CArray * a, MemoryPointer * out);
1212
CArray * CArray_Norm(CArray * a, int norm, MemoryPointer * out);
1313
CArray * CArray_Det(CArray * a, MemoryPointer * out);
1414
CArray * CArray_Vdot(CArray * a, CArray * b, MemoryPointer * out);
15+
CArray ** CArray_Svd(CArray * a, int full_matrices, int compute_uv, MemoryPointer * out);
1516
#endif

‎phpsci.c

Lines changed: 52 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1088,7 +1088,55 @@ PHP_METHOD(CArray, vdot)
10881088

10891089
RETURN_MEMORYPOINTER(return_value, &rtn_ptr);
10901090
}
1091+
PHP_METHOD(CArray, svd)
1092+
{
1093+
int full_matrices_int = 1, compute_uv_int = 1, i = 0;
1094+
MemoryPointer a_ptr, * out_ptr;
1095+
zval * a;
1096+
zval * tmp_zval;
1097+
zend_bool full_matrices, compute_uv;
1098+
CArray ** rtn, * a_ca;
1099+
ZEND_PARSE_PARAMETERS_START(1, 1)
1100+
Z_PARAM_ZVAL(a)
1101+
Z_PARAM_OPTIONAL
1102+
Z_PARAM_BOOL(full_matrices)
1103+
Z_PARAM_BOOL(compute_uv)
1104+
ZEND_PARSE_PARAMETERS_END();
1105+
ZVAL_TO_MEMORYPOINTER(a, &a_ptr, NULL);
1106+
a_ca = CArray_FromMemoryPointer(&a_ptr);
1107+
if(ZEND_NUM_ARGS() == 2) {
1108+
if (full_matrices == IS_FALSE) {
1109+
full_matrices_int = 0;
1110+
}
1111+
}
1112+
if(ZEND_NUM_ARGS() == 3) {
1113+
if (full_matrices == IS_FALSE) {
1114+
full_matrices_int = 0;
1115+
}
1116+
if (compute_uv == IS_FALSE) {
1117+
compute_uv_int = 0;
1118+
}
1119+
}
1120+
out_ptr = emalloc(3 * sizeof(MemoryPointer));
1121+
rtn = CArray_Svd(a_ca, full_matrices_int, compute_uv_int, out_ptr);
10911122

1123+
FREE_FROM_MEMORYPOINTER(&a_ptr);
1124+
1125+
if (rtn == NULL) {
1126+
efree(out_ptr);
1127+
return;
1128+
}
1129+
1130+
array_init_size(return_value, 3);
1131+
for (i = 0; i < 3;i ++) {
1132+
tmp_zval = MEMORYPOINTER_TO_ZVAL(&(out_ptr[i]));
1133+
zend_hash_next_index_insert_new(Z_ARRVAL_P(return_value), tmp_zval);
1134+
efree(tmp_zval);
1135+
}
1136+
1137+
efree(out_ptr);
1138+
efree(rtn);
1139+
}
10921140

10931141

10941142

@@ -1138,6 +1186,7 @@ PHP_METHOD(CArray, ones)
11381186
dtype = emalloc(sizeof(char));
11391187
*dtype = 'd';
11401188
}
1189+
11411190
shape = ZVAL_TO_TUPLE(zshape, &ndim);
11421191
CArray_Ones(shape, ndim, dtype, &order, &ptr);
11431192
efree(shape);
@@ -1165,7 +1214,6 @@ PHP_METHOD(CArray, add)
11651214
target_ca2 = CArray_FromMemoryPointer(&target2_ptr);
11661215
output_ca = CArray_Add(target_ca1, target_ca2, &result_ptr);
11671216

1168-
11691217
FREE_FROM_MEMORYPOINTER(&target1_ptr);
11701218
FREE_FROM_MEMORYPOINTER(&target2_ptr);
11711219
if (output_ca != NULL) {
@@ -2512,6 +2560,9 @@ static zend_function_entry carray_class_methods[] =
25122560
PHP_ME(CArray, det, NULL, ZEND_ACC_PUBLIC | ZEND_ACC_STATIC)
25132561
PHP_ME(CArray, vdot, NULL, ZEND_ACC_PUBLIC | ZEND_ACC_STATIC)
25142562

2563+
// DECOMPOSITIONS
2564+
PHP_ME(CArray, svd, NULL, ZEND_ACC_PUBLIC | ZEND_ACC_STATIC)
2565+
25152566
// ARITHMETIC
25162567
PHP_ME(CArray, add, NULL, ZEND_ACC_PUBLIC | ZEND_ACC_STATIC)
25172568
PHP_ME(CArray, subtract, NULL, ZEND_ACC_PUBLIC | ZEND_ACC_STATIC)

0 commit comments

Comments
 (0)