diff options
Diffstat (limited to 'sys/libmath/matrix.c')
| -rw-r--r-- | sys/libmath/matrix.c | 176 |
1 files changed, 0 insertions, 176 deletions
diff --git a/sys/libmath/matrix.c b/sys/libmath/matrix.c deleted file mode 100644 index e8bca0b..0000000 --- a/sys/libmath/matrix.c +++ /dev/null @@ -1,176 +0,0 @@ -#include <u.h> -#include <libn.h> -#include <libmath.h> - -/* TODO: replace (incrementally) with native C version! */ -#include <vendor/blas/cblas.h> -#include <vendor/blas/lapacke.h> - -// ----------------------------------------------------------------------- -// level 1 - -error -la·vecslice(math·Vector *x, int min, int max, int inc) -{ - if (max > x->len || min < 0) { - errorf("out of bounds: attempted to access vector past length"); - return 1; - } - x->len = (max - min) / inc; - x->d += x->inc * min; - x->inc *= inc; - - return 0; -} - -/* simple blas wrappers */ -void -la·veccopy(math·Vector *dst, math·Vector *src) -{ - return cblas_dcopy(src->len, src->d, src->inc, dst->d, dst->inc); -} - -double -la·vecnorm(math·Vector *x) -{ - return cblas_dnrm2(x->len, x->d, x->inc); -} - -void -la·vecscale(math·Vector *x, double a) -{ - return cblas_dscal(x->len, a, x->d, x->inc); -} - -double -la·vecdot(math·Vector *x, math·Vector *y) -{ - return cblas_ddot(x->len, x->d, x->inc, y->d, y->inc); -} - -// ----------------------------------------------------------------------- -// level 2 - -error -la·vecmat(math·Vector *x, math·Matrix *M) -{ - if (M->dim[1] != x->len) { - errorf("incompatible matrix dimensions"); - return 1; - } - if (M->state & ~mat·trans) - cblas_dgemv(CblasRowMajor,CblasNoTrans,M->dim[0],M->dim[1],1.,M->d,M->inc,x->d,x->inc,0.,x->d,x->inc); - else - cblas_dgemv(CblasRowMajor,CblasTrans,M->dim[0],M->dim[1],1.,M->d,M->inc,x->d,x->inc,0.,x->d,x->inc); - - return 0; -} - -// ----------------------------------------------------------------------- -// level 3 - -void -la·transpose(math·Matrix *X) -{ - int tmp; - X->state ^= mat·trans; - tmp = X->dim[0], X->dim[0] = X->dim[1], X->dim[1] = tmp; -} - -error -la·matrow(math·Matrix *X, int r, math·Vector *row) -{ - if (r < 0 || r >= X->dim[0]) { - errorf("out of bounds"); - return 1; - } - - row->len = X->dim[1]; - row->inc = 1; - row->d = X->d + X->dim[1] * r; - - return 0; -} - -error -la·matcol(math·Matrix *X, int c, math·Vector *col) -{ - if (c < 0 || c >= X->dim[1]) { - errorf("out of bounds"); - return 1; - } - - col->len = X->dim[0]; - col->inc = X->dim[1]; - col->d = X->d + c; - - return 0; -} - -error -la·matslice(math·Matrix *X, int r[3], int c[3]) -{ - /* TODO */ - return 0; -} - -error -la·eig(math·Matrix *X) -{ - -} - -/* X = A*B */ -error -la·matmul(math·Matrix *X, math·Matrix *A, math·Matrix *B) -{ - if (A->dim[1] != B->dim[0]) { - errorf("number of interior dimensions of A '%d' not equal to that of B '%d'", A->dim[1], B->dim[0]); - return 1; - } - if (X->dim[0] != A->dim[0]) { - errorf("number of exterior dimensions of X '%d' not equal to that of A '%d'", X->dim[0], A->dim[0]); - return 1; - } - if (X->dim[1] != B->dim[1]) { - errorf("number of exterior dimensions of X '%d' not equal to that of B '%d'", X->dim[1], B->dim[1]); - return 1; - } - - if (X->state & ~mat·trans) - if (A->state & ~mat·trans) - cblas_dgemm(CblasRowMajor,CblasNoTrans,CblasNoTrans,A->dim[0],B->dim[1],A->dim[1],1.,A->d,A->inc,B->d,B->inc,0.,X->d,X->inc); - else - cblas_dgemm(CblasRowMajor,CblasNoTrans,CblasTrans,A->dim[0],B->dim[1],A->dim[1],1.,A->d,A->inc,B->d,B->inc,0.,X->d,X->inc); - else - if (A->state & ~mat·trans) - cblas_dgemm(CblasRowMajor,CblasTrans,CblasNoTrans,A->dim[0],B->dim[1],A->dim[1],1.,A->d,A->inc,B->d,B->inc,0.,X->d,X->inc); - else - cblas_dgemm(CblasRowMajor,CblasTrans,CblasTrans,A->dim[0],B->dim[1],A->dim[1],1.,A->d,A->inc,B->d,B->inc,0.,X->d,X->inc); - - return 0; -} - -/* - * solves A*X=B - * pass in B via X - */ -error -la·solve(math·Matrix *X, math·Matrix *A) -{ - error err; - int n, *ipv; - static int buf[512]; - if (n = A->dim[0], n < arrlen(buf)) { - ipv = buf; - n = 0; - } else - ipv = malloc(n*sizeof(*ipv)); - - /* TODO: utilize more specific regimes if applicable */ - err = LAPACKE_dgesv(LAPACK_ROW_MAJOR,A->dim[0],X->dim[1],A->d,A->inc,ipv,X->d,X->inc); - - if (n) - free(ipv); - return err; -} |