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
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
|
/******************************************************************
* *
* File : dense.h *
* Programmers : Scott D. Cohen and Alan C. Hindmarsh @ LLNL *
* Last Modified : 1 September 1994 *
*----------------------------------------------------------------*
* This is the header file for a generic DENSE linear solver *
* package. There are two sets of dense solver routines listed in *
* this file: one set uses type DenseMat defined below and the *
* other set uses the type real ** for dense matrix arguments. *
* The two sets of dense solver routines make it easy to work *
* with two types of dense matrices: *
* *
* (1) The DenseMat type is intended for use with large dense *
* matrices whose elements/columns may be stored in *
* non-contiguous memory locations or even distributed into *
* different processor memories. This type may be modified to *
* include such distribution information. If this is done, *
* then all the routines that use DenseMat must be modified *
* to reflect the new data structure. *
* *
* (2) The set of routines that use real ** (and NOT the DenseMat *
* type) is intended for use with small matrices which can *
* easily be allocated within a contiguous block of memory *
* on a single processor. *
* *
* Routines that work with the type DenseMat begin with "Dense". *
* The DenseAllocMat function allocates a dense matrix for use in *
* the other DenseMat routines listed in this file. Matrix *
* storage details are given in the documentation for the type *
* DenseMat. The DenseAllocPiv function allocates memory for *
* pivot information. The storage allocated by DenseAllocMat and *
* DenseAllocPiv is deallocated by the routines DenseFreeMat and *
* DenseFreePiv, respectively. The DenseFactor and DenseBacksolve *
* routines perform the actual solution of a dense linear system. *
* Note that the DenseBacksolve routine has a parameter b of type *
* N_Vector. The current implementation makes use of a machine *
* environment-specific macro (N_VDATA) which may not exist for *
* other implementations of the type N_Vector. Thus, the *
* implementation of DenseBacksolve may need to change if the *
* type N_Vector is changed. *
* *
* Routines that work with real ** begin with "den" (except for *
* the factor and solve routines which are called gefa and gesl, *
* respectively). The underlying matrix storage is described in *
* the documentation for denalloc. *
* *
******************************************************************/
#ifndef _dense_h
#define _dense_h
#include "llnltyps.h"
#include "vector.h"
/******************************************************************
* *
* Type: DenseMat *
*----------------------------------------------------------------*
* The type DenseMat is defined to be a pointer to a structure *
* with a size and a data field. The size field indicates the *
* number of columns (== number of rows) of a dense matrix, while *
* the data field is a two dimensional array used for component *
* storage. The elements of a dense matrix are stored columnwise *
* (i.e columns are stored one on top of the other in memory). If *
* A is of type DenseMat, then the (i,j)th element of A (with *
* 0 <= i,j <= size-1) is given by the expression (A->data)[j][i] *
* or by the expression (A->data)[0][j*n+i]. The macros below *
* allow a user to access efficiently individual matrix *
* elements without writing out explicit data structure *
* references and without knowing too much about the underlying *
* element storage. The only storage assumption needed is that *
* elements are stored columnwise and that a pointer to the jth *
* column of elements can be obtained via the DENSE_COL macro. *
* Users should use these macros whenever possible. *
* *
******************************************************************/
typedef struct {
integer size;
real **data;
} *DenseMat;
/* DenseMat accessor macros */
/******************************************************************
* *
* Macro : DENSE_ELEM *
* Usage : DENSE_ELEM(A,i,j) = a_ij; OR *
* a_ij = DENSE_ELEM(A,i,j); *
*----------------------------------------------------------------*
* DENSE_ELEM(A,i,j) references the (i,j)th element of the N by N *
* DenseMat A, 0 <= i,j <= N-1. *
* *
******************************************************************/
#define DENSE_ELEM(A,i,j) ((A->data)[j][i])
/******************************************************************
* *
* Macro : DENSE_COL *
* Usage : col_j = DENSE_COL(A,j); *
*----------------------------------------------------------------*
* DENSE_COL(A,j) references the jth column of the N by N *
* DenseMat A, 0 <= j <= N-1. The type of the expression *
* DENSE_COL(A,j) is real *. After the assignment in the usage *
* above, col_j may be treated as an array indexed from 0 to N-1. *
* The (i,j)th element of A is referenced by col_j[i]. *
* *
******************************************************************/
#define DENSE_COL(A,j) ((A->data)[j])
/* Functions that use the DenseMat representation for a dense matrix */
/******************************************************************
* *
* Function : DenseAllocMat *
* Usage : A = DenseAllocMat(N); *
* if (A == NULL) ... memory request failed *
*----------------------------------------------------------------*
* DenseAllocMat allocates memory for an N by N dense matrix and *
* returns the storage allocated (type DenseMat). DenseAllocMat *
* returns NULL if the request for matrix storage cannot be *
* satisfied. See the above documentation for the type DenseMat *
* for matrix storage details. *
* *
******************************************************************/
DenseMat DenseAllocMat(integer N);
/******************************************************************
* *
* Function : DenseAllocPiv *
* Usage : p = DenseAllocPiv(N); *
* if (p == NULL) ... memory request failed *
*----------------------------------------------------------------*
* DenseAllocPiv allocates memory for pivot information to be *
* filled in by the DenseFactor routine during the factorization *
* of an N by N dense matrix. The underlying type for pivot *
* information is an array of N integers and this routine returns *
* the pointer to the memory it allocates. If the request for *
* pivot storage cannot be satisfied, DenseAllocPiv returns NULL. *
* *
******************************************************************/
integer *DenseAllocPiv(integer N);
/******************************************************************
* *
* Function : DenseFactor *
* Usage : ier = DenseFactor(A, p); *
* if (ier != 0) ... A is singular *
*----------------------------------------------------------------*
* DenseFactor performs the LU factorization of the N by N dense *
* matrix A. This is done using standard Gaussian elimination *
* with partial pivoting. *
* *
* A successful LU factorization leaves the matrix A and the *
* pivot array p with the following information: *
* *
* (1) p[k] contains the row number of the pivot element chosen *
* at the beginning of elimination step k, k=0, 1, ..., N-1. *
* *
* (2) If the unique LU factorization of A is given by PA = LU, *
* where P is a permutation matrix, L is a lower triangular *
* matrix with all 1's on the diagonal, and U is an upper *
* triangular matrix, then the upper triangular part of A *
* (including its diagonal) contains U and the strictly lower *
* triangular part of A contains the multipliers, I-L. *
* *
* DenseFactor returns 0 if successful. Otherwise it encountered *
* a zero diagonal element during the factorization. In this case *
* it returns the column index (numbered from one) at which *
* it encountered the zero. *
* *
******************************************************************/
integer DenseFactor(DenseMat A, integer *p);
/******************************************************************
* *
* Function : DenseBacksolve *
* Usage : DenseBacksolve(A, p, b); *
*----------------------------------------------------------------*
* DenseBacksolve solves the N-dimensional system A x = b using *
* the LU factorization in A and the pivot information in p *
* computed in DenseFactor. The solution x is returned in b. This *
* routine cannot fail if the corresponding call to DenseFactor *
* did not fail. *
* *
******************************************************************/
void DenseBacksolve(DenseMat A, integer *p, N_Vector b);
/******************************************************************
* *
* Function : DenseZero *
* Usage : DenseZero(A); *
*----------------------------------------------------------------*
* DenseZero sets all the elements of the N by N matrix A to 0.0. *
* *
******************************************************************/
void DenseZero(DenseMat A);
/******************************************************************
* *
* Function : DenseCopy *
* Usage : DenseCopy(A, B); *
*----------------------------------------------------------------*
* DenseCopy copies the contents of the N by N matrix A into the *
* N by N matrix B. *
* *
******************************************************************/
void DenseCopy(DenseMat A, DenseMat B);
/******************************************************************
* *
* Function: DenseScale *
* Usage : DenseScale(c, A); *
*----------------------------------------------------------------*
* DenseScale scales the elements of the N by N matrix A by the *
* constant c and stores the result back in A. *
* *
******************************************************************/
void DenseScale(real c, DenseMat A);
/******************************************************************
* *
* Function : DenseAddI *
* Usage : DenseAddI(A); *
*----------------------------------------------------------------*
* DenseAddI adds the identity matrix to A and stores the result *
* back in A. *
* *
******************************************************************/
void DenseAddI(DenseMat A);
/******************************************************************
* *
* Function : DenseFreeMat *
* Usage : DenseFreeMat(A); *
*----------------------------------------------------------------*
* DenseFreeMat frees the memory allocated by DenseAllocMat for *
* the N by N matrix A. *
* *
******************************************************************/
void DenseFreeMat(DenseMat A);
/******************************************************************
* *
* Function : DenseFreePiv *
* Usage : DenseFreePiv(p); *
*----------------------------------------------------------------*
* DenseFreePiv frees the memory allocated by DenseAllocPiv for *
* the pivot information array p. *
* *
******************************************************************/
void DenseFreePiv(integer *p);
/******************************************************************
* *
* Function : DensePrint *
* Usage : DensePrint(A); *
*----------------------------------------------------------------*
* This routine prints the N by N dense matrix A to standard *
* output as it would normally appear on paper. It is intended *
* as a debugging tool with small values of N. The elements are *
* printed using the %g option. A blank line is printed before *
* and after the matrix. *
* *
******************************************************************/
void DensePrint(DenseMat A);
/* Functions that use the real ** representation for a dense matrix */
/******************************************************************
* *
* Function : denalloc *
* Usage : real **a; *
* a = denalloc(n); *
* if (a == NULL) ... memory request failed *
*----------------------------------------------------------------*
* denalloc(n) allocates storage for an n by n dense matrix. It *
* returns a pointer to the newly allocated storage if *
* successful. If the memory request cannot be satisfied, then *
* denalloc returns NULL. The underlying type of the dense matrix *
* returned is real **. If we allocate a dense matrix real **a by *
* a = denalloc(n), then a[j][i] references the (i,j)th element *
* of the matrix a, 0 <= i,j <= n-1, and a[j] is a pointer to the *
* first element in the jth column of a. The location a[0] *
* contains a pointer to n^2 contiguous locations which contain *
* the elements of a. *
* *
******************************************************************/
real **denalloc(integer n);
/******************************************************************
* *
* Function : denallocpiv *
* Usage : integer *pivot; *
* pivot = denallocpiv(n); *
* if (pivot == NULL) ... memory request failed *
*----------------------------------------------------------------*
* denallocpiv(n) allocates an array of n integers. It returns a *
* pointer to the first element in the array if successful. It *
* returns NULL if the memory request could not be satisfied. *
* *
******************************************************************/
integer *denallocpiv(integer n);
/******************************************************************
* *
* Function : gefa *
* Usage : integer ier; *
* ier = gefa(a,n,p); *
* if (ier > 0) ... zero element encountered during *
* the factorization *
*----------------------------------------------------------------*
* gefa(a,n,p) factors the n by n dense matrix a. It overwrites *
* the elements of a with its LU factors and keeps track of the *
* pivot rows chosen in the pivot array p. *
* *
* A successful LU factorization leaves the matrix a and the *
* pivot array p with the following information: *
* *
* (1) p[k] contains the row number of the pivot element chosen *
* at the beginning of elimination step k, k=0, 1, ..., n-1. *
* *
* (2) If the unique LU factorization of a is given by Pa = LU, *
* where P is a permutation matrix, L is a lower triangular *
* matrix with all 1's on the diagonal, and U is an upper *
* triangular matrix, then the upper triangular part of a *
* (including its diagonal) contains U and the strictly lower *
* triangular part of a contains the multipliers, I-L. *
* *
* gefa returns 0 if successful. Otherwise it encountered a zero *
* diagonal element during the factorization. In this case it *
* returns the column index (numbered from one) at which it *
* encountered the zero. *
* *
******************************************************************/
integer gefa(real **a, integer n, integer *p);
/******************************************************************
* *
* Function : gesl *
* Usage : real *b; *
* ier = gefa(a,n,p); *
* if (ier == 0) gesl(a,n,p,b); *
*----------------------------------------------------------------*
* gesl(a,n,p,b) solves the n by n linear system ax = b. It *
* assumes that a has been LU factored and the pivot array p has *
* been set by a successful call to gefa(a,n,p). The solution x *
* is written into the b array. *
* *
******************************************************************/
void gesl(real **a, integer n, integer *p, real *b);
/******************************************************************
* *
* Function : denzero *
* Usage : denzero(a,n); *
*----------------------------------------------------------------*
* denzero(a,n) sets all the elements of the n by n dense matrix *
* a to be 0.0. *
* *
******************************************************************/
void denzero(real **a, integer n);
/******************************************************************
* *
* Function : dencopy *
* Usage : dencopy(a,b,n); *
*----------------------------------------------------------------*
* dencopy(a,b,n) copies the n by n dense matrix a into the *
* n by n dense matrix b. *
* *
******************************************************************/
void dencopy(real **a, real **b, integer n);
/******************************************************************
* *
* Function : denscale *
* Usage : denscale(c,a,n); *
*----------------------------------------------------------------*
* denscale(c,a,n) scales every element in the n by n dense *
* matrix a by c. *
* *
******************************************************************/
void denscale(real c, real **a, integer n);
/******************************************************************
* *
* Function : denaddI *
* Usage : denaddI(a,n); *
*----------------------------------------------------------------*
* denaddI(a,n) increments the n by n dense matrix a by the *
* identity matrix. *
* *
******************************************************************/
void denaddI(real **a, integer n);
/******************************************************************
* *
* Function : denfreepiv *
* Usage : denfreepiv(p); *
*----------------------------------------------------------------*
* denfreepiv(p) frees the pivot array p allocated by *
* denallocpiv. *
* *
******************************************************************/
void denfreepiv(integer *p);
/******************************************************************
* *
* Function : denfree *
* Usage : denfree(a); *
*----------------------------------------------------------------*
* denfree(a) frees the dense matrix a allocated by denalloc. *
* *
******************************************************************/
void denfree(real **a);
/******************************************************************
* *
* Function : denprint *
* Usage : denprint(a,n); *
*----------------------------------------------------------------*
* denprint(a,n) prints the n by n dense matrix a to standard *
* output as it would normally appear on paper. It is intended as *
* a debugging tool with small values of n. The elements are *
* printed using the %g option. A blank line is printed before *
* and after the matrix. *
* *
******************************************************************/
void denprint(real **a, integer n);
#endif
|
