From 08eec7a6bd5edc0fd1cae7f98d7788ea3905a467 Mon Sep 17 00:00:00 2001 From: cnlohr 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') 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.2.3