55 if (rows == 0 || cols == 0)
58 matd_t *m = calloc(1,
sizeof(
matd_t) + (rows*cols*
sizeof(
double)));
77 if (rows == 0 || cols == 0)
81 for (
int i = 0; i < rows * cols; i++)
89 if (rows == 0 || cols == 0)
93 for (
int i = 0; i < rows * cols; i++)
94 m->
data[i] = (
double)data[i];
105 for (
int i = 0; i < dim; i++)
117 assert(row < m->nrows);
119 assert(col < m->ncols);
135 assert(row < m->nrows);
137 assert(col < m->ncols);
175 assert(r0 >= 0 && r0 < a->nrows);
176 assert(c0 >= 0 && c0 < a->ncols);
178 int nrows = r1 - r0 + 1;
179 int ncols = c1 - c0 + 1;
183 for (
int row = r0; row <= r1; row++)
184 for (
int col = c0; col <= c1; col++)
199 for (
int i = 0; i < m->
nrows; i++) {
200 for (
int j = 0; j < m->
ncols; j++) {
217 for (
int j = 0; j < m->
ncols; j++) {
218 for (
int i = 0; i < m->
nrows; i++) {
248 for (
int i = 0; i < m->
nrows; i++) {
249 for (
int j = 0; j < m->
ncols; j++) {
251 for (
int k = 0; k < a->
ncols; k++) {
270 for (
int i = 0; i < m->
nrows; i++) {
271 for (
int j = 0; j < m->
ncols; j++) {
288 for (
int i = 0; i < a->
nrows; i++) {
289 for (
int j = 0; j < a->
ncols; j++) {
307 for (
int i = 0; i < m->
nrows; i++) {
308 for (
int j = 0; j < m->
ncols; j++) {
328 for (
int i = 0; i < a->
nrows; i++) {
329 for (
int j = 0; j < a->
ncols; j++) {
348 for (
int i = 0; i < m->
nrows; i++) {
349 for (
int j = 0; j < m->
ncols; j++) {
369 for (
int i = 0; i < a->
nrows; i++) {
370 for (
int j = 0; j < a->
ncols; j++) {
386 for (
int i = 0; i < a->
nrows; i++) {
387 for (
int j = 0; j < a->
ncols; j++) {
404 double detL = 1;
double detU = 1;
405 for (
int i = 0; i < a->
nrows; i++) {
415 double det = mlu->
pivsign * detL * detU;
433 assert(a->
nrows > 0);
455 double m00 =
MATD_EL(a,0,0), m01 =
MATD_EL(a,0,1), m02 =
MATD_EL(a,0,2), m03 =
MATD_EL(a,0,3);
456 double m10 =
MATD_EL(a,1,0), m11 =
MATD_EL(a,1,1), m12 =
MATD_EL(a,1,2), m13 =
MATD_EL(a,1,3);
457 double m20 =
MATD_EL(a,2,0), m21 =
MATD_EL(a,2,1), m22 =
MATD_EL(a,2,2), m23 =
MATD_EL(a,2,3);
458 double m30 =
MATD_EL(a,3,0), m31 =
MATD_EL(a,3,1), m32 =
MATD_EL(a,3,2), m33 =
MATD_EL(a,3,3);
460 return m00 * m11 * m22 * m33 - m00 * m11 * m23 * m32 -
461 m00 * m21 * m12 * m33 + m00 * m21 * m13 * m32 + m00 * m31 * m12 * m23 -
462 m00 * m31 * m13 * m22 - m10 * m01 * m22 * m33 +
463 m10 * m01 * m23 * m32 + m10 * m21 * m02 * m33 -
464 m10 * m21 * m03 * m32 - m10 * m31 * m02 * m23 +
465 m10 * m31 * m03 * m22 + m20 * m01 * m12 * m33 -
466 m20 * m01 * m13 * m32 - m20 * m11 * m02 * m33 +
467 m20 * m11 * m03 * m32 + m20 * m31 * m02 * m13 -
468 m20 * m31 * m03 * m12 - m30 * m01 * m12 * m23 +
469 m30 * m01 * m13 * m22 + m30 * m11 * m02 * m23 -
470 m30 * m11 * m03 * m22 - m30 * m21 * m02 * m13 +
471 m30 * m21 * m03 * m12;
501 double det = x->
data[0];
505 double invdet = 1.0 / det;
508 MATD_EL(m, 0, 0) = 1.0 * invdet;
517 double invdet = 1.0 / det;
558 while (expr[*pos] != 0) {
560 switch (expr[*pos]) {
565 garb[*garbpos] = res;
576 assert(expr[*pos+1] ==
'-');
577 assert(expr[*pos+2] ==
'1');
580 garb[*garbpos] = res;
599 matd_t **garb,
int *garbpos,
int oneterm)
601 while (expr[*pos] != 0) {
603 switch (expr[*pos]) {
606 if (oneterm && acc != NULL)
616 garb[*garbpos] = res;
642 garb[*garbpos] = res;
651 matd_t *rhs = args[*argpos];
652 garb[*garbpos] = rhs;
664 garb[*garbpos] = res;
673 matd_t *rhs = args[*argpos];
684 garb[*garbpos] = res;
714 const char *start = &expr[*pos];
716 double s = strtod(start, &end);
717 (*pos) += (end - start);
719 garb[*garbpos] = rhs;
728 garb[*garbpos] = res;
737 if (oneterm && acc != NULL)
748 garb[*garbpos] = res;
755 if (oneterm && acc != NULL)
765 garb[*garbpos] = res;
775 garb[*garbpos] = res;
789 fprintf(stderr,
"matd_op(): Unknown character: '%c'\n", expr[*pos]);
790 assert(expr[*pos] != expr[*pos]);
803 assert(expr != NULL);
805 for (
const char *p = expr; *p != 0; p++) {
806 if (*p ==
'M' || *p ==
'F')
820 for (
int i = 0; i < nargs; i++) {
821 args[i] = va_arg(ap,
matd_t*);
840 for (
int i = 0; i < garbpos; i++) {
854 for (
int i = 0; i < len; i++)
879 assert(n <= lena && n <= lenb);
882 for (
int i = 0; i < n; i++)
893 for (
int i = 0; i < maxcol; i++) {
896 double v = fabs(
MATD_EL(A, row, i));
913 assert(adim == bdim);
916 for (
int i = 0; i < adim; i++) {
934 for(
int i = 0; i < len; i++)
962 for (
int i = 0; i < a->
nrows; i++) {
963 for (
int j = 0; j < a->
ncols; j++) {
967 TYPE err = fabs(av - bv);
968 maxf = fmax(maxf, err);
1007 for (
int hhidx = 0; hhidx < A->
nrows; hhidx++) {
1009 if (hhidx < A->ncols) {
1027 int vlen = A->
nrows - hhidx;
1032 for (
int i = 0; i < vlen; i++) {
1033 v[i] =
MATD_EL(B, hhidx+i, hhidx);
1037 double oldv0 = v[0];
1043 mag2 += -oldv0*oldv0 + v[0]*v[0];
1046 double mag = sqrt(mag2);
1052 for (
int i = 0; i < vlen; i++)
1065 for (
int i = 0; i < LS->
nrows; i++) {
1067 for (
int j = 0; j < vlen; j++)
1068 dot +=
MATD_EL(LS, i, hhidx+j) * v[j];
1069 for (
int j = 0; j < vlen; j++)
1070 MATD_EL(LS, i, hhidx+j) -= 2*dot*v[j];
1074 for (
int i = 0; i < B->
ncols; i++) {
1076 for (
int j = 0; j < vlen; j++)
1077 dot +=
MATD_EL(B, hhidx+j, i) * v[j];
1078 for (
int j = 0; j < vlen; j++)
1079 MATD_EL(B, hhidx+j, i) -= 2*dot*v[j];
1083 if (hhidx+2 < A->
ncols) {
1084 int vlen = A->
ncols - hhidx - 1;
1089 for (
int i = 0; i < vlen; i++) {
1090 v[i] =
MATD_EL(B, hhidx, hhidx+i+1);
1094 double oldv0 = v[0];
1100 mag2 += -oldv0*oldv0 + v[0]*v[0];
1103 double mag = sqrt(mag2);
1109 for (
int i = 0; i < vlen; i++)
1119 for (
int i = 0; i < RS->
nrows; i++) {
1121 for (
int j = 0; j < vlen; j++)
1122 dot +=
MATD_EL(RS, i, hhidx+1+j) * v[j];
1123 for (
int j = 0; j < vlen; j++)
1124 MATD_EL(RS, i, hhidx+1+j) -= 2*dot*v[j];
1128 for (
int i = 0; i < B->
nrows; i++) {
1130 for (
int j = 0; j < vlen; j++)
1131 dot +=
MATD_EL(B, i, hhidx+1+j) * v[j];
1132 for (
int j = 0; j < vlen; j++)
1133 MATD_EL(B, i, hhidx+1+j) -= 2*dot*v[j];
1141 int maxiters = 1UL << 30;
1142 assert(maxiters > 0);
1151 const int find_max_method = 1;
1155 int maxrowidx[B->
ncols];
1156 int lastmaxi, lastmaxj;
1158 if (find_max_method == 1) {
1159 for (
int i = 2; i < B->
ncols; i++)
1165 lastmaxi = 0, lastmaxj = 1;
1168 for (iter = 0; iter < maxiters; iter++) {
1178 if (find_max_method == 1) {
1192 for (
int rowi = 0; rowi < B->
ncols; rowi++) {
1200 if (rowi == lastmaxi || rowi == lastmaxj) {
1202 thismaxv = fabs(
MATD_EL(B, rowi, maxrowidx[rowi]));
1210 if (maxrowidx[rowi] == lastmaxi || maxrowidx[rowi] == lastmaxj) {
1212 thismaxv = fabs(
MATD_EL(B, rowi, maxrowidx[rowi]));
1222 thismaxv = fabs(
MATD_EL(B, rowi, maxrowidx[rowi]));
1225 if (lastmaxi != rowi) {
1226 double v = fabs(
MATD_EL(B, rowi, lastmaxi));
1229 maxrowidx[rowi] = lastmaxi;
1234 if (lastmaxj != rowi) {
1235 double v = fabs(
MATD_EL(B, rowi, lastmaxj));
1238 maxrowidx[rowi] = lastmaxj;
1244 if (thismaxv > maxv) {
1251 maxj = maxrowidx[maxi];
1260 }
else if (find_max_method == 2) {
1265 for (
int i = 0; i < B->
ncols; i++) {
1266 for (
int j = 0; j < B->
ncols; j++) {
1270 double v = fabs(
MATD_EL(B, i, j));
1292 double A0 =
MATD_EL(B, maxi, maxi);
1293 double A1 =
MATD_EL(B, maxi, maxj);
1294 double A2 =
MATD_EL(B, maxj, maxi);
1295 double A3 =
MATD_EL(B, maxj, maxj);
1304 double U[4], S[2], V[4];
1331 for (
int i = 0; i < LS->
nrows; i++) {
1332 double vi =
MATD_EL(LS, i, maxi);
1333 double vj =
MATD_EL(LS, i, maxj);
1335 MATD_EL(LS, i, maxi) = U[0]*vi + U[2]*vj;
1336 MATD_EL(LS, i, maxj) = U[1]*vi + U[3]*vj;
1340 for (
int i = 0; i < RS->
nrows; i++) {
1341 double vi =
MATD_EL(RS, i, maxi);
1342 double vj =
MATD_EL(RS, i, maxj);
1344 MATD_EL(RS, i, maxi) = V[0]*vi + V[2]*vj;
1345 MATD_EL(RS, i, maxj) = V[1]*vi + V[3]*vj;
1350 for (
int i = 0; i < B->
ncols; i++) {
1351 double vi =
MATD_EL(B, maxi, i);
1352 double vj =
MATD_EL(B, maxj, i);
1354 MATD_EL(B, maxi, i) = U[0]*vi + U[2]*vj;
1355 MATD_EL(B, maxj, i) = U[1]*vi + U[3]*vj;
1359 for (
int i = 0; i < B->
nrows; i++) {
1360 double vi =
MATD_EL(B, i, maxi);
1361 double vj =
MATD_EL(B, i, maxj);
1363 MATD_EL(B, i, maxi) = V[0]*vi + V[2]*vj;
1364 MATD_EL(B, i, maxj) = V[1]*vi + V[3]*vj;
1370 printf(
"WARNING: maximum iters (maximum = %d, matrix %d x %d, max=%.15f)\n",
1379 double vals[A->
ncols];
1380 for (
int i = 0; i < A->
ncols; i++) {
1390 for (
int i = 0; i + 1 < A->
ncols; i++) {
1391 if (fabs(vals[i+1]) > fabs(vals[i])) {
1393 idxs[i] = idxs[i+1];
1396 double tmpv = vals[i];
1397 vals[i] = vals[i+1];
1408 for (
int i = 0; i < A->
ncols; i++) {
1409 MATD_EL(LP, idxs[i], idxs[i]) = 0;
1410 MATD_EL(RP, idxs[i], idxs[i]) = 0;
1412 MATD_EL(LP, idxs[i], i) = vals[i] < 0 ? -1 : 1;
1420 B =
matd_op(
"M'*F*M", LP, B, RP);
1430 memset(&res, 0,
sizeof(res));
1434 for (
int i = 0; i < B->
nrows; i++) {
1435 for (
int j = 0; j < B->
ncols; j++) {
1467 memset(&res, 0,
sizeof(res));
1503 unsigned int *piv = calloc(a->
nrows,
sizeof(
unsigned int));
1512 for (
int i = 0; i < a->
nrows; i++)
1515 for (
int j = 0; j < a->
ncols; j++) {
1516 for (
int i = 0; i < a->
nrows; i++) {
1517 int kmax = i < j ? i : j;
1521 for (
int k = 0; k < kmax; k++)
1530 for (
int i = j+1; i < lu->
nrows; i++) {
1549 double LUjj =
MATD_EL(lu, j, j);
1566 if (j < lu->ncols && j < lu->nrows && LUjj != 0) {
1568 for (
int i = j+1; i < lu->
nrows; i++)
1594 for (
int i = 0; i < lu->
ncols; i++)
1606 for (
int i = 0; i < lu->
nrows; i++) {
1618 for (
int i = 0; i < lu->
nrows; i++) {
1621 for (
int j = 0; j < i; j++) {
1634 for (
int i = 0; i < lu->
ncols; i++) {
1635 for (
int j = 0; j < lu->
ncols; j++) {
1653 for (
int i = 0; i < mlu->
lu->
nrows; i++)
1657 for (
int k = 0; k < mlu->
lu->
nrows; k++) {
1658 for (
int i = k+1; i < mlu->
lu->
nrows; i++) {
1660 for (
int t = 0; t < b->
ncols; t++)
1666 for (
int k = mlu->
lu->
ncols-1; k >= 0; k--) {
1667 double LUkk = 1.0 /
MATD_EL(mlu->
lu, k, k);
1668 for (
int t = 0; t < b->
ncols; t++)
1671 for (
int i = 0; i < k; i++) {
1673 for (
int t = 0; t < b->
ncols; t++)
1694 int v = random()&31;
1699 static double randf()
1701 double v = 1.0 *random() / RAND_MAX;
1705 int main(
int argc,
char *argv[])
1711 for (
int iter = 0; 1; iter++) {
1714 if (iter % 1000 == 0)
1715 printf(
"%d\n", iter);
1717 int m = 1 + (random()%(maxdim-1));
1718 int n = 1 + (random()%(maxdim-1));
1720 for (
int i = 0; i < m*n; i++)
1721 A->
data[i] = randi();
1759 for (
int iter = 0; 1; iter++) {
1760 if (iter % 100000 == 0)
1761 printf(
"%d\n", iter);
1768 matd_svd22_impl(A->
data, &s);
1770 memcpy(U->
data, s.U, 4*
sizeof(
double));
1771 memcpy(V->
data, s.V, 4*
sizeof(
double));
1775 assert(s.S[0] >= s.S[1]);
1776 assert(s.S[0] >= 0);
1777 assert(s.S[1] >= 0);
1789 for (
int i = 0; i < 4; i++)
1790 maxerr = fmax(maxerr, fabs(USV->
data[i] - A->
data[i]));
1793 printf(
"------------------------------------\n");
1796 printf(
"\nUSV':\n");
1798 printf(
"maxerr: %.15f\n", maxerr);
1804 assert(maxerr < 0.00001);
1858 for (
int i = 0; i < N; i++) {
1866 for (
int j = i; j < N; j++)
1869 for (
int j = i+1; j < N; j++) {
1875 for (
int k = j; k < N; k++) {
1897 memcpy(x, b, n*
sizeof(
TYPE));
1898 for (
int i = 0; i < n; i++) {
1901 for (
int j = i+1; j < u->
ncols; j++) {
1902 x[j] -= x[i] *
MATD_EL(u, i, j);
1912 for (
int i = 0; i < n; i++) {
1915 for (
int j = 0; j < i; j++) {
1919 x[i] = acc /
MATD_EL(L, i, i);
1926 for (
int i = u->
ncols-1; i >= 0; i--) {
1929 double diag =
MATD_EL(u, i, i);
1931 for (
int j = i+1; j < u->
ncols; j++)
1948 for (
int i = 0; i < u->
nrows; i++) {
1949 for (
int j = 0; j < i; j++) {
1953 for (
int k = 0; k < b->
ncols; k++) {
1958 for (
int k = 0; k < b->
ncols; k++) {
1964 for (
int k = u->
ncols-1; k >= 0; k--) {
1965 double LUkk = 1.0 /
MATD_EL(u, k, k);
1966 for (
int t = 0; t < b->
ncols; t++)
1969 for (
int i = 0; i < k; i++) {
1970 double LUik = -
MATD_EL(u, i, k);
1971 for (
int t = 0; t < b->
ncols; t++)
2006 double d = -DBL_MAX;
2007 for(
int x=0; x<m->
nrows; x++) {
2008 for(
int y=0; y<m->
ncols; y++) {
matd_svd_t matd_svd_flags(matd_t *A, int flags)
matd_t * matd_plu_l(const matd_plu_t *mlu)
matd_plu_t * matd_plu(const matd_t *a)
double matd_det(const matd_t *a)
matd_t * matd_select(const matd_t *a, int r0, int r1, int c0, int c1)
#define MATD_SVD_NO_WARNINGS
double matd_vec_dist_n(const matd_t *a, const matd_t *b, int n)
void matd_chol_destroy(matd_chol_t *chol)
void matd_utriangle_solve(matd_t *u, const TYPE *b, TYPE *x)
TYPE matd_get_scalar(const matd_t *m)
matd_t * matd_add(const matd_t *a, const matd_t *b)
void matd_ltriangle_solve(matd_t *L, const TYPE *b, TYPE *x)
matd_t * matd_chol_solve(const matd_chol_t *chol, const matd_t *b)
matd_t * matd_create_scalar(TYPE v)
void matd_put_scalar(matd_t *m, TYPE value)
#define MATD_EL(m, row, col)
matd_svd_t matd_svd(matd_t *A)
void matd_add_inplace(matd_t *a, const matd_t *b)
matd_chol_t * matd_chol(matd_t *A)
TYPE matd_get(const matd_t *m, int row, int col)
matd_t * matd_transpose(const matd_t *a)
static int matd_is_scalar(const matd_t *a)
matd_t * matd_vec_normalize(const matd_t *a)
matd_t * matd_solve(matd_t *A, matd_t *b)
static double sq(double v)
double matd_max(matd_t *m)
static int max_idx(const matd_t *A, int row, int maxcol)
matd_t * matd_create_data(int rows, int cols, const TYPE *data)
static int matd_is_vector(const matd_t *a)
static matd_t * matd_op_recurse(const char *expr, int *pos, matd_t *acc, matd_t **args, int *argpos, matd_t **garb, int *garbpos, int oneterm)
static double matd_det_general(const matd_t *a)
static matd_t * matd_op_gobble_right(const char *expr, int *pos, matd_t *acc, matd_t **garb, int *garbpos)
matd_t * matd_plu_u(const matd_plu_t *mlu)
void matd_scale_inplace(matd_t *a, double s)
matd_t * matd_scale(const matd_t *a, double s)
void matd_print_transpose(const matd_t *m, const char *fmt)
matd_t * matd_create_dataf(int rows, int cols, const float *data)
void matd_destroy(matd_t *m)
matd_t * matd_identity(int dim)
void matd_subtract_inplace(matd_t *a, const matd_t *b)
void matd_print(const matd_t *m, const char *fmt)
matd_t * matd_subtract(const matd_t *a, const matd_t *b)
matd_t * matd_create(int rows, int cols)
matd_t * matd_plu_solve(const matd_plu_t *mlu, const matd_t *b)
matd_t * matd_inverse(const matd_t *x)
void matd_put(matd_t *m, int row, int col, TYPE value)
matd_t * matd_plu_p(const matd_plu_t *mlu)
double matd_vec_dist(const matd_t *a, const matd_t *b)
double matd_vec_dot_product(const matd_t *a, const matd_t *b)
matd_t * matd_copy(const matd_t *m)
matd_t * matd_op(const char *expr,...)
void matd_ltransposetriangle_solve(matd_t *u, const TYPE *b, TYPE *x)
static matd_svd_t matd_svd_tall(matd_t *A, int flags)
matd_t * matd_multiply(const matd_t *a, const matd_t *b)
void matd_plu_destroy(matd_plu_t *mlu)
double matd_plu_det(const matd_plu_t *mlu)
matd_t * matd_crossproduct(const matd_t *a, const matd_t *b)
double matd_vec_mag(const matd_t *a)
static int matd_is_vector_len(const matd_t *a, int len)
matd_t * matd_chol_inverse(matd_t *a)
void svd22(const double A[4], double U[4], double S[2], double V[4])
TYPE matd_err_inf(const matd_t *a, const matd_t *b)