diff options
-rw-r--r-- | redist/dclapack.h | 78 | ||||
-rw-r--r-- | redist/dclhelpers.c | 63 | ||||
-rw-r--r-- | redist/dclhelpers.h | 61 |
3 files changed, 171 insertions, 31 deletions
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<n; i++) { \ for (j=0; j<m; j++) { \ - printf("%4.3f\t", A[i][j]); \ + printf("%4.3f\t", _(A,i,j)); \ } \ printf("\n"); \ } \ @@ -34,9 +37,9 @@ #define IDENTITY(I,n) { \ int i,j; \ for (i=0; i<n; i++) { \ - for (j=0; j<i; j++) { I[i][j]=0.0f; } \ - I[i][i] = 1.0f; \ - for (j=i+1; j<n; j++) { I[i][j]=0.0f; } \ + for (j=0; j<i; j++) { _(I,i,j)=0.0f; } \ + _(I,i,i) = 1.0f; \ + for (j=i+1; j<n; j++) { _(I,i,j)=0.0f; } \ } \ } @@ -49,7 +52,7 @@ int i,j; \ for (i=0; i<n; i++) { \ for (j=0; j<m; j++) { \ - B[j][i] = A[i][j]; \ + _(B,j,i) = _(A,i,j); \ } \ } \ } @@ -62,7 +65,7 @@ for (i=0; i<n; i++) { Piv[i]=i; } \ for (i=0; i<n; i++) { \ for (j=0; j<n; j++) { \ - U[i][j] = A[i][j]; \ + _(U,i,j) = _(A,i,j); \ } \ } \ IDENTITY(L,n); \ @@ -71,22 +74,22 @@ \ int max=i; \ for (j=i+1; j<n; j++) { \ - if (_ABS(U[j][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<n; k++) { \ - _tempf=U[i][k]; U[i][k]=U[max][k]; U[max][k]=_tempf; \ + _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; \ + _tempf=_(L,i,k); _(L,i,k)=_(L,max,k); _(L,max,k)=_tempf; \ } \ \ - FLOAT invDiag = 1.0 / U[i][i]; \ + 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; \ + 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; \ } \ } \ } @@ -100,7 +103,7 @@ int i,j; \ for (j=0; j<n; j++) { \ for (i=0; i<m; i++) { \ - B[j][i] = A[Piv[j]][i]; \ + _(B,j,i) = _(A,Piv[j],i); \ } \ } \ } @@ -115,8 +118,8 @@ for (i=0; i<m; i++) { \ for (j=0; j<n; j++) { \ float sum=0.0; \ - for (k=0; k<j; k++) { sum += L[j][k]*X[k][i]; } \ - X[j][i] = (B[j][i] - sum) / L[j][j]; \ + for (k=0; k<j; k++) { sum += _(L,j,k)*_(X,k,i); } \ + _(X,j,i) = (_(B,j,i) - sum) / _(L,j,j); \ } \ } \ } @@ -131,8 +134,8 @@ for (i=0; i<m; i++) { \ for (j=n-1; 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]; \ + 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<n; i++) { \ for (j=0; j<p; j++) { \ - C[i][j] = 0.0f; \ + _(C,i,j) = 0.0f; \ for (k=0; k<m; k++) { \ - C[i][j] += A[i][k] * B[k][j]; \ + _(C,i,j) += _(A,i,k) * _(B,k,j); \ } \ } \ } \ @@ -193,9 +209,9 @@ PRINT(Ainv,n,n); \ int i,j,k; \ for (i=0; i<n; i++) { \ for (j=0; j<p; j++) { \ - D[i][j] = C[i][j]; \ + _(D,i,j) = _(C,i,j); \ for (k=0; k<m; k++) { \ - D[i][j] += A[i][k] * B[k][j]; \ + _(D,i,j) += _(A,i,k) * _(B,k,j); \ } \ } \ } \ @@ -215,9 +231,9 @@ PRINT(Ainv,n,n); \ for (j=0; j<p; j++) { \ sum = 0.0f; \ for (k=0; k<m; k++) { \ - sum += A[i][k] * B[k][j]; \ + sum += _(A,i,k) * _(B,k,j); \ } \ - D[i][j] = alpha * sum + beta * C[i][j]; \ + _(D,i,j) = alpha * sum + beta * _(C,i,j); \ } \ } \ } 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 <stdio.h> +#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); +} + + diff --git a/redist/dclhelpers.h b/redist/dclhelpers.h new file mode 100644 index 0000000..fbaa7ec --- /dev/null +++ b/redist/dclhelpers.h @@ -0,0 +1,61 @@ +#ifndef _DCL_HELPERS_H +#define _DCL_HELPERS_H + +#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 n, int m ); + +/* Returns the identity matrix */ +void dclIdentity( DCL_FLOAT * A, int n ); + +/* B = Transpose(A) + A is (n by m) + B is (m by n) */ +void dclTransp( const DCL_FLOAT * A, DCL_FLOAT * B, int n, int m ); + +/* Calculate L,U of a matrix A with pivot table; the pivot table is output. */ +void dclLU( const DCL_FLOAT * A, DCL_FLOAT * L, DCL_FLOAT * U, int * Piv, int n ); + +/* Pivots a matrix to a different matrix + B = Pivot(A) given table 'Piv' + A and B are (n by m) */ +void dclPivot( const DCL_FLOAT * A, DCL_FLOAT * B, 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( const DCL_FLOAT * L, DCL_FLOAT * X, const DCL_FLOAT * B, 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( const DCL_FLOAT * U, DCL_FLOAT * X, const DCL_FLOAT * B, int n, int m ); + +/* Inverts a matrix X (n by n) using the method of LU decomposition */ +void dclInv( const DCL_FLOAT * A, DCL_FLOAT * Ainv, int n ); + +/* Matrix Multiply C = A * B + A (n by m) + B (m by p) + C (n by p) */ +void dclMul( const DCL_FLOAT * A, const DCL_FLOAT * B, DCL_FLOAT * C, int n, int m, int p ); + +/* Matrix Multiply D = A * B + C + A (n by m) + B (m by p) + C (n by p) + D (n by p) */ +void dclMulAdd( const DCL_FLOAT * A, const DCL_FLOAT * B, const DCL_FLOAT * C, DCL_FLOAT * D, int n, int m, int p ); + +/* Matrix Multiply D = alpha * A * B + beta * C + A (n by m) + B (m by p) + C (n by p) + D (n by 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 ); + +#endif + |