diff options
Diffstat (limited to 'redist/dclapack.h')
| -rw-r--r-- | redist/dclapack.h | 385 |
1 files changed, 221 insertions, 164 deletions
diff --git a/redist/dclapack.h b/redist/dclapack.h index 4e209d3..0950815 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,146 +9,207 @@ #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 _(AM, O, P) AM[O * AM##c + P] +#else +#define _(AM, O, P) AM[O][P] +#endif + /* - * Prints a matrix A (n by m) + * 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<n; i++) { \ - for (j=0; j<m; j++) { \ - printf("%4.3f\t", A[i][j]); \ - } \ - printf("\n"); \ - } \ - printf("\n"); \ -} +#define PRINT(A, m, 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)); \ + } \ + printf("\n"); \ + } \ + printf("\n"); \ + } /* - * Returns the identity matrix + * Returns the identity matrix (with size n x n, but width Ic) */ -#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; } \ - } \ -} +#define IDENTITY(I, m, n) \ + { \ + for (int _i = 0; _i < (m); _i++) { \ + for (int _j = 0; _j < _i; _j++) { \ + _(I, _i, _j) = 0.0f; \ + } \ + _(I, _i, _i) = 1.0f; \ + for (int _j = _i + 1; _j < (n); _j++) { \ + _(I, _i, _j) = 0.0f; \ + } \ + } \ + } /* - * B = Transpose(A) - * A is (n by m) - * B is (m by n) + * Returns the identity matrix (with size n x n, but width Ic) */ -#define TRANSP(A,B,n,m) { \ - int i,j; \ - for (i=0; i<n; i++) { \ - for (j=0; j<m; j++) { \ - B[j][i] = A[i][j]; \ - } \ - } \ -} +#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) + * R is (n by m) + */ +#define TRANSP(R, A, m, n) \ + { \ + for (int _i = 0; _i < (m); _i++) { \ + for (int _j = 0; _j < (n); _j++) { \ + _(R, _j, _i) = _(A, _i, _j); \ + } \ + } \ + } /* * Calculate L,U of a matrix A with pivot table */ -#define LU(A,L,U,Piv,n) { \ - int i,j,k,_tempi; float _tempf; \ - 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]; \ - } \ - } \ - IDENTITY(L,n); \ - \ - for (i=0; i<n-1; i++) { \ - \ - int max=i; \ - for (j=i+1; j<n; 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; \ - } \ - 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; \ - } \ - } \ -} +#define LU(L, U, A, Piv, n) \ + { \ + int _i, _j, _k, _tempi; \ + float _tempf; \ + 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); \ + } \ + } \ + IDENTITY(L, n, n); \ + \ + for (_i = 0; _i < (n)-1; _i++) { \ + \ + int max = _i; \ + for (_j = _i + 1; _j < (n); _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; \ + } \ + 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 - * B = Pivot(A) given table 'Piv' - * A and B are (n by m) + * R = Pivot(A) given table 'Piv' + * A and R are (m by n) */ -#define PIVOT(A,B,Piv,n,m) { \ - int i,j; \ - for (j=0; j<n; j++) { \ - for (i=0; i<m; i++) { \ - B[j][i] = A[Piv[j]][i]; \ - } \ - } \ -} +#define PIVOT(R, A, Piv, m, n) \ + { \ + for (int _i = 0; _i < (m); _i++) { \ + for (int _j = 0; _j < (n); _j++) { \ + _(R, _i, _j) = _(A, Piv[_i], _j); \ + } \ + } \ + } /* * Solve LX=B for matrix X and B - * L is n by n (lower triangular) - * B is n by m + * L is m by m (lower triangular) + * B is m by n */ -#define L_SUB(L,X,B,n,m) { \ - int i,j,k; \ - 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]; \ - } \ - } \ -} +#define L_SUB(X, L, B, m, n) \ + { \ + for (int _i = 0; _i < (n); _i++) { \ + for (int _j = 0; _j < (m); _j++) { \ + float sum = 0.0; \ + for (int _k = 0; _k < _j; _k++) { \ + sum += _(L, _j, _k) * _(X, _k, _i); \ + } \ + _(X, _j, _i) = (_(B, _j, _i) - sum) / _(L, _j, _j); \ + } \ + } \ + } /* * Solve UX=B for matrix X and B - * U is n by n (upper triangular) - * B is n by m + * U is m by m (upper triangular) + * B is m by n */ -#define U_SUB(U,X,B,n,m) { \ - int i,j,k; \ - 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]; \ - } \ - } \ -} +#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 */ -#define INV(A,Ainv,n) { \ - FLOAT Ipiv[ORDER][ORDER]; \ - FLOAT L[ORDER][ORDER]; \ - FLOAT U[ORDER][ORDER]; \ - FLOAT I[ORDER][ORDER]; \ - FLOAT C[ORDER][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); \ -} + +#ifdef DYNAMIC_INDEX +#define INV_SETUP(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]; \ + FLOAT L[ORDER][ORDER]; \ + FLOAT U[ORDER][ORDER]; \ + FLOAT I[ORDER][ORDER]; \ + FLOAT C[ORDER][ORDER]; +#endif + +#define INV(Ainv, A, n, ORDER) \ + { \ + INV_SETUP(ORDER) \ + int Piv[ORDER]; \ + IDENTITY(I, 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); \ + } /* PRINT(A,n,n); \ @@ -165,61 +222,61 @@ PRINT(Ainv,n,n); \ */ /* - * Matrix Multiply C = A * B - * A (n by m) + * Matrix Multiply R = A * B + * R (n by p) + * A (m by n) * B (m by p) - * C (n by p) */ -#define MUL(A,B,C,n,m,p) { \ - int i,j,k; \ - for (i=0; i<n; i++) { \ - for (j=0; j<p; j++) { \ - C[i][j] = 0.0f; \ - for (k=0; k<m; k++) { \ - C[i][j] += A[i][k] * B[k][j]; \ - } \ - } \ - } \ -} +#define MUL(R, A, B, m, n, p) \ + { \ + for (int _i = 0; _i < (m); _i++) { \ + for (int _j = 0; _j < p; _j++) { \ + _(R, _i, _j) = 0.0f; \ + for (int _k = 0; _k < (n); _k++) { \ + _(R, _i, _j) += _(A, _i, _k) * _(B, _k, _j); \ + } \ + } \ + } \ + } /* - * Matrix Multiply D = A * B + C - * A (n by m) - * B (m by p) - * C (n by p) - * D (n by p) + * Matrix Multiply R = A * B + C + * R (m by p) + * A (m by n) + * B (n by p) + * C (m by p) */ -#define MULADD(A,B,C,D,n,m,p) { \ - int i,j,k; \ - for (i=0; i<n; i++) { \ - for (j=0; j<p; j++) { \ - D[i][j] = C[i][j]; \ - for (k=0; k<m; k++) { \ - D[i][j] += A[i][k] * B[k][j]; \ - } \ - } \ - } \ -} +#define MULADD(R, A, B, C, m, n, p) \ + { \ + for (int _i = 0; _i < (m); _i++) { \ + for (int _j = 0; _j < p; _j++) { \ + _(R, _i, _j) = _(C, _i, _j); \ + for (int _k = 0; _k < (n); _k++) { \ + _(R, _i, _j) += _(A, _i, _k) * _(B, _k, _j); \ + } \ + } \ + } \ + } /* - * Matrix Multiply D = alpha * A * B + beta * C - * A (n by m) - * B (m by p) - * C (n by p) - * D (n by p) + * Matrix Multiply R = alpha * A * B + beta * C + * R (m by p) + * A (m by n) + * B (n by p) + * C (m by p) */ -#define GMULADD(A,B,C,D,alpha,beta,n,m,p) { \ - int i,j,k; \ - float sum; \ - for (i=0; i<n; i++) { \ - for (j=0; j<p; j++) { \ - sum = 0.0f; \ - for (k=0; k<m; k++) { \ - sum += A[i][k] * B[k][j]; \ - } \ - D[i][j] = alpha * sum + beta * C[i][j]; \ - } \ - } \ -} +#define GMULADD(R, A, B, C, alpha, beta, m, n, p) \ + { \ + float sum; \ + for (int _i = 0; _i < m; _i++) { \ + for (int _j = 0; _j < p; _j++) { \ + sum = 0.0f; \ + for (int _k = 0; _k < n; _k++) { \ + sum += _(A, _i, _k) * _(B, _k, _j); \ + } \ + _(R, _i, _j) = alpha * sum + beta * _(C, _i, _j); \ + } \ + } \ + } #endif |