Mercurial > forge

comparison main/sparse/SuperLU/SRC/zgssv.c @ 0:6b33357c7561 octave-forge

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Initial revision
author pkienzle
date Wed, 10 Oct 2001 19:54:49 +0000
parents
children 7dad48fc229c
comparison
equal deleted inserted replaced
-1:000000000000 0:6b33357c7561
1
2
3 /*
4 * -- SuperLU routine (version 2.0) --
5 * Univ. of California Berkeley, Xerox Palo Alto Research Center,
6 * and Lawrence Berkeley National Lab.
7 * November 15, 1997
8 *
9 */
10 #include "zsp_defs.h"
11 #include "util.h"
12
13 void
14 zgssv(SuperMatrix *A, int *perm_c, int *perm_r, SuperMatrix *L,
15 SuperMatrix *U, SuperMatrix *B, int *info )
16 {
17 /*
18 * Purpose
19 * =======
20 *
21 * ZGSSV solves the system of linear equations A*X=B, using the
22 * LU factorization from ZGSTRF. It performs the following steps:
23 *
24 * 1. If A is stored column-wise (A->Stype = NC):
25 *
26 * 1.1. Permute the columns of A, forming A*Pc, where Pc
27 * is a permutation matrix. For more details of this step,
28 * see sp_preorder.c.
29 *
30 * 1.2. Factor A as Pr*A*Pc=L*U with the permutation Pr determined
31 * by Gaussian elimination with partial pivoting.
32 * L is unit lower triangular with offdiagonal entries
33 * bounded by 1 in magnitude, and U is upper triangular.
34 *
35 * 1.3. Solve the system of equations A*X=B using the factored
36 * form of A.
37 *
38 * 2. If A is stored row-wise (A->Stype = NR), apply the
39 * above algorithm to the transpose of A:
40 *
41 * 2.1. Permute columns of transpose(A) (rows of A),
42 * forming transpose(A)*Pc, where Pc is a permutation matrix.
43 * For more details of this step, see sp_preorder.c.
44 *
45 * 2.2. Factor A as Pr*transpose(A)*Pc=L*U with the permutation Pr
46 * determined by Gaussian elimination with partial pivoting.
47 * L is unit lower triangular with offdiagonal entries
48 * bounded by 1 in magnitude, and U is upper triangular.
49 *
50 * 2.3. Solve the system of equations A*X=B using the factored
51 * form of A.
52 *
53 * See supermatrix.h for the definition of 'SuperMatrix' structure.
54 *
55 * Arguments
56 * =========
57 *
58 * A (input) SuperMatrix*
59 * Matrix A in A*X=B, of dimension (A->nrow, A->ncol). The number
60 * of linear equations is A->nrow. Currently, the type of A can be:
61 * Stype = NC or NR; Dtype = _Z; Mtype = GE. In the future, more
62 * general A will be handled.
63 *
64 * perm_c (input/output) int*
65 * If A->Stype = NC, column permutation vector of size A->ncol
66 * which defines the permutation matrix Pc; perm_c[i] = j means
67 * column i of A is in position j in A*Pc.
68 * On exit, perm_c may be overwritten by the product of the input
69 * perm_c and a permutation that postorders the elimination tree
70 * of Pc'*A'*A*Pc; perm_c is not changed if the elimination tree
71 * is already in postorder.
72 *
73 * If A->Stype = NR, column permutation vector of size A->nrow
74 * which describes permutation of columns of transpose(A)
75 * (rows of A) as described above.
76 *
77 * perm_r (output) int*
78 * If A->Stype = NC, row permutation vector of size A->nrow,
79 * which defines the permutation matrix Pr, and is determined
80 * by partial pivoting. perm_r[i] = j means row i of A is in
81 * position j in Pr*A.
82 *
83 * If A->Stype = NR, permutation vector of size A->ncol, which
84 * determines permutation of rows of transpose(A)
85 * (columns of A) as described above.
86 *
87 * L (output) SuperMatrix*
88 * The factor L from the factorization
89 * Pr*A*Pc=L*U (if A->Stype = NC) or
90 * Pr*transpose(A)*Pc=L*U (if A->Stype = NR).
91 * Uses compressed row subscripts storage for supernodes, i.e.,
92 * L has types: Stype = SC, Dtype = _Z, Mtype = TRLU.
93 *
94 * U (output) SuperMatrix*
95 * The factor U from the factorization
96 * Pr*A*Pc=L*U (if A->Stype = NC) or
97 * Pr*transpose(A)*Pc=L*U (if A->Stype = NR).
98 * Uses column-wise storage scheme, i.e., U has types:
99 * Stype = NC, Dtype = _Z, Mtype = TRU.
100 *
101 * B (input/output) SuperMatrix*
102 * B has types: Stype = DN, Dtype = _Z, Mtype = GE.
103 * On entry, the right hand side matrix.
104 * On exit, the solution matrix if info = 0;
105 *
106 * info (output) int*
107 * = 0: successful exit
108 * > 0: if info = i, and i is
109 * <= A->ncol: U(i,i) is exactly zero. The factorization has
110 * been completed, but the factor U is exactly singular,
111 * so the solution could not be computed.
112 * > A->ncol: number of bytes allocated when memory allocation
113 * failure occurred, plus A->ncol.
114 *
115 */
116 double t1; /* Temporary time */
117 char refact[1], trans[1];
118 DNformat *Bstore;
119 SuperMatrix *AA; /* A in NC format used by the factorization routine.*/
120 SuperMatrix AC; /* Matrix postmultiplied by Pc */
121 int lwork = 0, *etree, i;
122
123 /* Set default values for some parameters */
124 double diag_pivot_thresh = 1.0;
125 double drop_tol = 0;
126 int panel_size; /* panel size */
127 int relax; /* no of columns in a relaxed snodes */
128 double *utime;
129 extern SuperLUStat_t SuperLUStat;
130
131 /* Test the input parameters ... */
132 *info = 0;
133 Bstore = B->Store;
134 if ( A->nrow != A->ncol || A->nrow < 0 ||
135 (A->Stype != NC && A->Stype != NR) ||
136 A->Dtype != _Z || A->Mtype != GE )
137 *info = -1;
138 else if ( B->ncol < 0 || Bstore->lda < MAX(0, A->nrow) ||
139 B->Stype != DN || B->Dtype != _Z || B->Mtype != GE )
140 *info = -6;
141 if ( *info != 0 ) {
142 i = -(*info);
143 xerbla_("zgssv", &i);
144 return;
145 }
146
147 *refact = 'N';
148 *trans = 'N';
149 panel_size = sp_ienv(1);
150 relax = sp_ienv(2);
151
152 StatInit(panel_size, relax);
153 utime = SuperLUStat.utime;
154
155 /* Convert A to NC format when necessary. */
156 if ( A->Stype == NR ) {
157 NRformat *Astore = A->Store;
158 AA = (SuperMatrix *) SUPERLU_MALLOC( sizeof(SuperMatrix) );
159 zCreate_CompCol_Matrix(AA, A->ncol, A->nrow, Astore->nnz,
160 Astore->nzval, Astore->colind, Astore->rowptr,
161 NC, A->Dtype, A->Mtype);
162 *trans = 'T';
163 } else if ( A->Stype == NC ) AA = A;
164
165 etree = intMalloc(A->ncol);
166
167 t1 = SuperLU_timer_();
168 sp_preorder(refact, AA, perm_c, etree, &AC);
169 utime[ETREE] = SuperLU_timer_() - t1;
170
171 /*printf("Factor PA = LU ... relax %d\tw %d\tmaxsuper %d\trowblk %d\n",
172 relax, panel_size, sp_ienv(3), sp_ienv(4));*/
173 t1 = SuperLU_timer_();
174 /* Compute the LU factorization of A. */
175 zgstrf(refact, &AC, diag_pivot_thresh, drop_tol, relax, panel_size,
176 etree, NULL, lwork, perm_r, perm_c, L, U, info);
177 utime[FACT] = SuperLU_timer_() - t1;
178
179 t1 = SuperLU_timer_();
180 if ( *info == 0 ) {
181 /* Solve the system A*X=B, overwriting B with X. */
182 zgstrs (trans, L, U, perm_r, perm_c, B, info);
183 }
184 utime[SOLVE] = SuperLU_timer_() - t1;
185
186 SUPERLU_FREE (etree);
187 Destroy_CompCol_Permuted(&AC);
188 if ( A->Stype == NR ) {
189 Destroy_SuperMatrix_Store(AA);
190 SUPERLU_FREE(AA);
191 }
192
193 PrintStat( &SuperLUStat );
194 StatFree();
195
196 }