aboutsummaryrefslogtreecommitdiff
path: root/redist/dclapack.h
blob: 8d2aac6eccc97dde29eb10c4f504dcb94a9d7a63 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
#ifndef __DCLAPACK_H__
#define __DCLAPACK_H__

#ifndef FLOAT
#define FLOAT float
#endif

#include<stdio.h>

#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 (m by n)  (with witdth Ac)
 */
#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 (with size n x n, but width Ic)
 */
#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;                                                                                   \
			}                                                                                                          \
		}                                                                                                              \
	}

/*
 * Returns the identity matrix (with size n x n, but width Ic)
 */
#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(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
 *  R = Pivot(A) given table 'Piv'
 *  A and R are (m by n)
 */
#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 m by m (lower triangular)
 *  B is m by n
 */
#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 m by m (upper triangular)
 *  B is m by n
 */
#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
 */

#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); \
PRINT(L,n,n); \
PRINT(U,n,n); \
MUL(L,U,LU,n,n,n);\
PRINT(LU,n,n);\
PRINT(C,n,n); \
PRINT(Ainv,n,n); \
*/

/*
 * Matrix Multiply R = A * B
 *  R (n by p)
 *  A (m by n)
 *  B (m by p)
 */
#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 (m by p)
 *  A (m by n)
 *  B (n by p)
 *  C (m by p)
 */
#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 (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)                                                                      \
	{                                                                                                                  \
		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