lapack

Author: ktye

k+/blas.c

@@ -35,4 +35,4 @@ void blas(void){
  reg(rot,  "rot",4);
  reg(gemv,"gemv",7);
 }
-__attribute((section("reg")))void*rblas=blas;
+__attribute((section("rek")))void*rblas=blas;

k+/curl.c

@@ -20,4 +20,4 @@ uint64_t kurl(uint64_t x){
  dx(x);return r;}
 
 void curl(void){reg(kurl,"curl",1);}
-__attribute((section("reg")))void*rcurl=curl;
+__attribute((section("rek")))void*rcurl=curl;

k+/lapack.c

@@ -1,3 +1,48 @@
+int LAPACKE_dgesv(int,int,int,double*,int,int*,double*,int);
+int LAPACKE_zgesv(int,int,int,double*,int,int*,double*,int);
+int LAPACKE_dgels(int,char,int,int,int,double*,int,double*,int);
+int LAPACKE_zgels(int,char,int,int,int,double*,int,double*,int);
+int LAPACKE_dsyev(int,char,char,int,double*,int,double*);
+int LAPACKE_zheev(int,char,char,int,double*,int,double*);
+int LAPACKE_dgeev(int,char,char,int,double*,int,double*,double*,double*,int,double*,int);
+int LAPACKE_zgeev(int,char,char,int,double*,int,double*,        double*,int,double*,int);
+int LAPACKE_dgesvd(int,char,char,int,int,double*,int,double*,double*,int,double*,int,double*);
+int LAPACKE_zgesvd(int,char,char,int,int,double*,int,double*,double*,int,double*,int,double*);
+
+uint64_t gesv(uint64_t x){ //[A;B]
+ uint64_t A=FZ(rx(x0(x))),B=FZ(ati(x,1));
+ int32_t n=(int32_t)sqrt(nn(A));if(n*n!=nn(A))trap();
+ int32_t nb=nn(B),k=nb/n;if(nb!=k*n)trap();
+ int32_t t=maxtype(A,B);A=use(uptype(A,t));B=use(uptype(B,t));
+ int*p=malloc(4*n);int r;
+ if(t==zt)r=LAPACKE_zgesv(102,n,k,FK(A),n,p,FK(B),n);
+ else     r=LAPACKE_dgesv(102,n,k,FK(A),n,p,FK(B),n);
+ free(p);
+ if(r)trap();
+ dx(A);
+ return B;
+}
+uint64_t gels(uint64_t x){ //[rows;A;B]
+ uint64_t M=x0(x),A=FZ(rx(x1(x))),B=FZ(ati(x,2));
+ int32_t m=(int32_t)M;if(tp(M)!=it)trap();
+ int32_t n=nn(A)/m;if(n*m!=nn(A)||m<n)trap();
+ int32_t nb=nn(B),k=nb/m;if(n*k!=nb)trap();
+ int32_t r,t=maxtype(A,B);A=use(uptype(A,t));B=use(uptype(B,t));
+ if(t==zt)r=LAPACKE_zgels(102,'N',m,n,k,FK(A),m,FK(B),m);
+ else     r=LAPACKE_dgels(102,'N',m,n,k,FK(A),m,FK(B),m);
+ dx(A);return B;
+}
+void lapack(void){
+ reg(gesv,"gesv",2);
+ reg(gels,"gels",3);
+ //reg(geev,"geev",1);
+ //reg(geev,"geeV",1);
+ //reg(geev,"gesvd",1);
+ //reg(geev,"gesvD",1);
+}
+__attribute((section("rek")))void*rlapack=lapack;
+
+
 /*
 #include<stdlib.h>
 #include<string.h>

k+/readme

@@ -1,5 +1,5 @@
 # k+ extends k.c with some c-libraries using the native function type.
-# blas lapack fftw lp_solve ode specfun rmath curl sqlite..
+# blas lapack fftw lp_solve ode specfun curl sqlite..
 #
 # m.awk also replaces the builtin math functions with libm.
 #
@@ -19,14 +19,22 @@
 # dot dotc todo
 # gemv[a;A;rA;op;b;y;x] c:fz A:FZ b:fz y:FZ  a*A*x+b*y
 #
+# LAPACK
+# gesv[A;b]   or b,b,..  solve Ax=b, both col-major flat vectors
+# gels[n;A;b] or b,b,..  solve Ax=b least squares
+#
+# SPARSE  suitesparse:umfpack spqr 
+# spsolve[A;b] or (b;b;..) solve sparse system
+# spqr[A;b]    or  b,b,..  solve sparse least squares
+#
 # CURL
 # curl"https://ktye.github.io/k.png"
 #
-# RMATH (r-mathlib)
+# RMATH (r-mathlib standalone)
 
 # dependencies linux/windows, see deps
-# - linux(debian)    see deps
-# - windows(mingw64) see deps.$pacman -S mingw-w64-x86_64-lpsolve
+# - linux(debian13)  see deps
+# - windows(mingw64) see deps
 
 set -e
 set -x
@@ -43,10 +51,11 @@ cat       l.c>>k.c
 
 
 #these lines can be uncommented individually to exclude a package:
-#cat    blas.c>>k.c;L="$L -lopenblas"
-#cat    curl.c>>k.c;L="$L -lcurl"
+cat    blas.c>>k.c;L="$L -lopenblas"
+cat  lapack.c>>k.c;L="$L -llapacke"
+cat    curl.c>>k.c;L="$L -lcurl"
 cat  sparse.c>>k.c;L="$L -lumfpack -lspqr -lcholmod"
-#cat   rmath.c>>k.c;L="$L Rmath.dll"   #-lRmath
+#cat   rmath.c>>k.c;L="$L -lRmath"
 #              k.c
 
 #type help to see a the list of all registered functions and their arity

k+/rmath.c

@@ -6,7 +6,7 @@ rnorm(n, mean = 0, sd = 1)
 */
 
 //c/dat/r/R-4.5.1/src/include/Rmath.h
