3
* \brief Gives the approximate solutions of linear equations A*X=B or A'*X=B
6
* -- SuperLU routine (version 4.0) --
7
* Lawrence Berkeley National Laboratory.
11
#include "slu_cdefs.h"
19
* CGSISX gives the approximate solutions of linear equations A*X=B or A'*X=B,
20
* using the ILU factorization from cgsitrf(). An estimation of
21
* the condition number is provided. It performs the following steps:
23
* 1. If A is stored column-wise (A->Stype = SLU_NC):
25
* 1.1. If options->Equil = YES or options->RowPerm = LargeDiag, scaling
26
* factors are computed to equilibrate the system:
27
* options->Trans = NOTRANS:
28
* diag(R)*A*diag(C) *inv(diag(C))*X = diag(R)*B
29
* options->Trans = TRANS:
30
* (diag(R)*A*diag(C))**T *inv(diag(R))*X = diag(C)*B
31
* options->Trans = CONJ:
32
* (diag(R)*A*diag(C))**H *inv(diag(R))*X = diag(C)*B
33
* Whether or not the system will be equilibrated depends on the
34
* scaling of the matrix A, but if equilibration is used, A is
35
* overwritten by diag(R)*A*diag(C) and B by diag(R)*B
36
* (if options->Trans=NOTRANS) or diag(C)*B (if options->Trans
39
* 1.2. Permute columns of A, forming A*Pc, where Pc is a permutation
40
* matrix that usually preserves sparsity.
41
* For more details of this step, see sp_preorder.c.
43
* 1.3. If options->Fact != FACTORED, the LU decomposition is used to
44
* factor the matrix A (after equilibration if options->Equil = YES)
45
* as Pr*A*Pc = L*U, with Pr determined by partial pivoting.
47
* 1.4. Compute the reciprocal pivot growth factor.
49
* 1.5. If some U(i,i) = 0, so that U is exactly singular, then the
50
* routine fills a small number on the diagonal entry, that is
51
* U(i,i) = ||A(:,i)||_oo * options->ILU_FillTol ** (1 - i / n),
52
* and info will be increased by 1. The factored form of A is used
53
* to estimate the condition number of the preconditioner. If the
54
* reciprocal of the condition number is less than machine precision,
55
* info = A->ncol+1 is returned as a warning, but the routine still
56
* goes on to solve for X.
58
* 1.6. The system of equations is solved for X using the factored form
61
* 1.7. options->IterRefine is not used
63
* 1.8. If equilibration was used, the matrix X is premultiplied by
64
* diag(C) (if options->Trans = NOTRANS) or diag(R)
65
* (if options->Trans = TRANS or CONJ) so that it solves the
66
* original system before equilibration.
68
* 1.9. options for ILU only
69
* 1) If options->RowPerm = LargeDiag, MC64 is used to scale and
70
* permute the matrix to an I-matrix, that is Pr*Dr*A*Dc has
71
* entries of modulus 1 on the diagonal and off-diagonal entries
72
* of modulus at most 1. If MC64 fails, dgsequ() is used to
73
* equilibrate the system.
74
* 2) options->ILU_DropTol = tau is the threshold for dropping.
75
* For L, it is used directly (for the whole row in a supernode);
76
* For U, ||A(:,i)||_oo * tau is used as the threshold
77
* for the i-th column.
78
* If a secondary dropping rule is required, tau will
79
* also be used to compute the second threshold.
80
* 3) options->ILU_FillFactor = gamma, used as the initial guess
82
* If a secondary dropping rule is required, it will also
83
* be used as an upper bound of the memory.
84
* 4) options->ILU_DropRule specifies the dropping rule.
87
* DROP_BASIC: Basic dropping rule, supernodal based ILU.
88
* DROP_PROWS: Supernodal based ILUTP, p = gamma * nnz(A) / n.
89
* DROP_COLUMN: Variation of ILUTP, for j-th column,
90
* p = gamma * nnz(A(:,j)).
91
* DROP_AREA; Variation of ILUTP, for j-th column, use
92
* nnz(F(:,1:j)) / nnz(A(:,1:j)) to control the
94
* DROP_DYNAMIC: Modify the threshold tau during the
96
* If nnz(L(:,1:j)) / nnz(A(:,1:j)) < gamma
97
* tau_L(j) := MIN(1, tau_L(j-1) * 2);
99
* tau_L(j) := MIN(1, tau_L(j-1) * 2);
100
* tau_U(j) uses the similar rule.
101
* NOTE: the thresholds used by L and U are
103
* DROP_INTERP: Compute the second dropping threshold by
104
* interpolation instead of sorting (default).
105
* In this case, the actual fill ratio is not
106
* guaranteed smaller than gamma.
107
* DROP_PROWS, DROP_COLUMN and DROP_AREA are mutually exclusive.
108
* ( The default option is DROP_BASIC | DROP_AREA. )
109
* 5) options->ILU_Norm is the criterion of computing the average
110
* value of a row in L.
111
* options->ILU_Norm average(x[1:n])
112
* ================= ===============
113
* ONE_NORM ||x||_1 / n
114
* TWO_NORM ||x||_2 / sqrt(n)
115
* INF_NORM max{|x[i]|}
116
* 6) options->ILU_MILU specifies the type of MILU's variation.
117
* = SILU (default): do not perform MILU;
118
* = SMILU_1 (not recommended):
119
* U(i,i) := U(i,i) + sum(dropped entries);
121
* U(i,i) := U(i,i) + SGN(U(i,i)) * sum(dropped entries);
123
* U(i,i) := U(i,i) + SGN(U(i,i)) * sum(|dropped entries|);
124
* NOTE: Even SMILU_1 does not preserve the column sum because of
126
* 7) options->ILU_FillTol is used as the perturbation when
127
* encountering zero pivots. If some U(i,i) = 0, so that U is
128
* exactly singular, then
129
* U(i,i) := ||A(:,i)|| * options->ILU_FillTol ** (1 - i / n).
131
* 2. If A is stored row-wise (A->Stype = SLU_NR), apply the above algorithm
132
* to the transpose of A:
134
* 2.1. If options->Equil = YES or options->RowPerm = LargeDiag, scaling
135
* factors are computed to equilibrate the system:
136
* options->Trans = NOTRANS:
137
* diag(R)*A*diag(C) *inv(diag(C))*X = diag(R)*B
138
* options->Trans = TRANS:
139
* (diag(R)*A*diag(C))**T *inv(diag(R))*X = diag(C)*B
140
* options->Trans = CONJ:
141
* (diag(R)*A*diag(C))**H *inv(diag(R))*X = diag(C)*B
142
* Whether or not the system will be equilibrated depends on the
143
* scaling of the matrix A, but if equilibration is used, A' is
144
* overwritten by diag(R)*A'*diag(C) and B by diag(R)*B
145
* (if trans='N') or diag(C)*B (if trans = 'T' or 'C').
147
* 2.2. Permute columns of transpose(A) (rows of A),
148
* forming transpose(A)*Pc, where Pc is a permutation matrix that
149
* usually preserves sparsity.
150
* For more details of this step, see sp_preorder.c.
152
* 2.3. If options->Fact != FACTORED, the LU decomposition is used to
153
* factor the transpose(A) (after equilibration if
154
* options->Fact = YES) as Pr*transpose(A)*Pc = L*U with the
155
* permutation Pr determined by partial pivoting.
157
* 2.4. Compute the reciprocal pivot growth factor.
159
* 2.5. If some U(i,i) = 0, so that U is exactly singular, then the
160
* routine fills a small number on the diagonal entry, that is
161
* U(i,i) = ||A(:,i)||_oo * options->ILU_FillTol ** (1 - i / n).
162
* And info will be increased by 1. The factored form of A is used
163
* to estimate the condition number of the preconditioner. If the
164
* reciprocal of the condition number is less than machine precision,
165
* info = A->ncol+1 is returned as a warning, but the routine still
166
* goes on to solve for X.
168
* 2.6. The system of equations is solved for X using the factored form
171
* 2.7. If options->IterRefine is not used.
173
* 2.8. If equilibration was used, the matrix X is premultiplied by
174
* diag(C) (if options->Trans = NOTRANS) or diag(R)
175
* (if options->Trans = TRANS or CONJ) so that it solves the
176
* original system before equilibration.
178
* See supermatrix.h for the definition of 'SuperMatrix' structure.
183
* options (input) superlu_options_t*
184
* The structure defines the input parameters to control
185
* how the LU decomposition will be performed and how the
186
* system will be solved.
188
* A (input/output) SuperMatrix*
189
* Matrix A in A*X=B, of dimension (A->nrow, A->ncol). The number
190
* of the linear equations is A->nrow. Currently, the type of A can be:
191
* Stype = SLU_NC or SLU_NR, Dtype = SLU_C, Mtype = SLU_GE.
192
* In the future, more general A may be handled.
194
* On entry, If options->Fact = FACTORED and equed is not 'N',
195
* then A must have been equilibrated by the scaling factors in
197
* On exit, A is not modified if options->Equil = NO, or if
198
* options->Equil = YES but equed = 'N' on exit.
199
* Otherwise, if options->Equil = YES and equed is not 'N',
200
* A is scaled as follows:
201
* If A->Stype = SLU_NC:
202
* equed = 'R': A := diag(R) * A
203
* equed = 'C': A := A * diag(C)
204
* equed = 'B': A := diag(R) * A * diag(C).
205
* If A->Stype = SLU_NR:
206
* equed = 'R': transpose(A) := diag(R) * transpose(A)
207
* equed = 'C': transpose(A) := transpose(A) * diag(C)
208
* equed = 'B': transpose(A) := diag(R) * transpose(A) * diag(C).
210
* perm_c (input/output) int*
211
* If A->Stype = SLU_NC, Column permutation vector of size A->ncol,
212
* which defines the permutation matrix Pc; perm_c[i] = j means
213
* column i of A is in position j in A*Pc.
214
* On exit, perm_c may be overwritten by the product of the input
215
* perm_c and a permutation that postorders the elimination tree
216
* of Pc'*A'*A*Pc; perm_c is not changed if the elimination tree
217
* is already in postorder.
219
* If A->Stype = SLU_NR, column permutation vector of size A->nrow,
220
* which describes permutation of columns of transpose(A)
221
* (rows of A) as described above.
223
* perm_r (input/output) int*
224
* If A->Stype = SLU_NC, row permutation vector of size A->nrow,
225
* which defines the permutation matrix Pr, and is determined
226
* by partial pivoting. perm_r[i] = j means row i of A is in
227
* position j in Pr*A.
229
* If A->Stype = SLU_NR, permutation vector of size A->ncol, which
230
* determines permutation of rows of transpose(A)
231
* (columns of A) as described above.
233
* If options->Fact = SamePattern_SameRowPerm, the pivoting routine
234
* will try to use the input perm_r, unless a certain threshold
235
* criterion is violated. In that case, perm_r is overwritten by a
236
* new permutation determined by partial pivoting or diagonal
237
* threshold pivoting.
238
* Otherwise, perm_r is output argument.
240
* etree (input/output) int*, dimension (A->ncol)
241
* Elimination tree of Pc'*A'*A*Pc.
242
* If options->Fact != FACTORED and options->Fact != DOFACT,
243
* etree is an input argument, otherwise it is an output argument.
244
* Note: etree is a vector of parent pointers for a forest whose
245
* vertices are the integers 0 to A->ncol-1; etree[root]==A->ncol.
247
* equed (input/output) char*
248
* Specifies the form of equilibration that was done.
249
* = 'N': No equilibration.
250
* = 'R': Row equilibration, i.e., A was premultiplied by diag(R).
251
* = 'C': Column equilibration, i.e., A was postmultiplied by diag(C).
252
* = 'B': Both row and column equilibration, i.e., A was replaced
253
* by diag(R)*A*diag(C).
254
* If options->Fact = FACTORED, equed is an input argument,
255
* otherwise it is an output argument.
257
* R (input/output) float*, dimension (A->nrow)
258
* The row scale factors for A or transpose(A).
259
* If equed = 'R' or 'B', A (if A->Stype = SLU_NC) or transpose(A)
260
* (if A->Stype = SLU_NR) is multiplied on the left by diag(R).
261
* If equed = 'N' or 'C', R is not accessed.
262
* If options->Fact = FACTORED, R is an input argument,
263
* otherwise, R is output.
264
* If options->zFact = FACTORED and equed = 'R' or 'B', each element
265
* of R must be positive.
267
* C (input/output) float*, dimension (A->ncol)
268
* The column scale factors for A or transpose(A).
269
* If equed = 'C' or 'B', A (if A->Stype = SLU_NC) or transpose(A)
270
* (if A->Stype = SLU_NR) is multiplied on the right by diag(C).
271
* If equed = 'N' or 'R', C is not accessed.
272
* If options->Fact = FACTORED, C is an input argument,
273
* otherwise, C is output.
274
* If options->Fact = FACTORED and equed = 'C' or 'B', each element
275
* of C must be positive.
277
* L (output) SuperMatrix*
278
* The factor L from the factorization
279
* Pr*A*Pc=L*U (if A->Stype SLU_= NC) or
280
* Pr*transpose(A)*Pc=L*U (if A->Stype = SLU_NR).
281
* Uses compressed row subscripts storage for supernodes, i.e.,
282
* L has types: Stype = SLU_SC, Dtype = SLU_C, Mtype = SLU_TRLU.
284
* U (output) SuperMatrix*
285
* The factor U from the factorization
286
* Pr*A*Pc=L*U (if A->Stype = SLU_NC) or
287
* Pr*transpose(A)*Pc=L*U (if A->Stype = SLU_NR).
288
* Uses column-wise storage scheme, i.e., U has types:
289
* Stype = SLU_NC, Dtype = SLU_C, Mtype = SLU_TRU.
291
* work (workspace/output) void*, size (lwork) (in bytes)
292
* User supplied workspace, should be large enough
293
* to hold data structures for factors L and U.
294
* On exit, if fact is not 'F', L and U point to this array.
297
* Specifies the size of work array in bytes.
298
* = 0: allocate space internally by system malloc;
299
* > 0: use user-supplied work array of length lwork in bytes,
300
* returns error if space runs out.
301
* = -1: the routine guesses the amount of space needed without
302
* performing the factorization, and returns it in
303
* mem_usage->total_needed; no other side effects.
305
* See argument 'mem_usage' for memory usage statistics.
307
* B (input/output) SuperMatrix*
308
* B has types: Stype = SLU_DN, Dtype = SLU_C, Mtype = SLU_GE.
309
* On entry, the right hand side matrix.
310
* If B->ncol = 0, only LU decomposition is performed, the triangular
313
* if equed = 'N', B is not modified; otherwise
314
* if A->Stype = SLU_NC:
315
* if options->Trans = NOTRANS and equed = 'R' or 'B',
316
* B is overwritten by diag(R)*B;
317
* if options->Trans = TRANS or CONJ and equed = 'C' of 'B',
318
* B is overwritten by diag(C)*B;
319
* if A->Stype = SLU_NR:
320
* if options->Trans = NOTRANS and equed = 'C' or 'B',
321
* B is overwritten by diag(C)*B;
322
* if options->Trans = TRANS or CONJ and equed = 'R' of 'B',
323
* B is overwritten by diag(R)*B.
325
* X (output) SuperMatrix*
326
* X has types: Stype = SLU_DN, Dtype = SLU_C, Mtype = SLU_GE.
327
* If info = 0 or info = A->ncol+1, X contains the solution matrix
328
* to the original system of equations. Note that A and B are modified
329
* on exit if equed is not 'N', and the solution to the equilibrated
330
* system is inv(diag(C))*X if options->Trans = NOTRANS and
331
* equed = 'C' or 'B', or inv(diag(R))*X if options->Trans = 'T' or 'C'
332
* and equed = 'R' or 'B'.
334
* recip_pivot_growth (output) float*
335
* The reciprocal pivot growth factor max_j( norm(A_j)/norm(U_j) ).
336
* The infinity norm is used. If recip_pivot_growth is much less
337
* than 1, the stability of the LU factorization could be poor.
339
* rcond (output) float*
340
* The estimate of the reciprocal condition number of the matrix A
341
* after equilibration (if done). If rcond is less than the machine
342
* precision (in particular, if rcond = 0), the matrix is singular
343
* to working precision. This condition is indicated by a return
346
* mem_usage (output) mem_usage_t*
347
* Record the memory usage statistics, consisting of following fields:
349
* The amount of space used in bytes for L\U data structures.
350
* - total_needed (float)
351
* The amount of space needed in bytes to perform factorization.
353
* The number of memory expansions during the LU factorization.
355
* stat (output) SuperLUStat_t*
356
* Record the statistics on runtime and floating-point operation count.
357
* See slu_util.h for the definition of 'SuperLUStat_t'.
360
* = 0: successful exit
361
* < 0: if info = -i, the i-th argument had an illegal value
362
* > 0: if info = i, and i is
363
* <= A->ncol: number of zero pivots. They are replaced by small
364
* entries due to options->ILU_FillTol.
365
* = A->ncol+1: U is nonsingular, but RCOND is less than machine
366
* precision, meaning that the matrix is singular to
367
* working precision. Nevertheless, the solution and
368
* error bounds are computed because there are a number
369
* of situations where the computed solution can be more
370
* accurate than the value of RCOND would suggest.
371
* > A->ncol+1: number of bytes allocated when memory allocation
372
* failure occurred, plus A->ncol.
377
cgsisx(superlu_options_t *options, SuperMatrix *A, int *perm_c, int *perm_r,
378
int *etree, char *equed, float *R, float *C,
379
SuperMatrix *L, SuperMatrix *U, void *work, int lwork,
380
SuperMatrix *B, SuperMatrix *X,
381
float *recip_pivot_growth, float *rcond,
382
mem_usage_t *mem_usage, SuperLUStat_t *stat, int *info)
385
DNformat *Bstore, *Xstore;
386
complex *Bmat, *Xmat;
388
SuperMatrix *AA;/* A in SLU_NC format used by the factorization routine.*/
389
SuperMatrix AC; /* Matrix postmultiplied by Pc */
390
int colequ, equil, nofact, notran, rowequ, permc_spec, mc64;
394
float amax, anorm, bignum, smlnum, colcnd, rowcnd, rcmax, rcmin;
395
int relax, panel_size;
396
float diag_pivot_thresh;
397
double t0; /* temporary time */
402
/* External functions */
403
extern float clangs(char *, SuperMatrix *);
407
Bmat = Bstore->nzval;
408
Xmat = Xstore->nzval;
414
nofact = (options->Fact != FACTORED);
415
equil = (options->Equil == YES);
416
notran = (options->Trans == NOTRANS);
417
mc64 = (options->RowPerm == LargeDiag);
419
*(unsigned char *)equed = 'N';
423
rowequ = lsame_(equed, "R") || lsame_(equed, "B");
424
colequ = lsame_(equed, "C") || lsame_(equed, "B");
425
smlnum = slamch_("Safe minimum");
426
bignum = 1. / smlnum;
429
/* Test the input parameters */
430
if (!nofact && options->Fact != DOFACT && options->Fact != SamePattern &&
431
options->Fact != SamePattern_SameRowPerm &&
432
!notran && options->Trans != TRANS && options->Trans != CONJ &&
433
!equil && options->Equil != NO)
435
else if ( A->nrow != A->ncol || A->nrow < 0 ||
436
(A->Stype != SLU_NC && A->Stype != SLU_NR) ||
437
A->Dtype != SLU_C || A->Mtype != SLU_GE )
439
else if (options->Fact == FACTORED &&
440
!(rowequ || colequ || lsame_(equed, "N")))
446
for (j = 0; j < A->nrow; ++j) {
447
rcmin = SUPERLU_MIN(rcmin, R[j]);
448
rcmax = SUPERLU_MAX(rcmax, R[j]);
450
if (rcmin <= 0.) *info = -7;
451
else if ( A->nrow > 0)
452
rowcnd = SUPERLU_MAX(rcmin,smlnum) / SUPERLU_MIN(rcmax,bignum);
455
if (colequ && *info == 0) {
458
for (j = 0; j < A->nrow; ++j) {
459
rcmin = SUPERLU_MIN(rcmin, C[j]);
460
rcmax = SUPERLU_MAX(rcmax, C[j]);
462
if (rcmin <= 0.) *info = -8;
463
else if (A->nrow > 0)
464
colcnd = SUPERLU_MAX(rcmin,smlnum) / SUPERLU_MIN(rcmax,bignum);
468
if ( lwork < -1 ) *info = -12;
469
else if ( B->ncol < 0 || Bstore->lda < SUPERLU_MAX(0, A->nrow) ||
470
B->Stype != SLU_DN || B->Dtype != SLU_C ||
473
else if ( X->ncol < 0 || Xstore->lda < SUPERLU_MAX(0, A->nrow) ||
474
(B->ncol != 0 && B->ncol != X->ncol) ||
475
X->Stype != SLU_DN ||
476
X->Dtype != SLU_C || X->Mtype != SLU_GE )
482
xerbla_("cgsisx", &i);
486
/* Initialization for factor parameters */
487
panel_size = sp_ienv(1);
489
diag_pivot_thresh = options->DiagPivotThresh;
493
/* Convert A to SLU_NC format when necessary. */
494
if ( A->Stype == SLU_NR ) {
495
NRformat *Astore = A->Store;
496
AA = (SuperMatrix *) SUPERLU_MALLOC( sizeof(SuperMatrix) );
497
cCreate_CompCol_Matrix(AA, A->ncol, A->nrow, Astore->nnz,
498
Astore->nzval, Astore->colind, Astore->rowptr,
499
SLU_NC, A->Dtype, A->Mtype);
500
if ( notran ) { /* Reverse the transpose argument. */
507
} else { /* A->Stype == SLU_NC */
508
trant = options->Trans;
514
NCformat *Astore = AA->Store;
515
int nnz = Astore->nnz;
516
int *colptr = Astore->colptr;
517
int *rowind = Astore->rowind;
518
complex *nzval = (complex *)Astore->nzval;
524
t0 = SuperLU_timer_();
525
if ((perm = intMalloc(n)) == NULL)
526
ABORT("SUPERLU_MALLOC fails for perm[]");
528
info1 = cldperm(5, n, nnz, colptr, rowind, nzval, perm, R, C);
530
if (info1 > 0) { /* MC64 fails, call cgsequ() later */
535
for (i = 0; i < n; i++) {
539
/* permute and scale the matrix */
540
for (j = 0; j < n; j++) {
541
for (i = colptr[j]; i < colptr[j + 1]; i++) {
542
cs_mult(&nzval[i], &nzval[i], R[rowind[i]] * C[j]);
543
rowind[i] = perm[rowind[i]];
547
utime[EQUIL] = SuperLU_timer_() - t0;
549
if ( !mc64 & equil ) {
550
t0 = SuperLU_timer_();
551
/* Compute row and column scalings to equilibrate the matrix A. */
552
cgsequ(AA, R, C, &rowcnd, &colcnd, &amax, &info1);
555
/* Equilibrate matrix A. */
556
claqgs(AA, R, C, rowcnd, colcnd, amax, equed);
557
rowequ = lsame_(equed, "R") || lsame_(equed, "B");
558
colequ = lsame_(equed, "C") || lsame_(equed, "B");
560
utime[EQUIL] = SuperLU_timer_() - t0;
565
/* Scale the right hand side if equilibration was performed. */
568
for (j = 0; j < nrhs; ++j)
569
for (i = 0; i < A->nrow; ++i) {
570
cs_mult(&Bmat[i+j*ldb], &Bmat[i+j*ldb], R[i]);
573
} else if ( colequ ) {
574
for (j = 0; j < nrhs; ++j)
575
for (i = 0; i < A->nrow; ++i) {
576
cs_mult(&Bmat[i+j*ldb], &Bmat[i+j*ldb], C[i]);
583
t0 = SuperLU_timer_();
585
* Gnet column permutation vector perm_c[], according to permc_spec:
586
* permc_spec = NATURAL: natural ordering
587
* permc_spec = MMD_AT_PLUS_A: minimum degree on structure of A'+A
588
* permc_spec = MMD_ATA: minimum degree on structure of A'*A
589
* permc_spec = COLAMD: approximate minimum degree column ordering
590
* permc_spec = MY_PERMC: the ordering already supplied in perm_c[]
592
permc_spec = options->ColPerm;
593
if ( permc_spec != MY_PERMC && options->Fact == DOFACT )
594
get_perm_c(permc_spec, AA, perm_c);
595
utime[COLPERM] = SuperLU_timer_() - t0;
597
t0 = SuperLU_timer_();
598
sp_preorder(options, AA, perm_c, etree, &AC);
599
utime[ETREE] = SuperLU_timer_() - t0;
601
/* Compute the LU factorization of A*Pc. */
602
t0 = SuperLU_timer_();
603
cgsitrf(options, &AC, relax, panel_size, etree, work, lwork,
604
perm_c, perm_r, L, U, stat, info);
605
utime[FACT] = SuperLU_timer_() - t0;
608
mem_usage->total_needed = *info - A->ncol;
613
if ( options->PivotGrowth ) {
614
if ( *info > 0 ) return;
616
/* Compute the reciprocal pivot growth factor *recip_pivot_growth. */
617
*recip_pivot_growth = cPivotGrowth(A->ncol, AA, perm_c, L, U);
620
if ( options->ConditionNumber ) {
621
/* Estimate the reciprocal of the condition number of A. */
622
t0 = SuperLU_timer_();
624
*(unsigned char *)norm = '1';
626
*(unsigned char *)norm = 'I';
628
anorm = clangs(norm, AA);
629
cgscon(norm, L, U, anorm, rcond, stat, &info1);
630
utime[RCOND] = SuperLU_timer_() - t0;
634
/* Compute the solution matrix X. */
635
for (j = 0; j < nrhs; j++) /* Save a copy of the right hand sides */
636
for (i = 0; i < B->nrow; i++)
637
Xmat[i + j*ldx] = Bmat[i + j*ldb];
639
t0 = SuperLU_timer_();
640
cgstrs (trant, L, U, perm_c, perm_r, X, stat, &info1);
641
utime[SOLVE] = SuperLU_timer_() - t0;
643
/* Transform the solution matrix X to a solution of the original
647
for (j = 0; j < nrhs; ++j)
648
for (i = 0; i < A->nrow; ++i) {
649
cs_mult(&Xmat[i+j*ldx], &Xmat[i+j*ldx], C[i]);
658
if ((tmp = complexMalloc(n)) == NULL)
659
ABORT("SUPERLU_MALLOC fails for tmp[]");
660
for (j = 0; j < nrhs; j++) {
661
for (i = 0; i < n; i++)
662
tmp[i] = Xmat[i + j * ldx]; /*dcopy*/
663
for (i = 0; i < n; i++)
664
cs_mult(&Xmat[i+j*ldx], &tmp[perm[i]], R[i]);
668
for (j = 0; j < nrhs; ++j)
669
for (i = 0; i < A->nrow; ++i) {
670
cs_mult(&Xmat[i+j*ldx], &Xmat[i+j*ldx], R[i]);
675
} /* end if nrhs > 0 */
677
if ( options->ConditionNumber ) {
678
/* Set INFO = A->ncol+1 if the matrix is singular to working precision. */
679
if ( *rcond < slamch_("E") && *info == 0) *info = A->ncol + 1;
682
if (perm) SUPERLU_FREE(perm);
685
ilu_cQuerySpace(L, U, mem_usage);
686
Destroy_CompCol_Permuted(&AC);
688
if ( A->Stype == SLU_NR ) {
689
Destroy_SuperMatrix_Store(AA);
added
removed
scipy/sparse/linalg/dsolve/SuperLU/SRC/cgsisx.c
