aboutsummaryrefslogtreecommitdiff
path: root/redist
diff options
context:
space:
mode:
authorJustin Berger <j.david.berger@gmail.com>2018-03-17 10:33:14 -0600
committerJustin Berger <j.david.berger@gmail.com>2018-03-17 10:33:14 -0600
commit5d4f4068ae265271d0840ca1bddd7dfcddd948e9 (patch)
tree18a173a34c5f3a04983797b996fcb02b6379df2b /redist
parent0b9e66ad2ff686a4dcf8a6838f33edb203a1bff5 (diff)
downloadlibsurvive-5d4f4068ae265271d0840ca1bddd7dfcddd948e9.tar.gz
libsurvive-5d4f4068ae265271d0840ca1bddd7dfcddd948e9.tar.bz2
Made naming follow conventions
Diffstat (limited to 'redist')
-rw-r--r--redist/dclapack.h278
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); \