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
|
/*
Copyright (c) 2003, The Regents of the University of California, through
Lawrence Berkeley National Laboratory (subject to receipt of any required
approvals from U.S. Dept. of Energy)
All rights reserved.
The source code is distributed under BSD license, see the file License.txt
at the top-level directory.
*/
/*
* -- SuperLU routine (version 3.1) --
* Univ. of California Berkeley, Xerox Palo Alto Research Center,
* and Lawrence Berkeley National Lab.
* August 1, 2008
*
*/
/*! \file
* \brief LU factorization from ZGSTRF (ZGSSVX)
*
* \ingroup Example
*/
#include <getopt.h>
#include "slu_zdefs.h"
void parse_command_line(int argc, char *argv[], int_t *lwork,
double *u, yes_no_t *equil, trans_t *trans);
int main(int argc, char *argv[])
{
char equed[1];
yes_no_t equil;
trans_t trans;
SuperMatrix A, L, U;
SuperMatrix B, X;
NCformat *Astore;
NCformat *Ustore;
SCformat *Lstore;
GlobalLU_t Glu; /* facilitate multiple factorizations with
SamePattern_SameRowPerm */
doublecomplex *a;
int_t *asub, *xa;
int *perm_r; /* row permutations from partial pivoting */
int *perm_c; /* column permutation vector */
int *etree;
void *work = NULL;
int nrhs, ldx;
int_t info, lwork, nnz;
int m, n;
doublecomplex *rhsb, *rhsx, *xact;
double *R, *C;
double *ferr, *berr;
double u, rpg, rcond;
mem_usage_t mem_usage;
superlu_options_t options;
SuperLUStat_t stat;
FILE *fp = stdin;
#if ( DEBUGlevel>=1 )
CHECK_MALLOC("Enter main()");
#endif
/* Defaults */
lwork = 0;
nrhs = 1;
equil = YES;
u = 1.0;
trans = NOTRANS;
/* Set the default input options:
options.Fact = DOFACT;
options.Equil = YES;
options.ColPerm = COLAMD;
options.DiagPivotThresh = 1.0;
options.Trans = NOTRANS;
options.IterRefine = NOREFINE;
options.SymmetricMode = NO;
options.PivotGrowth = NO;
options.ConditionNumber = NO;
options.PrintStat = YES;
*/
set_default_options(&options);
/* Can use command line input to modify the defaults. */
parse_command_line(argc, argv, &lwork, &u, &equil, &trans);
options.Equil = equil;
options.DiagPivotThresh = u;
options.Trans = trans;
/* Add more functionalities that the defaults. */
options.PivotGrowth = YES; /* Compute reciprocal pivot growth */
options.ConditionNumber = YES;/* Compute reciprocal condition number */
options.IterRefine = SLU_DOUBLE; /* Perform double-precision refinement */
if ( lwork > 0 ) {
work = SUPERLU_MALLOC(lwork);
if ( !work ) {
ABORT("ZLINSOLX: cannot allocate work[]");
}
}
/* Read matrix A from a file in Harwell-Boeing format.*/
zreadhb(fp, &m, &n, &nnz, &a, &asub, &xa);
zCreate_CompCol_Matrix(&A, m, n, nnz, a, asub, xa, SLU_NC, SLU_Z, SLU_GE);
Astore = A.Store;
printf("Dimension %dx%d; # nonzeros %d\n", (int)A.nrow, (int)A.ncol, (int)Astore->nnz);
if ( !(rhsb = doublecomplexMalloc(m * nrhs)) ) ABORT("Malloc fails for rhsb[].");
if ( !(rhsx = doublecomplexMalloc(m * nrhs)) ) ABORT("Malloc fails for rhsx[].");
zCreate_Dense_Matrix(&B, m, nrhs, rhsb, m, SLU_DN, SLU_Z, SLU_GE);
zCreate_Dense_Matrix(&X, m, nrhs, rhsx, m, SLU_DN, SLU_Z, SLU_GE);
xact = doublecomplexMalloc(n * nrhs);
ldx = n;
zGenXtrue(n, nrhs, xact, ldx);
zFillRHS(trans, nrhs, xact, ldx, &A, &B);
if ( !(etree = int32Malloc(n)) ) ABORT("Malloc fails for etree[].");
if ( !(perm_r = int32Malloc(m)) ) ABORT("Malloc fails for perm_r[].");
if ( !(perm_c = int32Malloc(n)) ) ABORT("Malloc fails for perm_c[].");
if ( !(R = (double *) SUPERLU_MALLOC(A.nrow * sizeof(double))) )
ABORT("SUPERLU_MALLOC fails for R[].");
if ( !(C = (double *) SUPERLU_MALLOC(A.ncol * sizeof(double))) )
ABORT("SUPERLU_MALLOC fails for C[].");
if ( !(ferr = (double *) SUPERLU_MALLOC(nrhs * sizeof(double))) )
ABORT("SUPERLU_MALLOC fails for ferr[].");
if ( !(berr = (double *) SUPERLU_MALLOC(nrhs * sizeof(double))) )
ABORT("SUPERLU_MALLOC fails for berr[].");
/* Initialize the statistics variables. */
StatInit(&stat);
/* Solve the system and compute the condition number
and error bounds using dgssvx. */
zgssvx(&options, &A, perm_c, perm_r, etree, equed, R, C,
&L, &U, work, lwork, &B, &X, &rpg, &rcond, ferr, berr,
&Glu, &mem_usage, &stat, &info);
printf("zgssvx(): info %lld\n", (long long) info);
if ( info == 0 || info == n+1 ) {
/* This is how you could access the solution matrix. */
doublecomplex *sol = (doublecomplex*) ((DNformat*) X.Store)->nzval;
(void)sol; // suppress unused variable warning
if ( options.PivotGrowth == YES )
printf("Recip. pivot growth = %e\n", rpg);
if ( options.ConditionNumber == YES )
printf("Recip. condition number = %e\n", rcond);
if ( options.IterRefine != NOREFINE ) {
printf("Iterative Refinement:\n");
printf("%8s%8s%16s%16s\n", "rhs", "Steps", "FERR", "BERR");
for (int i = 0; i < nrhs; ++i)
printf("%8d%8d%16e%16e\n", i+1, stat.RefineSteps, ferr[i], berr[i]);
}
Lstore = (SCformat *) L.Store;
Ustore = (NCformat *) U.Store;
printf("No of nonzeros in factor L = %lld\n", (long long) Lstore->nnz);
printf("No of nonzeros in factor U = %lld\n", (long long) Ustore->nnz);
printf("No of nonzeros in L+U = %lld\n", (long long) Lstore->nnz + Ustore->nnz - n);
printf("FILL ratio = %.1f\n", (float)(Lstore->nnz + Ustore->nnz - n)/nnz);
printf("L\\U MB %.3f\ttotal MB needed %.3f\n",
mem_usage.for_lu/1e6, mem_usage.total_needed/1e6);
fflush(stdout);
} else if ( info > 0 && lwork == -1 ) {
printf("** Estimated memory: %lld bytes\n", (long long)info - n);
}
if ( options.PrintStat ) StatPrint(&stat);
StatFree(&stat);
SUPERLU_FREE (rhsb);
SUPERLU_FREE (rhsx);
SUPERLU_FREE (xact);
SUPERLU_FREE (etree);
SUPERLU_FREE (perm_r);
SUPERLU_FREE (perm_c);
SUPERLU_FREE (R);
SUPERLU_FREE (C);
SUPERLU_FREE (ferr);
SUPERLU_FREE (berr);
Destroy_CompCol_Matrix(&A);
Destroy_SuperMatrix_Store(&B);
Destroy_SuperMatrix_Store(&X);
if ( lwork == 0 ) {
Destroy_SuperNode_Matrix(&L);
Destroy_CompCol_Matrix(&U);
} else if ( lwork > 0 ) {
SUPERLU_FREE(work);
}
#if ( DEBUGlevel>=1 )
CHECK_MALLOC("Exit main()");
#endif
return EXIT_SUCCESS;
}
/*!
* \brief Parse command line inputs.
*/
void
parse_command_line(int argc, char *argv[], int_t *lwork,
double *u, yes_no_t *equil, trans_t *trans )
{
int c;
extern char *optarg;
while ( (c = getopt(argc, argv, "hl:w:r:u:f:t:p:e:")) != EOF ) {
switch (c) {
case 'h':
printf("Options:\n");
printf("\t-l <int> - length of work[*] array\n");
printf("\t-u <int> - pivoting threshold\n");
printf("\t-e <0 or 1> - equilibrate or not\n");
printf("\t-t <0 or 1> - solve transposed system or not\n");
exit(1);
break;
case 'l': *lwork = atoi(optarg);
break;
case 'u': *u = atof(optarg);
break;
case 'e': *equil = atoi(optarg);
break;
case 't': *trans = atoi(optarg);
break;
}
}
}
|
