1 /* zbicgstabl.c */
2
3 #include "../Iter.h"
4
5 /*--------------------------------------------------------------------*/
6 /*
7 ---------------------------------------------------------------------
8 purpose -- to solve a complex matrix equation
9
10 Ax=b
11 using left preconditioned BiCGSTABL
12
13 x -- Initial guess
14 A -- Input matrix
15 M -- Front Matrix as the preconditioner
16 b -- Right-hand side
17 tol -- Convergence tolerance
18 type -- type of entries
19 SPOOLES_REAL or SPOOLES_COMPLEX
20 symmetryflag -- symmetry of the matrix
21 SPOOLES_SYMMETRIC, SPOOLES_HERMITIAN or SPOOLES_NONSYMMETRIC
22 nrhs -- number of right hand sides
23 msglvl -- message level
24 msgFile -- message file
25
26 return value -- error flag
27
28 created -- Dec. 10, 1998 Wei-Pai Tang
29 ---------------------------------------------------------------------
30 */
31
32 int
zbicgstabl(int n_matrixSize,int type,int symmetryflag,InpMtx * mtxA,FrontMtx * Precond,DenseMtx * mtxX,DenseMtx * mtxB,int itermax,double convergetol,int msglvl,FILE * msgFile)33 zbicgstabl (
34 int n_matrixSize,
35 int type,
36 int symmetryflag,
37 InpMtx *mtxA,
38 FrontMtx *Precond,
39 DenseMtx *mtxX,
40 DenseMtx *mtxB,
41 int itermax,
42 double convergetol,
43 int msglvl,
44 FILE *msgFile
45 )
46 {
47 Chv *chv, *rootchv ;
48 ChvManager *chvmanager ;
49 DenseMtx *mtxZ ;
50 DenseMtx *vecP, *vecR, *vecR0, *vecT, *vecV, *vecW;
51 DenseMtx *vecX, *vecY, *vecZ ;
52 double Alpha[2], Beta[2], Init_norm, Omega[2], Rtmp[2];
53 double ratio, Res_norm, Rho[2], Rho_new[2], Ttmp[2];
54 double t1, t2, cpus[9] ;
55 double one[2] = {1.0, 0.0}, zero[2] = {0.0, 0.0} ;
56 double Tiny = 0.1e-28;
57 int Iter, Imv, neqns;
58 int stats[6] ;
59
60
61
62 neqns = n_matrixSize;
63
64 /*
65 --------------------
66 init the vectors in ZBICGSTABL
67 --------------------
68 */
69 vecP = DenseMtx_new() ;
70 DenseMtx_init(vecP, type, 0, 0, neqns, 1, 1, neqns) ;
71
72 vecR = DenseMtx_new() ;
73 DenseMtx_init(vecR, type, 0, 0, neqns, 1, 1, neqns) ;
74
75 vecR0 = DenseMtx_new() ;
76 DenseMtx_init(vecR0, type, 0, 0, neqns, 1, 1, neqns) ;
77
78 vecT = DenseMtx_new() ;
79 DenseMtx_init(vecT, type, 0, 0, neqns, 1, 1, neqns) ;
80
81 vecV = DenseMtx_new() ;
82 DenseMtx_init(vecV, type, 0, 0, neqns, 1, 1, neqns) ;
83
84 vecW = DenseMtx_new() ;
85 DenseMtx_init(vecW, type, 0, 0, neqns, 1, 1, neqns) ;
86
87 vecX = DenseMtx_new() ;
88 DenseMtx_init(vecX, type, 0, 0, neqns, 1, 1, neqns) ;
89
90 vecY = DenseMtx_new() ;
91 DenseMtx_init(vecY, type, 0, 0, neqns, 1, 1, neqns) ;
92
93 vecZ = DenseMtx_new() ;
94 DenseMtx_init(vecZ, type, 0, 0, neqns, 1, 1, neqns) ;
95
96
97
98 /*
99 --------------------------
100 Initialize the iterations
101 --------------------------
102 */
103 /* ---- Set initial guess as zero ---- */
104 DenseMtx_zero(vecX) ;
105
106 DenseMtx_colCopy (vecT, 0, mtxB, 0);
107 /* */
108 FrontMtx_solve(Precond, vecR, vecT, Precond->manager,
109 cpus, msglvl, msgFile) ;
110 /* */
111
112
113
114
115 Init_norm = DenseMtx_twoNormOfColumn(vecR, 0);
116 if ( Init_norm == 0.0 ){
117 Init_norm = 1.0;
118 };
119
120 ratio = 1.0;
121
122 DenseMtx_colCopy (vecR0, 0, vecR, 0);
123
124 fprintf(msgFile, "\n ZBiCGSTABL Initial norm: %6.2e ", Init_norm ) ;
125 fprintf(msgFile, "\n ZBiCGSTABL Conveg. Control: %6.2e ", convergetol ) ;
126
127 /*
128
129 */
130 Rho[0] = 1.0;
131 Alpha[0] = 1.0;
132 Omega[0] = 1.0;
133 Rho[1] = 0.0;
134 Alpha[1] = 0.0;
135 Omega[1] = 0.0;
136
137 DenseMtx_zero(vecV) ;
138 DenseMtx_zero(vecP) ;
139 /*
140 ------------------------------
141
142 ------------------------------
143 */
144
145 MARKTIME(t1) ;
146 Iter = 0;
147 Imv = 0;
148
149 /*
150 -----------------
151 factor the matrix
152 -----------------
153 */
154 while ( ratio > convergetol && Iter <= itermax )
155 {
156 Iter++;
157 DenseMtx_colDotProduct (vecR0, 0, vecR,0, Rho_new);
158 if ( zabs(Rho_new) == 0.0 || zabs(Omega)== 0.0 ){
159 fprintf(stderr, "\n breakdown in ZBiCGSTABL !! "
160 "\n Fatal error \n");
161 exit(-1) ;
162 }
163
164 /* Beta = (Rho_new / (Rho+Tiny)) * (Alpha / (Omega+Tiny)); */
165 zdiv(Rho_new, Rho, Beta);
166 zmul(Beta, Alpha, Rtmp);
167 zdiv(Rtmp, Omega, Beta);
168
169
170 Rho[0] = Rho_new[0];
171 Rho[1] = Rho_new[1];
172 /*
173 DenseMtx_axpy(vecP, vecV, -Omega);
174 DenseMtx_aypx(vecP, vecR, Beta);
175 */
176 zsub(zero, Omega, Rtmp);
177 DenseMtx_colGenAxpy (one, vecP, 0, Rtmp, vecV, 0 );
178 DenseMtx_colGenAxpy (Beta, vecP, 0, one, vecR, 0 );
179
180
181
182 switch ( symmetryflag ) {
183 case SPOOLES_SYMMETRIC :
184 InpMtx_sym_gmmm(mtxA, zero, vecY, one, vecP) ;
185 break ;
186 case SPOOLES_HERMITIAN :
187 InpMtx_herm_gmmm(mtxA, zero, vecY, one, vecP) ;
188 break ;
189 case SPOOLES_NONSYMMETRIC :
190 InpMtx_nonsym_gmmm(mtxA, zero, vecY, one, vecP) ;
191 break ;
192 default :
193 fprintf(msgFile, "\n BiCGSTABL Matrix type wrong");
194 fprintf(msgFile, "\n Fatal error");
195 goto end;
196 }
197
198 /* */
199 FrontMtx_solve(Precond, vecV, vecY, Precond->manager,
200 cpus, msglvl, msgFile) ;
201 Imv++;
202 /* */
203 /* Alpha = Rho / (DenseMtx_dot(vecR0, vecV)+Tiny); */
204
205 DenseMtx_colDotProduct (vecR0, 0, vecV,0, Rtmp);
206 if ( zabs(Rtmp) == 0.0 ){
207 fprintf(stderr, "\n breakdown in ZBiCGSTABL !! "
208 "\n Fatal error \n");
209 exit(-1) ;
210 }
211 zdiv(Rho, Rtmp, Alpha);
212
213
214 /* */
215 /* DenseMtx_axpy(vecR, vecV, -Alpha); */
216
217 zsub(zero, Alpha, Rtmp);
218 DenseMtx_colGenAxpy (one, vecR, 0, Rtmp, vecV, 0 );
219
220 switch ( symmetryflag ) {
221 case SPOOLES_SYMMETRIC :
222 InpMtx_sym_gmmm(mtxA, zero, vecZ, one, vecR) ;
223 break ;
224 case SPOOLES_HERMITIAN :
225 InpMtx_herm_gmmm(mtxA, zero, vecZ, one, vecR) ;
226 break ;
227 case SPOOLES_NONSYMMETRIC :
228 InpMtx_nonsym_gmmm(mtxA, zero, vecZ, one, vecR) ;
229 break ;
230 default :
231 fprintf(msgFile, "\n BiCGSTABL Matrix type wrong");
232 fprintf(msgFile, "\n Fatal error");
233 goto end;
234 }
235
236 /* */
237 FrontMtx_solve(Precond, vecT, vecZ, Precond->manager,
238 cpus, msglvl, msgFile) ;
239 Imv++;
240 /* */
241
242 /* Omega = DenseMtx_dot(vecT, vecR)/(DenseMtx_dot(vecT, vecT)+Tiny); */
243 DenseMtx_colDotProduct (vecT, 0, vecT,0, Ttmp);
244 if ( zabs(Ttmp) == 0.0 ){
245 fprintf(stderr, "\n breakdown in ZBiCGSTABL !! "
246 "\n Fatal error \n");
247 exit(-1) ;
248 };
249 DenseMtx_colDotProduct (vecT, 0, vecR, 0, Rtmp);
250 zdiv(Rtmp, Ttmp, Omega);
251
252
253 DenseMtx_colGenAxpy (one, vecX, 0, Alpha, vecP, 0);
254 DenseMtx_colGenAxpy (one, vecX, 0, Omega, vecR, 0);
255
256 /* */
257 /* DenseMtx_axpy(vecR, vecT, -Omega); */
258 zsub(zero, Omega, Rtmp);
259 DenseMtx_colGenAxpy (one, vecR, 0, Rtmp, vecT, 0);
260
261 Res_norm = DenseMtx_twoNormOfColumn (vecR, 0);
262 ratio = Res_norm/Init_norm;
263 fprintf(msgFile, "\n\n At iteration %d"
264 " the convergence ratio is %12.4e"
265 "\n Residual norm is %6.2e",
266 Imv, ratio, Res_norm) ;
267 }
268 /* End of while loop */
269 MARKTIME(t2) ;
270 fprintf(msgFile, "\n CPU : Converges in time: %8.3f ", t2 - t1) ;
271 fprintf(msgFile, "\n # iterations = %d", Imv) ;
272
273 fprintf(msgFile, "\n\n after ZBICGSTABL") ;
274 /* DenseMtx_copy(mtxX, vecX); */
275 DenseMtx_colCopy (mtxX, 0, vecX, 0);
276
277 /*
278
279 ------------------------
280 free the working storage
281 ------------------------
282 */
283 end:
284 DenseMtx_free(vecP) ;
285 DenseMtx_free(vecR) ;
286 DenseMtx_free(vecR0) ;
287 DenseMtx_free(vecT) ;
288 DenseMtx_free(vecV) ;
289 DenseMtx_free(vecW) ;
290 DenseMtx_free(vecX) ;
291 DenseMtx_free(vecY) ;
292 DenseMtx_free(vecZ) ;
293
294 fprintf(msgFile, "\n") ;
295
296 return(1) ; }
297
298 /*--------------------------------------------------------------------*/
299