diff options
Diffstat (limited to 'redist/dclapack.h')
-rw-r--r-- | redist/dclapack.h | 278 |
1 files changed, 149 insertions, 129 deletions
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<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 "\n"); \ + for (int _i = 0; _i < (m); _i++) { \ + for (int _j = 0; _j < (n); _j++) { \ + printf("%4.3f\t", _(A, _i, _j)); \ + } \ + printf("\n"); \ + } \ + printf("\n"); \ + } /* * 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, n) \ + { \ + for (int _i = 0; _i < (n); _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; \ + } \ + } \ + } /* * R = Transpose(A) - * A is (n by m) - * R is (m by n) + * A is (m by n) + * R is (n by m) */ -#define TRANSP(R,A,n,m) { \ - int i,j; \ - for (i=0; i<n; i++) { \ - for (j=0; j<m; j++) { \ - _(R,j,i) = _(A,i,j); \ - } \ - } \ -} +#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(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); \ - \ - 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); \ + \ + 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 * 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<n; j++) { \ - for (i=0; i<m; i++) { \ - _(R,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(X,L,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(X,U,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 @@ -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<n; i++) { \ - for (j=0; j<p; j++) { \ - _(R,i,j) = 0.0f; \ - for (k=0; k<m; k++) { \ - _(R,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 R = A * B + C - * R (n by p) - * A (n by m) - * B (m by p) - * C (n by p) + * R (m by p) + * A (m by n) + * B (n by p) + * C (m by p) */ -#define MULADD(R,A,B,C,n,m,p) { \ - int i,j,k; \ - for (i=0; i<n; i++) { \ - for (j=0; j<p; j++) { \ - _(R,i,j) = _(C,i,j); \ - for (k=0; k<m; k++) { \ - _(R,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 R = alpha * A * B + beta * C - * R (n by p) + * R (m by p) * A (m by n) * B (n by p) * C (m by p) */ #define GMULADD(R, A, B, C, alpha, beta, m, n, p) \ { \ - int _i, _j, _k; \ float sum; \ - for (_i = 0; _i < m; _i++) { \ - for (_j = 0; _j < p; _j++) { \ + for (int _i = 0; _i < m; _i++) { \ + for (int _j = 0; _j < p; _j++) { \ sum = 0.0f; \ - for (_k = 0; _k < n; _k++) { \ + for (int _k = 0; _k < n; _k++) { \ sum += _(A, _i, _k) * _(B, _k, _j); \ } \ _(R, _i, _j) = alpha * sum + beta * _(C, _i, _j); \ |