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/dclapack.h | 78 +++++++++++++++++++++++++++++++++---------------------- 1 file changed, 47 insertions(+), 31 deletions(-) (limited to 'redist/dclapack.h') diff --git a/redist/dclapack.h b/redist/dclapack.h index 4e209d3..62f3e9d 100644 --- a/redist/dclapack.h +++ b/redist/dclapack.h @@ -1,10 +1,6 @@ #ifndef __DCLAPACK_H__ #define __DCLAPACK_H__ -#ifndef ORDER -#define ORDER 50 -#endif - #ifndef FLOAT #define FLOAT float #endif @@ -13,6 +9,13 @@ #define _ABS(a) ( (a)<=0 ? 0-(a) : (a) ) +#ifdef DYNAMIC_INDEX + #define _(A,O,P) A[O*n+P] +#else + #define _(A,O,P) A[O][P] +#endif + + /* * Prints a matrix A (n by m) */ @@ -21,7 +24,7 @@ printf(#A "\n"); \ for (i=0; i _ABS(U[max][i])) { max = j; } \ + if (_ABS( _(U,j,i)) > _ABS(_(U,max,i))) { max = j; } \ } \ _tempi=Piv[i]; Piv[i]=Piv[max]; Piv[max]=_tempi; \ for (k=i; k=0; j--) { \ float sum=0.0; \ - for (k=n-1; k>j; k--) { sum += U[j][k]*X[k][i]; } \ - X[j][i] = (B[j][i] - sum) / U[j][j]; \ + for (k=n-1; k>j; k--) { sum += _(U,j,k)*_(X,k,i); } \ + _(X,j,i) = (_(B,j,i) - sum) / _(U,j,j); \ } \ } \ } @@ -140,12 +143,25 @@ /* * Inverts a matrix X (n by n) using the method of LU decomposition */ -#define INV(A,Ainv,n) { \ + +#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]; +#else + #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 C[ORDER][ORDER]; +#endif + +#define INV(A,Ainv,n,ORDER) { \ + INV_SETUP(ORDER) \ int Piv[ORDER]; \ IDENTITY(I,n); \ LU(A,L,U,Piv,n); \ @@ -174,9 +190,9 @@ PRINT(Ainv,n,n); \ int i,j,k; \ for (i=0; i Date: Sat, 17 Mar 2018 01:13:28 -0400 Subject: Update variable ordering... --- redist/dclapack.h | 48 ++++++++++++++++++++++++------------------------ 1 file changed, 24 insertions(+), 24 deletions(-) (limited to 'redist/dclapack.h') 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 01:34:36 -0400 Subject: Fix comments in dclapack code to reflect current behavior --- redist/dclapack.h | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) (limited to 'redist/dclapack.h') diff --git a/redist/dclapack.h b/redist/dclapack.h index 5aba214..bebda2d 100644 --- a/redist/dclapack.h +++ b/redist/dclapack.h @@ -181,10 +181,10 @@ PRINT(Ainv,n,n); \ */ /* - * Matrix Multiply C = A * B + * Matrix Multiply R = A * B + * R (n by p) * A (n by m) * B (m by p) - * C (n by p) */ #define MUL(R,A,B,n,m,p) { \ int i,j,k; \ @@ -199,11 +199,11 @@ PRINT(Ainv,n,n); \ } /* - * Matrix Multiply D = A * B + C + * Matrix Multiply R = A * B + C + * R (n by p) * A (n by m) * B (m by p) * C (n by p) - * D (n by p) */ #define MULADD(R,A,B,C,n,m,p) { \ int i,j,k; \ @@ -218,11 +218,11 @@ PRINT(Ainv,n,n); \ } /* - * Matrix Multiply D = alpha * A * B + beta * C + * Matrix Multiply R = alpha * A * B + beta * C + * R (n by p) * A (n by m) * B (m by p) * C (n by p) - * D (n by p) */ #define GMULADD(R,A,B,C,alpha,beta,n,m,p) { \ int i,j,k; \ -- cgit v1.2.3 From bcf08b95ab6daa7ac7bffe1449fa8a11cad2a02a Mon Sep 17 00:00:00 2001 From: cnlohr Date: Sat, 17 Mar 2018 02:28:11 -0400 Subject: Get closer to functional. --- redist/dclapack.h | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) (limited to 'redist/dclapack.h') 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 +++++++++++------------- 1 file changed, 11 insertions(+), 13 deletions(-) (limited to 'redist/dclapack.h') 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 Date: Sat, 17 Mar 2018 03:15:46 -0400 Subject: Update dclapack to do automatic size-of-array'ing --- redist/dclapack.h | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) (limited to 'redist/dclapack.h') 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 -- cgit v1.2.3 From 0b9e66ad2ff686a4dcf8a6838f33edb203a1bff5 Mon Sep 17 00:00:00 2001 From: Justin Berger Date: Sat, 17 Mar 2018 09:32:22 -0600 Subject: Fixed gemm --- redist/dclapack.h | 33 +++++++++++++++++---------------- 1 file changed, 17 insertions(+), 16 deletions(-) (limited to 'redist/dclapack.h') 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 Date: Sat, 17 Mar 2018 10:33:14 -0600 Subject: Made naming follow conventions --- redist/dclapack.h | 278 +++++++++++++++++++++++++++++------------------------- 1 file changed, 149 insertions(+), 129 deletions(-) (limited to 'redist/dclapack.h') diff --git a/redist/dclapack.h b/redist/dclapack.h index af5035b..2823c5d 100644 --- a/redist/dclapack.h +++ b/redist/dclapack.h @@ -19,130 +19,151 @@ #define _(AM,O,P) AM[O][P] #endif - /* - * Prints a matrix A (n by m) (with witdth Ac) + * Prints a matrix A (m by n) (with witdth Ac) */ -#define PRINT(A,n,m) { \ - int i,j; \ - printf(#A "\n"); \ - for (i=0; i _ABS(_(U,max,i))) { max = j; } \ - } \ - _tempi=Piv[i]; Piv[i]=Piv[max]; Piv[max]=_tempi; \ - for (k=i; k _ABS(_(U, max, _i))) { \ + max = _j; \ + } \ + } \ + _tempi = Piv[_i]; \ + Piv[_i] = Piv[max]; \ + Piv[max] = _tempi; \ + for (_k = _i; _k < (n); _k++) { \ + _tempf = _(U, _i, _k); \ + _(U, _i, _k) = _(U, max, _k); \ + _(U, max, _k) = _tempf; \ + } \ + for (_k = 0; _k < _i; _k++) { \ + _tempf = _(L, _i, _k); \ + _(L, _i, _k) = _(L, max, _k); \ + _(L, max, _k) = _tempf; \ + } \ + \ + FLOAT invDiag = 1.0 / _(U, _i, _i); \ + for (_j = _i + 1; _j < (n); _j++) { \ + FLOAT scale = _(U, _j, _i) * invDiag; \ + _(U, _j, _i) = 0.0; \ + for (_k = _i + 1; _k < (n); _k++) { \ + _(U, _j, _k) -= _(U, _i, _k) * scale; \ + } \ + _(L, _j, _i) = scale; \ + } \ + } \ + } /* * Pivots a matrix to a different matrix * R = Pivot(A) given table 'Piv' - * A and R are (n by m) + * A and R are (m by n) */ -#define PIVOT(R,A,Piv,n,m) { \ - int i,j; \ - for (j=0; j=0; j--) { \ - float sum=0.0; \ - for (k=n-1; k>j; k--) { sum += _(U,j,k)*_(X,k,i); } \ - _(X,j,i) = (_(B,j,i) - sum) / _(U,j,j); \ - } \ - } \ -} +#define U_SUB(X, U, B, m, n) \ + { \ + for (int _i = 0; _i < (n); _i++) { \ + for (int _j = m - 1; _j >= 0; _j--) { \ + float sum = 0.0; \ + for (int _k = n - 1; _k > _j; _k--) { \ + sum += _(U, _j, _k) * _(X, _k, _i); \ + } \ + _(X, _j, _i) = (_(B, _j, _i) - sum) / _(U, _j, _j); \ + } \ + } \ + } /* * Inverts a matrix X (n by n) using the method of LU decomposition @@ -187,55 +208,54 @@ PRINT(Ainv,n,n); \ /* * Matrix Multiply R = A * B * R (n by p) - * A (n by m) + * A (m by n) * B (m by p) */ -#define MUL(R,A,B,n,m,p) { \ - int i,j,k; \ - for (i=0; i Date: Sat, 17 Mar 2018 14:21:20 -0400 Subject: Update dcl and test. --- redist/dclapack.h | 20 ++++++++++++++++---- 1 file changed, 16 insertions(+), 4 deletions(-) (limited to 'redist/dclapack.h') 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); \ } \ } \ -- cgit v1.2.3 From f6586a91478b1ca5920c58f90925deaafab81465 Mon Sep 17 00:00:00 2001 From: cnlohr Date: Sat, 17 Mar 2018 14:44:10 -0400 Subject: Update tests... seems to work with both cblas and dclapack... Still haven't checked outputs. --- redist/dclapack.h | 1 - 1 file changed, 1 deletion(-) (limited to 'redist/dclapack.h') diff --git a/redist/dclapack.h b/redist/dclapack.h index cae377b..d4634ac 100644 --- a/redist/dclapack.h +++ b/redist/dclapack.h @@ -269,7 +269,6 @@ 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); \ } \ } \ -- cgit v1.2.3 From 4cb248d46808c07f029c0c595af8f14517925816 Mon Sep 17 00:00:00 2001 From: Justin Berger Date: Sat, 17 Mar 2018 14:00:32 -0600 Subject: Fixed test --- redist/dclapack.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'redist/dclapack.h') diff --git a/redist/dclapack.h b/redist/dclapack.h index d4634ac..e43a4f9 100644 --- a/redist/dclapack.h +++ b/redist/dclapack.h @@ -27,7 +27,7 @@ printf(#A "\n"); \ for (int _i = 0; _i < (m); _i++) { \ for (int _j = 0; _j < (n); _j++) { \ - printf("%4.3f\t", _(A, _i, _j)); \ + printf("%4.3f ", _(A, _i, _j)); \ } \ printf("\n"); \ } \ -- cgit v1.2.3 From 728936f7dce1bc545b5136590a6eb660771266d4 Mon Sep 17 00:00:00 2001 From: Justin Berger Date: Sat, 17 Mar 2018 16:08:49 -0600 Subject: Added transpose timing tests --- redist/dclapack.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'redist/dclapack.h') diff --git a/redist/dclapack.h b/redist/dclapack.h index e43a4f9..7f30187 100644 --- a/redist/dclapack.h +++ b/redist/dclapack.h @@ -24,7 +24,7 @@ */ #define PRINT(A, m, n) \ { \ - printf(#A "\n"); \ + printf(#A " %dx%d\n", m, n); \ for (int _i = 0; _i < (m); _i++) { \ for (int _j = 0; _j < (n); _j++) { \ printf("%4.3f ", _(A, _i, _j)); \ -- cgit v1.2.3