1 | /* ----------------------------------------------------------------- |
---|
2 | * Programmer: Radu Serban @ LLNL |
---|
3 | * ----------------------------------------------------------------- |
---|
4 | * LLNS Copyright Start |
---|
5 | * Copyright (c) 2014, Lawrence Livermore National Security |
---|
6 | * This work was performed under the auspices of the U.S. Department |
---|
7 | * of Energy by Lawrence Livermore National Laboratory in part under |
---|
8 | * Contract W-7405-Eng-48 and in part under Contract DE-AC52-07NA27344. |
---|
9 | * Produced at the Lawrence Livermore National Laboratory. |
---|
10 | * All rights reserved. |
---|
11 | * For details, see the LICENSE file. |
---|
12 | * LLNS Copyright End |
---|
13 | * ----------------------------------------------------------------- |
---|
14 | * This header file contains definitions and declarations for use by |
---|
15 | * generic direct linear solvers for Ax = b. It defines types for |
---|
16 | * dense and banded matrices and corresponding accessor macros. |
---|
17 | * -----------------------------------------------------------------*/ |
---|
18 | |
---|
19 | #ifndef _SUNDIALS_DIRECT_H |
---|
20 | #define _SUNDIALS_DIRECT_H |
---|
21 | |
---|
22 | #include <stdio.h> |
---|
23 | #include <sundials/sundials_types.h> |
---|
24 | |
---|
25 | #ifdef __cplusplus /* wrapper to enable C++ usage */ |
---|
26 | extern "C" { |
---|
27 | #endif |
---|
28 | |
---|
29 | /* |
---|
30 | * ================================================================= |
---|
31 | * C O N S T A N T S |
---|
32 | * ================================================================= |
---|
33 | */ |
---|
34 | |
---|
35 | /* |
---|
36 | * SUNDIALS_DENSE: dense matrix |
---|
37 | * SUNDIALS_BAND: banded matrix |
---|
38 | */ |
---|
39 | |
---|
40 | #define SUNDIALS_DENSE 1 |
---|
41 | #define SUNDIALS_BAND 2 |
---|
42 | |
---|
43 | /* |
---|
44 | * ================================================================== |
---|
45 | * Type definitions |
---|
46 | * ================================================================== |
---|
47 | */ |
---|
48 | |
---|
49 | /* |
---|
50 | * ----------------------------------------------------------------- |
---|
51 | * Type : DlsMat |
---|
52 | * ----------------------------------------------------------------- |
---|
53 | * The type DlsMat is defined to be a pointer to a structure |
---|
54 | * with various sizes, a data field, and an array of pointers to |
---|
55 | * the columns which defines a dense or band matrix for use in |
---|
56 | * direct linear solvers. The M and N fields indicates the number |
---|
57 | * of rows and columns, respectively. The data field is a one |
---|
58 | * dimensional array used for component storage. The cols field |
---|
59 | * stores the pointers in data for the beginning of each column. |
---|
60 | * ----------------------------------------------------------------- |
---|
61 | * For DENSE matrices, the relevant fields in DlsMat are: |
---|
62 | * type = SUNDIALS_DENSE |
---|
63 | * M - number of rows |
---|
64 | * N - number of columns |
---|
65 | * ldim - leading dimension (ldim >= M) |
---|
66 | * data - pointer to a contiguous block of realtype variables |
---|
67 | * ldata - length of the data array =ldim*N |
---|
68 | * cols - array of pointers. cols[j] points to the first element |
---|
69 | * of the j-th column of the matrix in the array data. |
---|
70 | * |
---|
71 | * The elements of a dense matrix are stored columnwise (i.e. columns |
---|
72 | * are stored one on top of the other in memory). |
---|
73 | * If A is of type DlsMat, then the (i,j)th element of A (with |
---|
74 | * 0 <= i < M and 0 <= j < N) is given by (A->data)[j*n+i]. |
---|
75 | * |
---|
76 | * The DENSE_COL and DENSE_ELEM macros below allow a user to access |
---|
77 | * efficiently individual matrix elements without writing out explicit |
---|
78 | * data structure references and without knowing too much about the |
---|
79 | * underlying element storage. The only storage assumption needed is |
---|
80 | * that elements are stored columnwise and that a pointer to the |
---|
81 | * jth column of elements can be obtained via the DENSE_COL macro. |
---|
82 | * ----------------------------------------------------------------- |
---|
83 | * For BAND matrices, the relevant fields in DlsMat are: |
---|
84 | * type = SUNDIALS_BAND |
---|
85 | * M - number of rows |
---|
86 | * N - number of columns |
---|
87 | * mu - upper bandwidth, 0 <= mu <= min(M,N) |
---|
88 | * ml - lower bandwidth, 0 <= ml <= min(M,N) |
---|
89 | * s_mu - storage upper bandwidth, mu <= s_mu <= N-1. |
---|
90 | * The dgbtrf routine writes the LU factors into the storage |
---|
91 | * for A. The upper triangular factor U, however, may have |
---|
92 | * an upper bandwidth as big as MIN(N-1,mu+ml) because of |
---|
93 | * partial pivoting. The s_mu field holds the upper |
---|
94 | * bandwidth allocated for A. |
---|
95 | * ldim - leading dimension (ldim >= s_mu) |
---|
96 | * data - pointer to a contiguous block of realtype variables |
---|
97 | * ldata - length of the data array =ldim*(s_mu+ml+1) |
---|
98 | * cols - array of pointers. cols[j] points to the first element |
---|
99 | * of the j-th column of the matrix in the array data. |
---|
100 | * |
---|
101 | * The BAND_COL, BAND_COL_ELEM, and BAND_ELEM macros below allow a |
---|
102 | * user to access individual matrix elements without writing out |
---|
103 | * explicit data structure references and without knowing too much |
---|
104 | * about the underlying element storage. The only storage assumption |
---|
105 | * needed is that elements are stored columnwise and that a pointer |
---|
106 | * into the jth column of elements can be obtained via the BAND_COL |
---|
107 | * macro. The BAND_COL_ELEM macro selects an element from a column |
---|
108 | * which has already been isolated via BAND_COL. The macro |
---|
109 | * BAND_COL_ELEM allows the user to avoid the translation |
---|
110 | * from the matrix location (i,j) to the index in the array returned |
---|
111 | * by BAND_COL at which the (i,j)th element is stored. |
---|
112 | * ----------------------------------------------------------------- |
---|
113 | */ |
---|
114 | |
---|
115 | typedef struct _DlsMat { |
---|
116 | int type; |
---|
117 | sunindextype M; |
---|
118 | sunindextype N; |
---|
119 | sunindextype ldim; |
---|
120 | sunindextype mu; |
---|
121 | sunindextype ml; |
---|
122 | sunindextype s_mu; |
---|
123 | realtype *data; |
---|
124 | sunindextype ldata; |
---|
125 | realtype **cols; |
---|
126 | } *DlsMat; |
---|
127 | |
---|
128 | /* |
---|
129 | * ================================================================== |
---|
130 | * Data accessor macros |
---|
131 | * ================================================================== |
---|
132 | */ |
---|
133 | |
---|
134 | /* |
---|
135 | * ----------------------------------------------------------------- |
---|
136 | * DENSE_COL and DENSE_ELEM |
---|
137 | * ----------------------------------------------------------------- |
---|
138 | * |
---|
139 | * DENSE_COL(A,j) references the jth column of the M-by-N dense |
---|
140 | * matrix A, 0 <= j < N. The type of the expression DENSE_COL(A,j) |
---|
141 | * is (realtype *). After the assignment col_j = DENSE_COL(A,j), |
---|
142 | * col_j may be treated as an array indexed from 0 to M-1. |
---|
143 | * The (i,j)-th element of A is thus referenced by col_j[i]. |
---|
144 | * |
---|
145 | * DENSE_ELEM(A,i,j) references the (i,j)th element of the dense |
---|
146 | * M-by-N matrix A, 0 <= i < M ; 0 <= j < N. |
---|
147 | * |
---|
148 | * ----------------------------------------------------------------- |
---|
149 | */ |
---|
150 | |
---|
151 | #define DENSE_COL(A,j) ((A->cols)[j]) |
---|
152 | #define DENSE_ELEM(A,i,j) ((A->cols)[j][i]) |
---|
153 | |
---|
154 | /* |
---|
155 | * ----------------------------------------------------------------- |
---|
156 | * BAND_COL, BAND_COL_ELEM, and BAND_ELEM |
---|
157 | * ----------------------------------------------------------------- |
---|
158 | * |
---|
159 | * BAND_COL(A,j) references the diagonal element of the jth column |
---|
160 | * of the N by N band matrix A, 0 <= j <= N-1. The type of the |
---|
161 | * expression BAND_COL(A,j) is realtype *. The pointer returned by |
---|
162 | * the call BAND_COL(A,j) can be treated as an array which is |
---|
163 | * indexed from -(A->mu) to (A->ml). |
---|
164 | * |
---|
165 | * BAND_COL_ELEM references the (i,j)th entry of the band matrix A |
---|
166 | * when used in conjunction with BAND_COL. The index (i,j) should |
---|
167 | * satisfy j-(A->mu) <= i <= j+(A->ml). |
---|
168 | * |
---|
169 | * BAND_ELEM(A,i,j) references the (i,j)th element of the M-by-N |
---|
170 | * band matrix A, where 0 <= i,j <= N-1. The location (i,j) should |
---|
171 | * further satisfy j-(A->mu) <= i <= j+(A->ml). |
---|
172 | * |
---|
173 | * ----------------------------------------------------------------- |
---|
174 | */ |
---|
175 | |
---|
176 | #define BAND_COL(A,j) (((A->cols)[j])+(A->s_mu)) |
---|
177 | #define BAND_COL_ELEM(col_j,i,j) (col_j[(i)-(j)]) |
---|
178 | #define BAND_ELEM(A,i,j) ((A->cols)[j][(i)-(j)+(A->s_mu)]) |
---|
179 | |
---|
180 | /* |
---|
181 | * ================================================================== |
---|
182 | * Exported function prototypes (functions working on dlsMat) |
---|
183 | * ================================================================== |
---|
184 | */ |
---|
185 | |
---|
186 | /* |
---|
187 | * ----------------------------------------------------------------- |
---|
188 | * Function: NewDenseMat |
---|
189 | * ----------------------------------------------------------------- |
---|
190 | * NewDenseMat allocates memory for an M-by-N dense matrix and |
---|
191 | * returns the storage allocated (type DlsMat). NewDenseMat |
---|
192 | * returns NULL if the request for matrix storage cannot be |
---|
193 | * satisfied. See the above documentation for the type DlsMat |
---|
194 | * for matrix storage details. |
---|
195 | * ----------------------------------------------------------------- |
---|
196 | */ |
---|
197 | |
---|
198 | SUNDIALS_EXPORT DlsMat NewDenseMat(sunindextype M, sunindextype N); |
---|
199 | |
---|
200 | /* |
---|
201 | * ----------------------------------------------------------------- |
---|
202 | * Function: NewBandMat |
---|
203 | * ----------------------------------------------------------------- |
---|
204 | * NewBandMat allocates memory for an M-by-N band matrix |
---|
205 | * with upper bandwidth mu, lower bandwidth ml, and storage upper |
---|
206 | * bandwidth smu. Pass smu as follows depending on whether A will |
---|
207 | * be LU factored: |
---|
208 | * |
---|
209 | * (1) Pass smu = mu if A will not be factored. |
---|
210 | * |
---|
211 | * (2) Pass smu = MIN(N-1,mu+ml) if A will be factored. |
---|
212 | * |
---|
213 | * NewBandMat returns the storage allocated (type DlsMat) or |
---|
214 | * NULL if the request for matrix storage cannot be satisfied. |
---|
215 | * See the documentation for the type DlsMat for matrix storage |
---|
216 | * details. |
---|
217 | * ----------------------------------------------------------------- |
---|
218 | */ |
---|
219 | |
---|
220 | SUNDIALS_EXPORT DlsMat NewBandMat(sunindextype N, sunindextype mu, sunindextype ml, sunindextype smu); |
---|
221 | |
---|
222 | /* |
---|
223 | * ----------------------------------------------------------------- |
---|
224 | * Functions: DestroyMat |
---|
225 | * ----------------------------------------------------------------- |
---|
226 | * DestroyMat frees the memory allocated by NewDenseMat or NewBandMat |
---|
227 | * ----------------------------------------------------------------- |
---|
228 | */ |
---|
229 | |
---|
230 | SUNDIALS_EXPORT void DestroyMat(DlsMat A); |
---|
231 | |
---|
232 | /* |
---|
233 | * ----------------------------------------------------------------- |
---|
234 | * Function: NewIntArray |
---|
235 | * ----------------------------------------------------------------- |
---|
236 | * NewIntArray allocates memory an array of N int's and returns |
---|
237 | * the pointer to the memory it allocates. If the request for |
---|
238 | * memory storage cannot be satisfied, it returns NULL. |
---|
239 | * ----------------------------------------------------------------- |
---|
240 | */ |
---|
241 | |
---|
242 | SUNDIALS_EXPORT int *NewIntArray(int N); |
---|
243 | |
---|
244 | /* |
---|
245 | * ----------------------------------------------------------------- |
---|
246 | * Function: NewIndexArray |
---|
247 | * ----------------------------------------------------------------- |
---|
248 | * NewIndexArray allocates memory an array of N sunindextype's and |
---|
249 | * returns the pointer to the memory it allocates. If the request |
---|
250 | * for memory storage cannot be satisfied, it returns NULL. |
---|
251 | * ----------------------------------------------------------------- |
---|
252 | */ |
---|
253 | |
---|
254 | SUNDIALS_EXPORT sunindextype *NewIndexArray(sunindextype N); |
---|
255 | |
---|
256 | /* |
---|
257 | * ----------------------------------------------------------------- |
---|
258 | * Function: NewRealArray |
---|
259 | * ----------------------------------------------------------------- |
---|
260 | * NewRealArray allocates memory an array of N realtype and returns |
---|
261 | * the pointer to the memory it allocates. If the request for |
---|
262 | * memory storage cannot be satisfied, it returns NULL. |
---|
263 | * ----------------------------------------------------------------- |
---|
264 | */ |
---|
265 | |
---|
266 | SUNDIALS_EXPORT realtype *NewRealArray(sunindextype N); |
---|
267 | |
---|
268 | /* |
---|
269 | * ----------------------------------------------------------------- |
---|
270 | * Function: DestroyArray |
---|
271 | * ----------------------------------------------------------------- |
---|
272 | * DestroyArray frees memory allocated by NewIntArray, NewIndexArray, |
---|
273 | * or NewRealArray. |
---|
274 | * ----------------------------------------------------------------- |
---|
275 | */ |
---|
276 | |
---|
277 | SUNDIALS_EXPORT void DestroyArray(void *p); |
---|
278 | |
---|
279 | /* |
---|
280 | * ----------------------------------------------------------------- |
---|
281 | * Function : AddIdentity |
---|
282 | * ----------------------------------------------------------------- |
---|
283 | * AddIdentity adds 1.0 to the main diagonal (A_ii, i=0,1,...,N-1) of |
---|
284 | * the M-by-N matrix A (M>= N) and stores the result back in A. |
---|
285 | * AddIdentity is typically used with square matrices. |
---|
286 | * AddIdentity does not check for M >= N and therefore a segmentation |
---|
287 | * fault will occur if M < N! |
---|
288 | * ----------------------------------------------------------------- |
---|
289 | */ |
---|
290 | |
---|
291 | SUNDIALS_EXPORT void AddIdentity(DlsMat A); |
---|
292 | |
---|
293 | /* |
---|
294 | * ----------------------------------------------------------------- |
---|
295 | * Function : SetToZero |
---|
296 | * ----------------------------------------------------------------- |
---|
297 | * SetToZero sets all the elements of the M-by-N matrix A to 0.0. |
---|
298 | * ----------------------------------------------------------------- |
---|
299 | */ |
---|
300 | |
---|
301 | SUNDIALS_EXPORT void SetToZero(DlsMat A); |
---|
302 | |
---|
303 | /* |
---|
304 | * ----------------------------------------------------------------- |
---|
305 | * Functions: PrintMat |
---|
306 | * ----------------------------------------------------------------- |
---|
307 | * This function prints the M-by-N (dense or band) matrix A to |
---|
308 | * outfile as it would normally appear on paper. |
---|
309 | * It is intended as debugging tools with small values of M and N. |
---|
310 | * The elements are printed using the %g/%lg/%Lg option. |
---|
311 | * A blank line is printed before and after the matrix. |
---|
312 | * ----------------------------------------------------------------- |
---|
313 | */ |
---|
314 | |
---|
315 | SUNDIALS_EXPORT void PrintMat(DlsMat A, FILE *outfile); |
---|
316 | |
---|
317 | |
---|
318 | /* |
---|
319 | * ================================================================== |
---|
320 | * Exported function prototypes (functions working on realtype**) |
---|
321 | * ================================================================== |
---|
322 | */ |
---|
323 | |
---|
324 | SUNDIALS_EXPORT realtype **newDenseMat(sunindextype m, sunindextype n); |
---|
325 | SUNDIALS_EXPORT realtype **newBandMat(sunindextype n, sunindextype smu, sunindextype ml); |
---|
326 | SUNDIALS_EXPORT void destroyMat(realtype **a); |
---|
327 | SUNDIALS_EXPORT int *newIntArray(int n); |
---|
328 | SUNDIALS_EXPORT sunindextype *newIndexArray(sunindextype n); |
---|
329 | SUNDIALS_EXPORT realtype *newRealArray(sunindextype m); |
---|
330 | SUNDIALS_EXPORT void destroyArray(void *v); |
---|
331 | |
---|
332 | |
---|
333 | #ifdef __cplusplus |
---|
334 | } |
---|
335 | #endif |
---|
336 | |
---|
337 | #endif |
---|