1 /*
2 * -----------------------------------------------------------------
3 * $Revision: 4272 $
4 * $Date: 2014-12-02 11:19:41 -0800 (Tue, 02 Dec 2014) $
5 * -----------------------------------------------------------------
6 * Programmer: Radu Serban @ LLNL
7 * -----------------------------------------------------------------
8 * LLNS Copyright Start
9 * Copyright (c) 2014, Lawrence Livermore National Security
10 * This work was performed under the auspices of the U.S. Department
11 * of Energy by Lawrence Livermore National Laboratory in part under
12 * Contract W-7405-Eng-48 and in part under Contract DE-AC52-07NA27344.
13 * Produced at the Lawrence Livermore National Laboratory.
14 * All rights reserved.
15 * For details, see the LICENSE file.
16 * LLNS Copyright End
17 * -----------------------------------------------------------------
18 * This is the implementation file for operations to be used by a
19 * generic direct linear solver.
20 * -----------------------------------------------------------------
21 */
22
23 #include <stdio.h>
24 #include <stdlib.h>
25
26 #include <sundials/sundials_direct.h>
27 #include <sundials/sundials_math.h>
28
29 #define ZERO RCONST(0.0)
30 #define ONE RCONST(1.0)
31
NewDenseMat(long int M,long int N)32 DlsMat NewDenseMat(long int M, long int N)
33 {
34 DlsMat A;
35 long int j;
36
37 if ( (M <= 0) || (N <= 0) ) return(NULL);
38
39 A = NULL;
40 A = (DlsMat) malloc(sizeof *A);
41 if (A==NULL) return (NULL);
42
43 A->data = (realtype *) malloc(M * N * sizeof(realtype));
44 if (A->data == NULL) {
45 free(A); A = NULL;
46 return(NULL);
47 }
48 A->cols = (realtype **) malloc(N * sizeof(realtype *));
49 if (A->cols == NULL) {
50 free(A->data); A->data = NULL;
51 free(A); A = NULL;
52 return(NULL);
53 }
54
55 for (j=0; j < N; j++) A->cols[j] = A->data + j * M;
56
57 A->M = M;
58 A->N = N;
59 A->ldim = M;
60 A->ldata = M*N;
61
62 A->type = SUNDIALS_DENSE;
63
64 return(A);
65 }
66
newDenseMat(long int m,long int n)67 realtype **newDenseMat(long int m, long int n)
68 {
69 long int j;
70 realtype **a;
71
72 if ( (n <= 0) || (m <= 0) ) return(NULL);
73
74 a = NULL;
75 a = (realtype **) malloc(n * sizeof(realtype *));
76 if (a == NULL) return(NULL);
77
78 a[0] = NULL;
79 a[0] = (realtype *) malloc(m * n * sizeof(realtype));
80 if (a[0] == NULL) {
81 free(a); a = NULL;
82 return(NULL);
83 }
84
85 for (j=1; j < n; j++) a[j] = a[0] + j * m;
86
87 return(a);
88 }
89
90
NewBandMat(long int N,long int mu,long int ml,long int smu)91 DlsMat NewBandMat(long int N, long int mu, long int ml, long int smu)
92 {
93 DlsMat A;
94 long int j, colSize;
95
96 if (N <= 0) return(NULL);
97
98 A = NULL;
99 A = (DlsMat) malloc(sizeof *A);
100 if (A == NULL) return (NULL);
101
102 colSize = smu + ml + 1;
103 A->data = NULL;
104 A->data = (realtype *) malloc(N * colSize * sizeof(realtype));
105 if (A->data == NULL) {
106 free(A); A = NULL;
107 return(NULL);
108 }
109
110 A->cols = NULL;
111 A->cols = (realtype **) malloc(N * sizeof(realtype *));
112 if (A->cols == NULL) {
113 free(A->data);
114 free(A); A = NULL;
115 return(NULL);
116 }
117
118 for (j=0; j < N; j++) A->cols[j] = A->data + j * colSize;
119
120 A->M = N;
121 A->N = N;
122 A->mu = mu;
123 A->ml = ml;
124 A->s_mu = smu;
125 A->ldim = colSize;
126 A->ldata = N * colSize;
127
128 A->type = SUNDIALS_BAND;
129
130 return(A);
131 }
132
newBandMat(long int n,long int smu,long int ml)133 realtype **newBandMat(long int n, long int smu, long int ml)
134 {
135 realtype **a;
136 long int j, colSize;
137
138 if (n <= 0) return(NULL);
139
140 a = NULL;
141 a = (realtype **) malloc(n * sizeof(realtype *));
142 if (a == NULL) return(NULL);
143
144 colSize = smu + ml + 1;
145 a[0] = NULL;
146 a[0] = (realtype *) malloc(n * colSize * sizeof(realtype));
147 if (a[0] == NULL) {
148 free(a); a = NULL;
149 return(NULL);
150 }
151
152 for (j=1; j < n; j++) a[j] = a[0] + j * colSize;
153
154 return(a);
155 }
156
DestroyMat(DlsMat A)157 void DestroyMat(DlsMat A)
158 {
159 free(A->data); A->data = NULL;
160 free(A->cols);
161 free(A); A = NULL;
162 }
163
destroyMat(realtype ** a)164 void destroyMat(realtype **a)
165 {
166 free(a[0]); a[0] = NULL;
167 free(a); a = NULL;
168 }
169
NewIntArray(int N)170 int *NewIntArray(int N)
171 {
172 int *vec;
173
174 if (N <= 0) return(NULL);
175
176 vec = NULL;
177 vec = (int *) malloc(N * sizeof(int));
178
179 return(vec);
180 }
181
newIntArray(int n)182 int *newIntArray(int n)
183 {
184 int *v;
185
186 if (n <= 0) return(NULL);
187
188 v = NULL;
189 v = (int *) malloc(n * sizeof(int));
190
191 return(v);
192 }
193
NewLintArray(long int N)194 long int *NewLintArray(long int N)
195 {
196 long int *vec;
197
198 if (N <= 0) return(NULL);
199
200 vec = NULL;
201 vec = (long int *) malloc(N * sizeof(long int));
202
203 return(vec);
204 }
205
newLintArray(long int n)206 long int *newLintArray(long int n)
207 {
208 long int *v;
209
210 if (n <= 0) return(NULL);
211
212 v = NULL;
213 v = (long int *) malloc(n * sizeof(long int));
214
215 return(v);
216 }
217
NewRealArray(long int N)218 realtype *NewRealArray(long int N)
219 {
220 realtype *vec;
221
222 if (N <= 0) return(NULL);
223
224 vec = NULL;
225 vec = (realtype *) malloc(N * sizeof(realtype));
226
227 return(vec);
228 }
229
newRealArray(long int m)230 realtype *newRealArray(long int m)
231 {
232 realtype *v;
233
234 if (m <= 0) return(NULL);
235
236 v = NULL;
237 v = (realtype *) malloc(m * sizeof(realtype));
238
239 return(v);
240 }
241
DestroyArray(void * V)242 void DestroyArray(void *V)
243 {
244 free(V);
245 V = NULL;
246 }
247
destroyArray(void * v)248 void destroyArray(void *v)
249 {
250 free(v);
251 v = NULL;
252 }
253
254
AddIdentity(DlsMat A)255 void AddIdentity(DlsMat A)
256 {
257 long int i;
258
259 switch (A->type) {
260
261 case SUNDIALS_DENSE:
262 for (i=0; i<A->N; i++) A->cols[i][i] += ONE;
263 break;
264
265 case SUNDIALS_BAND:
266 for (i=0; i<A->M; i++) A->cols[i][A->s_mu] += ONE;
267 break;
268
269 }
270
271 }
272
273
SetToZero(DlsMat A)274 void SetToZero(DlsMat A)
275 {
276 long int i, j, colSize;
277 realtype *col_j;
278
279 switch (A->type) {
280
281 case SUNDIALS_DENSE:
282
283 for (j=0; j<A->N; j++) {
284 col_j = A->cols[j];
285 for (i=0; i<A->M; i++)
286 col_j[i] = ZERO;
287 }
288
289 break;
290
291 case SUNDIALS_BAND:
292
293 colSize = A->mu + A->ml + 1;
294 for (j=0; j<A->M; j++) {
295 col_j = A->cols[j] + A->s_mu - A->mu;
296 for (i=0; i<colSize; i++)
297 col_j[i] = ZERO;
298 }
299
300 break;
301
302 }
303
304 }
305
306
PrintMat(DlsMat A)307 void PrintMat(DlsMat A)
308 {
309 long int i, j, start, finish;
310 realtype **a;
311
312 switch (A->type) {
313
314 case SUNDIALS_DENSE:
315
316 printf("\n");
317 for (i=0; i < A->M; i++) {
318 for (j=0; j < A->N; j++) {
319 #if defined(SUNDIALS_EXTENDED_PRECISION)
320 printf("%12Lg ", DENSE_ELEM(A,i,j));
321 #elif defined(SUNDIALS_DOUBLE_PRECISION)
322 printf("%12g ", DENSE_ELEM(A,i,j));
323 #else
324 printf("%12g ", DENSE_ELEM(A,i,j));
325 #endif
326 }
327 printf("\n");
328 }
329 printf("\n");
330
331 break;
332
333 case SUNDIALS_BAND:
334
335 a = A->cols;
336 printf("\n");
337 for (i=0; i < A->N; i++) {
338 start = SUNMAX(0,i-A->ml);
339 finish = SUNMIN(A->N-1,i+A->mu);
340 for (j=0; j < start; j++) printf("%12s ","");
341 for (j=start; j <= finish; j++) {
342 #if defined(SUNDIALS_EXTENDED_PRECISION)
343 printf("%12Lg ", a[j][i-j+A->s_mu]);
344 #elif defined(SUNDIALS_DOUBLE_PRECISION)
345 printf("%12g ", a[j][i-j+A->s_mu]);
346 #else
347 printf("%12g ", a[j][i-j+A->s_mu]);
348 #endif
349 }
350 printf("\n");
351 }
352 printf("\n");
353
354 break;
355
356 }
357
358 }
359
360
361