-double dnorm5(double,double,double,int);
+double dnorm4(double,double,double,int);
 double pnorm(double,double,double,int,int);
 double qnorm(double,double,double,int,int);
 double rnorm(double,double);
@@ -20,16 +20,15 @@ double rnorm(double,double);
 //t(2double) (student t dist)
 //cauchy
 //exp
+//..
 
-
-
-uint64_t dnrm(uint64_t x){x=fF(x);
- if(tp(x)<16){dx(x);return Kf(F64((int32_t)x));};
+uint64_t dnrm(uint64_t x){x=fF(Fst(x));
+ if(tp(x)<16){dx(x);return Kf(dnorm4(F64((int32_t)x),0,1,0));};
  x=use(x);int32_t e=ep(x);
- for(int32_t p=(int32_t)x;p<e;p+=8)SetF64(p,dnorm5(F64(p),0,1,0));
+ for(int32_t p=(int32_t)x;p<e;p+=8)SetF64(p,dnorm4(F64(p),0,1,0));
  return x;
 }
 void rmath(void){
  reg(dnrm,"dnorm",1);
 }
-__attribute((section("reg")))void*rrmath=rmath;
+__attribute((section("rek")))void*rrmath=rmath;

k+/sparse.c

@@ -1,4 +1,3 @@
-#include<suitesparse/SuiteSparseQR_C.h>
 //umfpack solve sparse Ax=b
 int  umfpack_di_symbolic(int,int,int*,int*,double*,        void**,double*,double*);
 int  umfpack_zi_symbolic(int,int,int*,int*,double*,double*,void**,double*,double*);
@@ -11,7 +10,6 @@ int  umfpack_zi_solve(int,int*,int*,double*,double*,double*,double*,double*,doub
 void umfpack_di_free_numeric(void**);
 void umfpack_zi_free_numeric(void**);
 
-/*
 //spqr sparse least squares (cholmod.h spqr.h)
 typedef struct choldmod_dense_struct {
  size_t nrow,ncol,nzmax,d; //A is nrow x ncol, nmax(entries) d(lda>=nrow)
@@ -25,8 +23,8 @@ typedef struct cholmod_sparse_struct{
 } cholmod_sparse;
 int cholmod_start(void*cc);
 int cholmod_finish(void*cc);
+int cholmod_free_dense(cholmod_dense**,void*);
 cholmod_dense*SuiteSparseQR_C_backslash_default(cholmod_sparse*A,cholmod_dense*B,void*cc);
-*/
 
 uint64_t spqr(uint64_t x){ //todo: b maybe m-by-k dense..
  uint64_t common[400];//enough to hold cholmod_common structure (sizeof 2680)
@@ -35,20 +33,22 @@ uint64_t spqr(uint64_t x){ //todo: b maybe m-by-k dense..
  uint64_t M=x0(A),N=x1(A),Ap=x2(A),Ai=x3(A),Ax=FZ(x4(A));dx(A);int z=(tp(Ax)==Zt);
  if(tp(M)!=it||tp(N)!=it||tp(Ap)!=It||tp(Ai)!=It)trap();
  int32_t m=(int32_t)M,n=(int32_t)N;
- if(nn(B)!=m)trap();if(z&&tp(B)!=Zt)B=uptype(B,zt);
+ if(nn(Ap)!=1+n||nn(Ai)!=nn(Ax))trap();
+ int32_t nb=nn(B),k=nb/m;if(nb!=k*m)trap(); //k>1: multirhs flat colmajor vector
+ if(z&&tp(B)!=Zt)B=uptype(B,zt);
  void *cc=&common;cholmod_start(cc);
  cholmod_sparse a;
- a.nrow=m;a.ncol=n;a.nzmax=nn(A);a.p=IK(Ap);a.i=IK(Ai);a.nz=NULL;a.x=FK(A);a.z=NULL;
+ a.nrow=m;a.ncol=n;a.nzmax=nn(Ax);a.p=IK(Ap);a.i=IK(Ai);a.nz=NULL;a.x=FK(Ax);a.z=NULL;
  a.stype=0;a.itype=0;a.xtype=1+z;a.dtype=0;a.sorted=0;a.packed=1;
  cholmod_dense b;
- b.nrow=m;b.ncol=1;b.nzmax=m;b.x=FK(B);b.z=NULL;b.xtype=1+z;b.dtype=0;
+ b.nrow=m;b.ncol=k;b.nzmax=m;b.x=FK(B);b.z=NULL;b.xtype=1+z;b.dtype=0;
  cholmod_dense*r=SuiteSparseQR_C_backslash_default(&a,&b,cc);
- //printf("r: %d %d #%d xt%d dt%d: [", r->nrow, r->ncol, r->nzmax, r->xtype, r->dtype);
- //printf("%g %g %g %g %g]\n", r->x[0], r->x[1], r->x[2], r->x[3], r->x[4]);
+ dx(Ap);dx(Ai);dx(Ax);dx(B);
+ x=mk(z?Zt:Ft,nb);
+ memcpy(M_+(int32_t)x,r->x,nb<<(z?4:3));
+ cholmod_free_dense(&r,cc);
  cholmod_finish(cc);
- //todo is A/B overwritten? where is the workspace, what must be freed?
- //cholmod_free_dense(
- return Ki(0);
+ return x;
 }
 
 //dense matrix representation
@@ -63,6 +63,7 @@ uint64_t spsolve(uint64_t x){ // spsolve[A;b]
  uint64_t N=x1(A),Ap=x2(A),Ai=x3(A),Ax=FZ(x4(A));dx(A);
  if(tp(N)!=it||tp(Ap)!=It||tp(Ai)!=It)trap();
  int z=0,n=(int32_t)N;int*ap=IK(Ap),*ai=IK(Ai);double*ax=FK(Ax);void*sy,*nu;
+ if(nn(Ap)!=1+n||nn(Ai)!=nn(Ax))trap();
  if(tp(Ax)==Ft){
   if(umfpack_di_symbolic(n,n,ap,ai,ax,&sy,NULL,NULL))trap();
   if(umfpack_di_numeric(ap,ai,ax,sy,&nu, NULL,NULL)){umfpack_di_free_symbolic(&sy);trap();}
@@ -96,7 +97,6 @@ uint64_t spsolve(uint64_t x){ // spsolve[A;b]
  return mrhs?r:Fst(r);
 }
 void umfpack(void){
- printf("%d\n", sizeof(cholmod_common)); //2680
  reg(spsolve,"spsolve",2);
  reg(spqr,   "spqr",   2);
 }
@@ -132,6 +132,6 @@ void *Symbolic, *Numeric ;
 #include<stdio.h>
 #include<cholmod.h>
 int main(){
- printf("%d\n", sizeof(cholmod_common)); //2680(win) 2664(debian)
+ printf("%d\n", sizeof(cholmod_common)); //should be 2680
 }
 */

k+/tests/sparse.t

@@ -3,4 +3,5 @@ b:8. 45. -3. 3. 19. /
 r:1 2 3 4 5. /
 1e-12>|/abs spsolve[A;b]-r /1
 _0.5++/+/spsolve[A;(b;b;b)] /45
-spqr[A;b] /3
+1e-12>|/abs spqr[A;b]-r /1
++/_0.1+spqr[A;b,b,b] /45