From 783a88f8904b5758195f96c65c7e0b6e87f8a51e Mon Sep 17 00:00:00 2001 From: cnlohr Date: Sat, 17 Mar 2018 01:00:07 -0400 Subject: Update DCLapack with the C helpers to cize stuff. --- redist/dclhelpers.c | 63 +++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 63 insertions(+) create mode 100644 redist/dclhelpers.c (limited to 'redist/dclhelpers.c') diff --git a/redist/dclhelpers.c b/redist/dclhelpers.c new file mode 100644 index 0000000..9d3d0eb --- /dev/null +++ b/redist/dclhelpers.c @@ -0,0 +1,63 @@ +#include "dclhelpers.h" +#define FLOAT DCL_FLOAT +#define DYNAMIC_INDEX +#include +#include "dclapack.h" + + +void dclPrint( const DCL_FLOAT * A, int n, int m ) +{ + PRINT( A, n, m ); +} + +void dclIdentity( DCL_FLOAT * A, int n ) +{ + IDENTITY( A, n ); +} + +void dclTransp( const DCL_FLOAT * A, DCL_FLOAT * B, int n, int m) +{ + TRANSP(A,B,n,m); +} + +void dclLU( const DCL_FLOAT * A, DCL_FLOAT * L, DCL_FLOAT * U, int * Piv, int n ) +{ + LU(A,L,U,Piv,n); +} + +void dclPivot( const DCL_FLOAT * A, DCL_FLOAT * B, int * Piv, int n, int m ) +{ + PIVOT(A,B,Piv,n,m); +} + +void dclLSub( const DCL_FLOAT * L, DCL_FLOAT * X, const DCL_FLOAT * B, int n, int m ) +{ + L_SUB(L,X,B,n,m); +} + +void dclUSub( const DCL_FLOAT * U, DCL_FLOAT * X, const DCL_FLOAT * B, int n, int m ) +{ + U_SUB(U,X,B,n,m); +} + +void dclInv( const DCL_FLOAT * A, DCL_FLOAT * Ainv, int n ) +{ + INV(A,Ainv,n,n); +} + +void dclMul( const DCL_FLOAT * A, const DCL_FLOAT * B, DCL_FLOAT * C, int n, int m, int p ) +{ + MUL(A,B,C,n,m,p); +} + +void dclMulAdd( const DCL_FLOAT * A, const DCL_FLOAT * B, const DCL_FLOAT * C, DCL_FLOAT * D, int n, int m, int p ) +{ + MULADD(A,B,C,D,n,m,p); +} + +void dclGMulAdd( const DCL_FLOAT * A, const DCL_FLOAT * B, const DCL_FLOAT * C, DCL_FLOAT * D, DCL_FLOAT alpha, DCL_FLOAT beta, int n, int m, int p ) +{ + GMULADD(A,B,C,D,alpha,beta,n,m,p); +} + + -- cgit v1.3.1 From 8b9a196232661bd6506a41b7ceca015d684086a8 Mon Sep 17 00:00:00 2001 From: cnlohr Date: Sat, 17 Mar 2018 01:13:28 -0400 Subject: Update variable ordering... --- redist/dclapack.h | 48 ++++++++++++++++++++++++------------------------ redist/dclhelpers.c | 36 ++++++++++++++++++------------------ redist/dclhelpers.h | 26 +++++++++++++------------- 3 files changed, 55 insertions(+), 55 deletions(-) (limited to 'redist/dclhelpers.c') diff --git a/redist/dclapack.h b/redist/dclapack.h index 62f3e9d..5aba214 100644 --- a/redist/dclapack.h +++ b/redist/dclapack.h @@ -44,15 +44,15 @@ } /* - * B = Transpose(A) + * R = Transpose(A) * A is (n by m) - * B is (m by n) + * R is (m by n) */ -#define TRANSP(A,B,n,m) { \ +#define TRANSP(R,A,n,m) { \ int i,j; \ for (i=0; i=0; j--) { \ @@ -160,14 +160,14 @@ FLOAT C[ORDER][ORDER]; #endif -#define INV(A,Ainv,n,ORDER) { \ +#define INV(Ainv,A,n,ORDER) { \ INV_SETUP(ORDER) \ int Piv[ORDER]; \ IDENTITY(I,n); \ - LU(A,L,U,Piv,n); \ - PIVOT(I,Ipiv,Piv,n,n); \ - L_SUB(L,C,Ipiv,n,n); \ - U_SUB(U,Ainv,C,n,n); \ + LU(L,U,A,Piv,n); \ + PIVOT(Ipiv,I,Piv,n,n); \ + L_SUB(C,L,Ipiv,n,n); \ + U_SUB(Ainv,U,C,n,n); \ } /* @@ -186,13 +186,13 @@ PRINT(Ainv,n,n); \ * B (m by p) * C (n by p) */ -#define MUL(A,B,C,n,m,p) { \ +#define MUL(R,A,B,n,m,p) { \ int i,j,k; \ for (i=0; i Date: Sat, 17 Mar 2018 02:28:11 -0400 Subject: Get closer to functional. --- redist/dclapack.h | 8 +++++--- redist/dclhelpers.c | 25 +++++++++++++++++++++++++ redist/dclhelpers.h | 24 ++++++++++++++++++++++++ 3 files changed, 54 insertions(+), 3 deletions(-) (limited to 'redist/dclhelpers.c') diff --git a/redist/dclapack.h b/redist/dclapack.h index bebda2d..76ce545 100644 --- a/redist/dclapack.h +++ b/redist/dclapack.h @@ -11,8 +11,10 @@ #ifdef DYNAMIC_INDEX #define _(A,O,P) A[O*n+P] + #define _I(A,O,P,n) A[O*n+P] #else #define _(A,O,P) A[O][P] + #define _I(A,O,P,n) A[O][P] #endif @@ -192,7 +194,7 @@ PRINT(Ainv,n,n); \ for (j=0; j Date: Sat, 17 Mar 2018 03:05:21 -0400 Subject: Now wrangled with ability to use submatricies. --- redist/dclapack.h | 24 +++++++++++------------- redist/dclhelpers.c | 46 ++++++++++++++++++++++++++++++---------------- redist/dclhelpers.h | 28 ++++++++++++++-------------- 3 files changed, 55 insertions(+), 43 deletions(-) (limited to 'redist/dclhelpers.c') diff --git a/redist/dclapack.h b/redist/dclapack.h index 76ce545..1dc7b77 100644 --- a/redist/dclapack.h +++ b/redist/dclapack.h @@ -10,16 +10,14 @@ #define _ABS(a) ( (a)<=0 ? 0-(a) : (a) ) #ifdef DYNAMIC_INDEX - #define _(A,O,P) A[O*n+P] - #define _I(A,O,P,n) A[O*n+P] + #define _(A,O,P) A[O*A##c+P] #else #define _(A,O,P) A[O][P] - #define _I(A,O,P,n) A[O][P] #endif /* - * Prints a matrix A (n by m) + * Prints a matrix A (n by m) (with witdth Ac) */ #define PRINT(A,n,m) { \ int i,j; \ @@ -34,7 +32,7 @@ } /* - * Returns the identity matrix + * Returns the identity matrix (with size n x n, but width Ic) */ #define IDENTITY(I,n) { \ int i,j; \ @@ -148,11 +146,11 @@ #ifdef DYNAMIC_INDEX #define INV_SETUP(ORDER) \ - FLOAT Ipiv[ORDER*ORDER]; \ - FLOAT L[ORDER*ORDER]; \ - FLOAT U[ORDER*ORDER]; \ - FLOAT I[ORDER*ORDER]; \ - FLOAT C[ORDER*ORDER]; + FLOAT Ipiv[ORDER*ORDER]; const int Ipivc = ORDER; \ + FLOAT L[ORDER*ORDER]; const int Lc = ORDER; \ + FLOAT U[ORDER*ORDER]; const int Uc = ORDER; \ + FLOAT I[ORDER*ORDER]; const int Ic = ORDER; \ + FLOAT C[ORDER*ORDER]; const int Cc = ORDER; #else #define INV_SETUP(ORDER) \ FLOAT Ipiv[ORDER][ORDER]; \ @@ -194,7 +192,7 @@ PRINT(Ainv,n,n); \ for (j=0; j #include "dclapack.h" +#include - -void dclPrint( const DCL_FLOAT * A, int n, int m ) +void dclPrint( const DCL_FLOAT * A, int Ac, int n, int m ) { PRINT( A, n, m ); } -void dclIdentity( DCL_FLOAT * A, int n ) +void dclIdentity( DCL_FLOAT * I, int Ic, int n ) { - IDENTITY( A, n ); + IDENTITY( I, n ); } -void dclTransp( DCL_FLOAT * R, const DCL_FLOAT * A, int n, int m ) +void dclTransp( DCL_FLOAT * R, int Rc, const DCL_FLOAT * A, int Ac, int n, int m ) { TRANSP(R,A,n,m); } -void dclLU( DCL_FLOAT * L, DCL_FLOAT * U, const DCL_FLOAT * A, int * Piv, int n ) +void dclLU( DCL_FLOAT * L, int Lc, DCL_FLOAT * U, int Uc, const DCL_FLOAT * A, int Ac, int * Piv, int n ) { LU(L,U,A,Piv,n); } -void dclPivot( DCL_FLOAT * R, const DCL_FLOAT * A, int * Piv, int n, int m ) +void dclPivot( DCL_FLOAT * R, int Rc, const DCL_FLOAT * A, int Ac, int * Piv, int n, int m ) { PIVOT(R,A,Piv,n,m); } -void dclLSub( DCL_FLOAT * X, const DCL_FLOAT * L, const DCL_FLOAT * B, int n, int m ) +void dclLSub( DCL_FLOAT * X, int Xc, const DCL_FLOAT * L, int Lc, const DCL_FLOAT * B, int Bc, int n, int m ) { L_SUB(X,L,B,n,m); } -void dclUSub( DCL_FLOAT * X, const DCL_FLOAT * U, const DCL_FLOAT * B, int n, int m ) +void dclUSub( DCL_FLOAT * X, int Xc, const DCL_FLOAT * U, int Uc, const DCL_FLOAT * B, int Bc, int n, int m ) { U_SUB(X,U,B,n,m); } -void dclInv( DCL_FLOAT * Ainv, const DCL_FLOAT * A, int n ) +void dclInv( DCL_FLOAT * Ainv, int Ainvc, const DCL_FLOAT * A, int Ac, int n ) { INV(Ainv,A,n,n); } -void dclMul( DCL_FLOAT * R, const DCL_FLOAT * A, const DCL_FLOAT * B, int n, int m, int p ) +void dclMul( DCL_FLOAT * R, int Rc, const DCL_FLOAT * A, int Ac, const DCL_FLOAT * B, int Bc, int n, int m, int p ) { MUL(R,A,B,n,m,p); } -void dclMulAdd( DCL_FLOAT * R, const DCL_FLOAT * A, const DCL_FLOAT * B, const DCL_FLOAT * C, int n, int m, int p ) +void dclMulAdd( DCL_FLOAT * R, int Rc, const DCL_FLOAT * A, int Ac, const DCL_FLOAT * B, int Bc, const DCL_FLOAT * C, int Cc, int n, int m, int p ) { MULADD(R,A,B,C,n,m,p); } -void dclGMulAdd( DCL_FLOAT * R, const DCL_FLOAT * A, const DCL_FLOAT * B, const DCL_FLOAT * C, DCL_FLOAT alpha, DCL_FLOAT beta, int n, int m, int p ) +void dclGMulAdd( DCL_FLOAT * R, int Rc, const DCL_FLOAT * A, int Ac, const DCL_FLOAT * B, int Bc, const DCL_FLOAT * C, int Cc, DCL_FLOAT alpha, DCL_FLOAT beta, int n, int m, int p ) { GMULADD(R,A,B,C,alpha,beta,n,m,p); } @@ -77,12 +77,26 @@ void dcldgemm( int ldc //must be n ) { - DCL_FLOAT * ta; - DCL_FLOAT * tb; + const DCL_FLOAT * ta; + const DCL_FLOAT * tb; if( transA ) { ta = alloca( sizeof( DCL_FLOAT ) * n * m ); - + //TRANSP( ta, A, n, m ); + } + else + ta = A; + + if( transB ) + { + tb = alloca( sizeof( DCL_FLOAT ) * n * m ); + //TRANSP( tb, B, n, m ); } + else + tb = B; + + + //XXX Tricky: In here, A is mxn + } diff --git a/redist/dclhelpers.h b/redist/dclhelpers.h index fd4e02d..5010e02 100644 --- a/redist/dclhelpers.h +++ b/redist/dclhelpers.h @@ -6,56 +6,56 @@ //XXX XXX XXX WARNING XXX XXX XXX The argument order may be changing!!! /* Prints matrix A of size[n][m] */ -void dclPrint( const DCL_FLOAT * A, int n, int m ); +void dclPrint( const DCL_FLOAT * A, int Ac, int n, int m ); /* Returns the identity matrix */ -void dclIdentity( DCL_FLOAT * A, int n ); +void dclIdentity( DCL_FLOAT * I, int Ic, int n ); /* R = Transpose(A) A is (n by m) R is (m by n) */ -void dclTransp( DCL_FLOAT * R, const DCL_FLOAT * A, int n, int m ); +void dclTransp( DCL_FLOAT * R, int Rc, const DCL_FLOAT * A, int Ac, int n, int m ); /* Calculate L,U of a matrix A with pivot table; the pivot table is output. */ -void dclLU( DCL_FLOAT * L, DCL_FLOAT * U, const DCL_FLOAT * A, int * Piv, int n ); +void dclLU( DCL_FLOAT * L, int Lc, DCL_FLOAT * U, int Uc, const DCL_FLOAT * A, int Ac, int * Piv, int n ); /* Pivots a matrix to a different matrix R = Pivot(A) given table 'Piv' A and R are (n by m) */ -void dclPivot( DCL_FLOAT * R, const DCL_FLOAT * A, int * Piv, int n, int m ); +void dclPivot( DCL_FLOAT * R, int Rc, const DCL_FLOAT * A, int Ac, int * Piv, int n, int m ); /* Solve LX=B for matrix X and B L is n by n (lower triangular) B is n by m */ -void dclLSub( DCL_FLOAT * X, const DCL_FLOAT * L, const DCL_FLOAT * B, int n, int m ); +void dclLSub( DCL_FLOAT * X, int Xc, const DCL_FLOAT * L, int Lc, const DCL_FLOAT * B, int Bc, int n, int m ); /* Solve UX=B for matrix X and B U is n by n (upper triangular) B is n by m */ -void dclUSub( DCL_FLOAT * X, const DCL_FLOAT * U, const DCL_FLOAT * B, int n, int m ); +void dclUSub( DCL_FLOAT * X, int Xc, const DCL_FLOAT * U, int Uc, const DCL_FLOAT * B, int Bc, int n, int m ); /* Inverts a matrix X (n by n) using the method of LU decomposition */ -void dclInv( DCL_FLOAT * Ainv, const DCL_FLOAT * A, int n ); +void dclInv( DCL_FLOAT * Ainv, int Ainvc, const DCL_FLOAT * A, int Ac, int n ); /* Matrix Multiply R = A * B A (n by m) B (m by p) R (n by p) */ -void dclMul( DCL_FLOAT * R, const DCL_FLOAT * A, const DCL_FLOAT * B, int n, int m, int p ); +void dclMul( DCL_FLOAT * R, int Rc, const DCL_FLOAT * A, int Ac, const DCL_FLOAT * B, int Bc, int n, int m, int p ); /* Matrix Multiply R = A * B + C A (n by m) B (m by p) C (n by p) R (n by p) */ -void dclMulAdd( DCL_FLOAT * R, const DCL_FLOAT * A, const DCL_FLOAT * B, const DCL_FLOAT * C, int n, int m, int p ); +void dclMulAdd( DCL_FLOAT * R, int Rc, const DCL_FLOAT * A, int Ac, const DCL_FLOAT * B, int Bc, const DCL_FLOAT * C, int Cc, int n, int m, int p ); /* Matrix Multiply R = alpha * A * B + beta * C A (n by m) B (m by p) C (n by p) R (n by p) */ -void dclGMulAdd( DCL_FLOAT * R, const DCL_FLOAT * A, const DCL_FLOAT * B, const DCL_FLOAT * C, DCL_FLOAT alpha, DCL_FLOAT beta, int n, int m, int p ); +void dclGMulAdd( DCL_FLOAT * R, int Rc, const DCL_FLOAT * A, int Ac, const DCL_FLOAT * B, int Bc, const DCL_FLOAT * C, int Cc, DCL_FLOAT alpha, DCL_FLOAT beta, int n, int m, int p ); /******************************** @@ -71,12 +71,12 @@ void dcldgemm( int k, DCL_FLOAT alpha, const DCL_FLOAT* A, - int lda, //must be n + int Ac, const DCL_FLOAT* B, - int ldb, //must be m + int Bc, DCL_FLOAT beta, const DCL_FLOAT * C, - int ldc //must be n + int Cc ); -- cgit v1.3.1 From f0c26bd1b0b8ffe05c0c5f04567a9b7aa47c3e6b Mon Sep 17 00:00:00 2001 From: cnlohr Date: Sat, 17 Mar 2018 03:15:46 -0400 Subject: Update dclapack to do automatic size-of-array'ing --- redist/dclapack.h | 8 ++++++-- redist/dclhelpers.c | 28 +++++++++++++++++----------- redist/dclhelpers.h | 4 +--- 3 files changed, 24 insertions(+), 16 deletions(-) (limited to 'redist/dclhelpers.c') diff --git a/redist/dclapack.h b/redist/dclapack.h index 1dc7b77..af8869c 100644 --- a/redist/dclapack.h +++ b/redist/dclapack.h @@ -9,10 +9,14 @@ #define _ABS(a) ( (a)<=0 ? 0-(a) : (a) ) +//Tricky: If you want to use this with pointers, instead of 2D arrays, you will +//need to #define DYNAMIC_INDEX, as well as, for all arrays, suffix their name +//with 'c' + #ifdef DYNAMIC_INDEX - #define _(A,O,P) A[O*A##c+P] + #define _(AM,O,P) AM[O*AM##c+P] #else - #define _(A,O,P) A[O][P] + #define _(AM,O,P) AM[O][P] #endif diff --git a/redist/dclhelpers.c b/redist/dclhelpers.c index 9fa9266..c3725f8 100644 --- a/redist/dclhelpers.c +++ b/redist/dclhelpers.c @@ -69,34 +69,40 @@ void dcldgemm( int k, DCL_FLOAT alpha, const DCL_FLOAT* A, - int lda, //must be n + int Ac, //must be n const DCL_FLOAT* B, - int ldb, //must be m + int Bc, //must be m DCL_FLOAT beta, - const DCL_FLOAT * C, - int ldc //must be n + DCL_FLOAT * C, + int Cc //must be n ) { const DCL_FLOAT * ta; const DCL_FLOAT * tb; + int tac = Ac; + int tbc = Bc; if( transA ) { - ta = alloca( sizeof( DCL_FLOAT ) * n * m ); - //TRANSP( ta, A, n, m ); + DCL_FLOAT * la = alloca( sizeof( DCL_FLOAT ) * n * m ); + const int lac = m; + TRANSP( la, A, n, m ); + ta = la; + tac = lac; } else ta = A; if( transB ) { - tb = alloca( sizeof( DCL_FLOAT ) * n * m ); - //TRANSP( tb, B, n, m ); + DCL_FLOAT * lb = alloca( sizeof( DCL_FLOAT ) * n * m ); + const int lbc = m; + TRANSP( lb, B, n, m ); + tb = lb; + tbc = lbc; } else tb = B; - - //XXX Tricky: In here, A is mxn - + GMULADD(C,ta,tb,C,alpha,beta,n,m,k); } diff --git a/redist/dclhelpers.h b/redist/dclhelpers.h index 5010e02..8779e1a 100644 --- a/redist/dclhelpers.h +++ b/redist/dclhelpers.h @@ -3,8 +3,6 @@ #define DCL_FLOAT FLT -//XXX XXX XXX WARNING XXX XXX XXX The argument order may be changing!!! - /* Prints matrix A of size[n][m] */ void dclPrint( const DCL_FLOAT * A, int Ac, int n, int m ); @@ -75,7 +73,7 @@ void dcldgemm( const DCL_FLOAT* B, int Bc, DCL_FLOAT beta, - const DCL_FLOAT * C, + DCL_FLOAT * C, int Cc ); -- cgit v1.3.1 From 9115ffd3138b460707dd1ba45dd7f6fccde87a46 Mon Sep 17 00:00:00 2001 From: cnlohr Date: Sat, 17 Mar 2018 03:28:23 -0400 Subject: Test DCL. --- redist/dclhelpers.c | 4 ++-- redist/test_dcl.c | 39 +++++++++++++++++++++++++++++++++++++++ 2 files changed, 41 insertions(+), 2 deletions(-) create mode 100644 redist/test_dcl.c (limited to 'redist/dclhelpers.c') diff --git a/redist/dclhelpers.c b/redist/dclhelpers.c index c3725f8..fb6aba6 100644 --- a/redist/dclhelpers.c +++ b/redist/dclhelpers.c @@ -5,9 +5,9 @@ #include "dclapack.h" #include -void dclPrint( const DCL_FLOAT * A, int Ac, int n, int m ) +void dclPrint( const DCL_FLOAT * PMATRIX, int PMATRIXc, int n, int m ) { - PRINT( A, n, m ); + PRINT( PMATRIX, n, m ); } void dclIdentity( DCL_FLOAT * I, int Ic, int n ) diff --git a/redist/test_dcl.c b/redist/test_dcl.c new file mode 100644 index 0000000..adea7b5 --- /dev/null +++ b/redist/test_dcl.c @@ -0,0 +1,39 @@ +#include "dclhelpers.h" +#include +#include + + +int main() +{ + FLT A[2][4] = { { 0, 1, 2, 3 }, { 4, 5, 6, 7} }; + FLT B[4][2]; + dclPrint( A[0], 4, 2, 4 ); + dclTransp( B[0], 2, A[0], 4, 2, 4 ); + dclPrint( B[0], 2, 4, 2 ); + + int i; + for( i = 0; i < 8; i++ ) + { + printf( "%f\n", ((float*)(B[0]))[i] ); + } + + FLT M[3][3] = { + { .32, 1, 0 }, + { 0, 1, 2 }, + { 1, 0, 1 } }; + FLT Mo[3][3]; + dclInv( Mo[0], 3, M[0], 3, 3 ); + dclPrint( Mo[0], 3, 3, 3 ); + + FLT MM[3][3]; + dclMul( MM[0], 3, M[0], 3, Mo[0], 3, 3, 3, 3 ); + + printf( "The following should be an identity matrix\n" ); + dclPrint( MM[0], 3, 3, 3 ); + +//void dclTransp( DCL_FLOAT * R, int Rc, const DCL_FLOAT * A, int Ac, int n, int m ); + +// dclIdentity( A[0], MATx, 5 ); +// dclPrint( A[0], MATx, MATx, MATy ); +} + -- cgit v1.3.1 From 0b9e66ad2ff686a4dcf8a6838f33edb203a1bff5 Mon Sep 17 00:00:00 2001 From: Justin Berger Date: Sat, 17 Mar 2018 09:32:22 -0600 Subject: Fixed gemm --- Makefile | 8 +++++++- redist/dclapack.h | 33 +++++++++++++++++---------------- redist/dclhelpers.c | 9 +++++---- redist/test_dcl.c | 23 +++++++++++++++-------- 4 files changed, 44 insertions(+), 29 deletions(-) (limited to 'redist/dclhelpers.c') diff --git a/Makefile b/Makefile index 9239c39..dfcb9f4 100644 --- a/Makefile +++ b/Makefile @@ -82,7 +82,13 @@ calibrate_client : calibrate_client.c ./lib/libsurvive.so redist/os_generic.c $ static_calibrate : calibrate.c redist/os_generic.c $(DRAWFUNCTIONS) $(LIBSURVIVE_C) tcc -o $@ $^ $(CFLAGS) $(LDFLAGS) -DTCC -test_dcl: ./redist/test_dcl.c ./redist/dclhelpers.c +./redist/dclhelpers_debuggable.c : ./redist/dclhelpers.c ./redist/dclhelpers.h ./redist/dclapack.h + gcc -E ./redist/dclhelpers.c > ./redist/dclhelpers_debuggable.c + clang-format -i ./redist/dclhelpers_debuggable.c + sed -i 's/#/\/\/#/g' ./redist/dclhelpers_debuggable.c + + +test_dcl: ./redist/test_dcl.c ./redist/dclhelpers_debuggable.c ./redist/dclhelpers.h ./redist/dclapack.h $(CC) -o $@ $^ $(LDFLAGS) $(CFLAGS) -DFLT=double -fsanitize=address -fsanitize=undefined test_minimal_cv: ./src/epnp/test_minimal_cv.c ./lib/libsurvive.so diff --git a/redist/dclapack.h b/redist/dclapack.h index af8869c..af5035b 100644 --- a/redist/dclapack.h +++ b/redist/dclapack.h @@ -224,22 +224,23 @@ PRINT(Ainv,n,n); \ /* * Matrix Multiply R = alpha * A * B + beta * C * R (n by p) - * A (n by m) - * B (m by p) - * C (n by p) + * A (m by n) + * B (n by p) + * C (m by p) */ -#define GMULADD(R,A,B,C,alpha,beta,n,m,p) { \ - int i,j,k; \ - float sum; \ - for (i=0; i #include "dclapack.h" #include +#include +#include void dclPrint( const DCL_FLOAT * PMATRIX, int PMATRIXc, int n, int m ) { @@ -77,7 +78,7 @@ void dcldgemm( int Cc //must be n ) { - const DCL_FLOAT * ta; + const DCL_FLOAT *ta; const DCL_FLOAT * tb; int tac = Ac; int tbc = Bc; @@ -102,7 +103,7 @@ void dcldgemm( } else tb = B; - - GMULADD(C,ta,tb,C,alpha,beta,n,m,k); + printf("%d %d %d\n", tac, tbc, Cc); + GMULADD(C, ta, tb, C, alpha, beta, m, n, k); } diff --git a/redist/test_dcl.c b/redist/test_dcl.c index 6d49548..42f4fd6 100644 --- a/redist/test_dcl.c +++ b/redist/test_dcl.c @@ -1,5 +1,6 @@ #include "dclhelpers.h" #include +#include #include #include @@ -36,18 +37,24 @@ int main() dclIdentity(A[0], 4, 3); dclPrint(A[0], 4, 3, 4); - FLT x[4] = {7, 8, 9, 10}; - FLT R[4]; + FLT x[4][2] = { + {7, -7}, {8, -8}, {9, -9}, {10, -10}, + }; + FLT R[4][2]; + printf("%p %p %p\n", A, x, R); // dclMul(R, 1, A[0], 4, x, 1, 4, 1, 3); - dcldgemm(0, 0, 4, 1, 3, 1, A[0], 4, x, 1, 0, R, 1); + dcldgemm(0, 0, 3, 4, 2, 1, A[0], 4, x[0], 2, 0, R[0], 2); - dclPrint(x, 1, 4, 1); - dclPrint(R, 1, 4, 1); + dclPrint(x[0], 2, 4, 2); + dclPrint(R[0], 2, 4, 2); - for (int i = 0; i < 3; i++) - assert(R[i] == x[i]); - assert(R[3] == 0.); + for (int j = 0; j < 2; j++) { + for (int i = 0; i < 3; i++) + assert(R[i][j] == x[i][j]); + + assert(fabs(R[3][j]) < .0000001); + } } // void dclTransp( DCL_FLOAT * R, int Rc, const DCL_FLOAT * A, int Ac, int n, int m ); -- cgit v1.3.1 From 04bd16aeb391e67716344268cf0f43d1f31f180a Mon Sep 17 00:00:00 2001 From: cnlohr Date: Sat, 17 Mar 2018 14:21:20 -0400 Subject: Update dcl and test. --- redist/dclapack.h | 20 ++++++++++++++++---- redist/dclhelpers.c | 14 +++++++++++--- redist/dclhelpers.h | 8 +++++++- redist/test_dcl.c | 41 ++++++++++++++++++++++++++--------------- 4 files changed, 60 insertions(+), 23 deletions(-) (limited to 'redist/dclhelpers.c') diff --git a/redist/dclapack.h b/redist/dclapack.h index 2823c5d..cae377b 100644 --- a/redist/dclapack.h +++ b/redist/dclapack.h @@ -37,9 +37,9 @@ /* * Returns the identity matrix (with size n x n, but width Ic) */ -#define IDENTITY(I, n) \ +#define IDENTITY(I, m, n) \ { \ - for (int _i = 0; _i < (n); _i++) { \ + for (int _i = 0; _i < (m); _i++) { \ for (int _j = 0; _j < _i; _j++) { \ _(I, _i, _j) = 0.0f; \ } \ @@ -50,6 +50,17 @@ } \ } + +/* + * Returns the identity matrix (with size n x n, but width Ic) + */ +#define ZERO(Z, m, n) \ + { \ + for (int _i = 0; _i < (m); _i++) \ + for (int _j = 0; _j < (n); _j++) \ + _(Z, _i, _j) = 0.0f; \ + } + /* * R = Transpose(A) * A is (m by n) @@ -79,7 +90,7 @@ _(U, _i, _j) = _(A, _i, _j); \ } \ } \ - IDENTITY(L, n); \ + IDENTITY(L, n, n); \ \ for (_i = 0; _i < (n)-1; _i++) { \ \ @@ -188,7 +199,7 @@ #define INV(Ainv,A,n,ORDER) { \ INV_SETUP(ORDER) \ int Piv[ORDER]; \ - IDENTITY(I,n); \ + IDENTITY(I,n,n); \ LU(L,U,A,Piv,n); \ PIVOT(Ipiv,I,Piv,n,n); \ L_SUB(C,L,Ipiv,n,n); \ @@ -258,6 +269,7 @@ PRINT(Ainv,n,n); \ for (int _k = 0; _k < n; _k++) { \ sum += _(A, _i, _k) * _(B, _k, _j); \ } \ + printf( "-> %d %d = %f %f %f %f[%d %d]\n", _i, _j, alpha,sum,beta, _(C,_i,_j), _i, _j ); \ _(R, _i, _j) = alpha * sum + beta * _(C, _i, _j); \ } \ } \ diff --git a/redist/dclhelpers.c b/redist/dclhelpers.c index 3e51fd2..5b956b0 100644 --- a/redist/dclhelpers.c +++ b/redist/dclhelpers.c @@ -5,15 +5,23 @@ #include #include #include +#include void dclPrint( const DCL_FLOAT * PMATRIX, int PMATRIXc, int n, int m ) { PRINT( PMATRIX, n, m ); } -void dclIdentity( DCL_FLOAT * I, int Ic, int n ) +void dclIdentity( DCL_FLOAT * I, int Ic, int m, int n ) { - IDENTITY( I, n ); + IDENTITY( I, m, n ); +} + + +/* Returns the zero matrix */ +void dclZero( DCL_FLOAT * Z, int Zc, int m, int n ) +{ + memset( Z, 0, m*n*sizeof(DCL_FLOAT) ); } void dclTransp( DCL_FLOAT * R, int Rc, const DCL_FLOAT * A, int Ac, int n, int m ) @@ -103,7 +111,7 @@ void dcldgemm( } else tb = B; - printf("%d %d %d\n", tac, tbc, Cc); + GMULADD(C, ta, tb, C, alpha, beta, m, n, k); } diff --git a/redist/dclhelpers.h b/redist/dclhelpers.h index 8779e1a..b4acec8 100644 --- a/redist/dclhelpers.h +++ b/redist/dclhelpers.h @@ -3,11 +3,17 @@ #define DCL_FLOAT FLT +//Use this macro to safely +#define DMS( m ) ((m)[0]), (sizeof((m)[0])/sizeof((m)[0][0])) + /* Prints matrix A of size[n][m] */ void dclPrint( const DCL_FLOAT * A, int Ac, int n, int m ); /* Returns the identity matrix */ -void dclIdentity( DCL_FLOAT * I, int Ic, int n ); +void dclIdentity( DCL_FLOAT * I, int Ic, int m, int n ); + +/* Returns the zero matrix */ +void dclZero( DCL_FLOAT * I, int Ic, int m, int n ); /* R = Transpose(A) A is (n by m) diff --git a/redist/test_dcl.c b/redist/test_dcl.c index 056acdf..6b0d870 100644 --- a/redist/test_dcl.c +++ b/redist/test_dcl.c @@ -12,7 +12,7 @@ int main() dclTransp( B[0], 2, A[0], 4, 2, 4 ); dclPrint( B[0], 2, 4, 2 ); - int i; + int i, j; for( i = 0; i < 8; i++ ) { printf( "%f\n", ((float*)(B[0]))[i] ); @@ -34,54 +34,65 @@ int main() { FLT A[3][4]; - dclIdentity(A[0], 4, 3); - dclPrint(A[0], 4, 3, 4); + dclIdentity( DMS(A), 3, 4); + dclPrint( DMS(A), 3, 4); FLT x[4][2] = { {7, -7}, {8, -8}, {9, -9}, {10, -10}, }; FLT R[4][2]; + dclZero( DMS(R), 4, 2 ); - printf("%p %p %p\n", A, x, R); // dclMul(R, 1, A[0], 4, x, 1, 4, 1, 3); dcldgemm(0, 0, 3, 4, 2, 1, A[0], 4, x[0], 2, 0, R[0], 2); - dclPrint(x[0], 2, 4, 2); - dclPrint(R[0], 2, 4, 2); + dclPrint(DMS(x), 4, 2); + dclPrint(DMS(R), 3, 2); for (int j = 0; j < 2; j++) { for (int i = 0; i < 3; i++) + { + printf( "[%d][%d]\n", i, j ); assert(R[i][j] == x[i][j]); + } assert(fabs(R[3][j]) < .0000001); } } - printf( "The following should be an identity matrix\n" ); - dclPrint( MM[0], 3, 3, 3 ); - +#if 0 //Currently failing test... { +// FLT em1[3][4]; +// FLT em2[4][2]; +// FLT emo[4][2]; + FLT em1[12][20]; FLT em2[20][12]; - FLT emo[12][12]; + FLT emo[20][12]; int x, y; for( y = 0; y < 12; y++ ) for( x = 0; x < 20; x++ ) - { em1[y][x] = (rand()%1000)/1000.0; - em2[x][y] = (rand()%1000)/1000.0; - } + + for( y = 0; y < 20; y++ ) + for( x = 0; x < 12; x++ ) + em2[y][x] = (rand()%1000)/1000.0; int m = 12; int n = 20; int k = 12; - dcldgemm( 0, 0, m, n, k, 0.1, em1[0], 20, em2[0], 12, .1, emo[0], 12 ); - dclPrint( emo[0], 12, 12, 12 ); + + dclPrint( em1[0], 20, 12, 20 ); + dclPrint( em2[0], 12, 20, 12 ); + + dcldgemm( 0, 0, m, n, k, 1.0, em1[0], 20, em2[0], 12, .1, emo[0], 12 ); + dclPrint( emo[0], 12, 20, 12 ); } +#endif } -- cgit v1.3.1