|
0
|
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 } |
