1 /* ZMLBiCGSTABRg.c */
2
3
4 #include "../Iter.h"
5
6 /*--------------------------------------------------------------------*/
7 /*
8 ---------------------------------------------------------------------
9 purpose -- to solve a complex matrix equation
10
11 Ax=b
12
13 using right preconditioned ML(k)BiCGSTABR method
14 by M-C Yeung and T. Chan
15
16
17 x -- Initial guess as zeros
18 A -- Input matrix
19 Precond -- Front Matrix as the preconditioner
20 Q -- Starting vectors
21 b -- Right-hand side
22 tol -- Convergence tolerance
23 type -- type of entries
24 SPOOLES_REAL or SPOOLES_COMPLEX
25 symmetryflag -- symmetry of the matrix
26 SPOOLES_SYMMETRIC, SPOOLES_HERMITIAN or SPOOLES_NONSYMMETRIC
27 nrhs -- number of right hand sides
28 msglvl -- message level
29 msgFile -- message file
30
31 return value -- error flag
32
33 created -- Dec. 10, 1998 Wei-Pai Tang
34 ---------------------------------------------------------------------
35 */
36
37 int
zmlbicgstabr(int n_matrixSize,int type,int symmetryflag,InpMtx * mtxA,FrontMtx * Precond,DenseMtx * mtxQ,DenseMtx * mtxX,DenseMtx * mtxB,int itermax,double convergetol,int msglvl,FILE * msgFile)38 zmlbicgstabr (
39 int n_matrixSize,
40 int type,
41 int symmetryflag,
42 InpMtx *mtxA,
43 FrontMtx *Precond,
44 DenseMtx *mtxQ,
45 DenseMtx *mtxX,
46 DenseMtx *mtxB,
47 int itermax,
48 double convergetol,
49 int msglvl,
50 FILE *msgFile
51 )
52 {
53 Chv *chv, *rootchv ;
54 ChvManager *chvmanager ;
55 DenseMtx *mtxD, *mtxG, *mtxU, *mtxW;
56 DenseMtx *vecGt, *vecR, *vecT, *vecT1;
57 DenseMtx *vecX, *vecZD, *vecZG, *vecZW ;
58 double Alpha[2], Beta[2], Rho[2] ;
59 double c[200] ;
60 double Init_norm, ratio, Res_norm;
61 double error_trol, m, Rtmp[2], Ttmp[2];
62 double t1, t2, cpus[9] ;
63 double one[2] = {1.0, 0.0} ;
64 double zero[2] = {0.0, 0.0} ;
65 double minusone[2] = {-1.0, 0.0};
66 double Tiny = 0.1e-28;
67 int Iter, Imv, neqns, Ik, ii, is;
68 int stats[6] ;
69 int return_flag;
70
71
72
73 neqns = n_matrixSize;
74 Ik = mtxQ->ncol;
75
76 if (Ik > 99){
77 fprintf(msgFile, "\n\n Fatal Error, \n"
78 " Too many starting vectors in Q !!") ;
79 return(-1);
80 };
81
82 return_flag = 1;
83
84 /*
85 --------------------
86 init the vectors in ZMLBiCGSTABR
87 --------------------
88 */
89 mtxD = DenseMtx_new() ;
90 DenseMtx_init(mtxD, type, 0, 0, neqns, Ik, 1, neqns) ;
91
92 mtxG = DenseMtx_new() ;
93 DenseMtx_init(mtxG, type, 0, 0, neqns, Ik, 1, neqns) ;
94
95 mtxU = DenseMtx_new() ;
96 DenseMtx_init(mtxU, type, 0, 0, neqns, 2, 1, neqns) ;
97
98 mtxW = DenseMtx_new() ;
99 DenseMtx_init(mtxW, type, 0, 0, neqns, Ik, 1, neqns) ;
100
101 vecGt = DenseMtx_new() ;
102 DenseMtx_init(vecGt, type, 0, 0, neqns, 1, 1, neqns) ;
103
104 vecR = DenseMtx_new() ;
105 DenseMtx_init(vecR, type, 0, 0, neqns, 1, 1, neqns) ;
106
107 vecT = DenseMtx_new() ;
108 DenseMtx_init(vecT, type, 0, 0, neqns, 1, 1, neqns) ;
109
110 vecT1 = DenseMtx_new() ;
111 DenseMtx_init(vecT1, type, 0, 0, neqns, 1, 1, neqns) ;
112
113 vecX = DenseMtx_new() ;
114 DenseMtx_init(vecX, type, 0, 0, neqns, 1, 1, neqns) ;
115
116
117 vecZD = DenseMtx_new() ;
118 DenseMtx_init(vecZD, type, 0, 0, neqns, 1, 1, neqns) ;
119
120 vecZG = DenseMtx_new() ;
121 DenseMtx_init(vecZG, type, 0, 0, neqns, 1, 1, neqns) ;
122
123 vecZW = DenseMtx_new() ;
124 DenseMtx_init(vecZW, type, 0, 0, neqns, 1, 1, neqns) ;
125
126
127 for ( ii = 0; ii <200; ii++){
128 c[ii] = 0;
129 }
130 /*
131 --------------------------
132 Initialize the iterations
133 --------------------------
134 */
135 /* ---- Set initial guess as zero ---- */
136 DenseMtx_zero(vecX) ;
137
138 /* ---- If x_0 is not zero ---- */
139 /*
140 DenseMtx_colCopy (vecX, 0, mtxX, 0);
141 */
142
143
144 switch ( symmetryflag ) {
145 case SPOOLES_SYMMETRIC :
146 InpMtx_sym_gmmm(mtxA, zero, vecT, one, vecX) ;
147 break ;
148 case SPOOLES_HERMITIAN :
149 InpMtx_herm_gmmm(mtxA, zero, vecT, one, vecX) ;
150 break ;
151 case SPOOLES_NONSYMMETRIC :
152 InpMtx_nonsym_gmmm(mtxA, zero, vecT, one, vecX) ;
153 break ;
154 default :
155 fprintf(msgFile, "\n ZTFQMR Matrix type wrong");
156 fprintf(msgFile, "\n Fatal error");
157 goto end;
158 }
159
160 DenseMtx_colCopy(vecR, 0, mtxB, 0);
161 DenseMtx_sub(vecR, vecT) ;
162
163
164 Init_norm = DenseMtx_twoNormOfColumn(vecR, 0);
165 if ( Init_norm == 0.0 ){
166 Init_norm = 1.0;
167 };
168 error_trol = Init_norm * convergetol ;
169
170 fprintf(msgFile, "\n ZMLBiCGSTABR Initial norml: %6.2e ", Init_norm ) ;
171 fprintf(msgFile, "\n ZMLBiCGSTABR Conveg. Control: %7.3e ", convergetol ) ;
172 fprintf(msgFile, "\n ZMLBiCGSTABR Convergen Control: %7.3e ",error_trol ) ;
173
174 DenseMtx_zero(mtxG) ;
175 DenseMtx_zero(mtxD) ;
176 DenseMtx_zero(mtxW) ;
177
178
179 Iter = 0;
180 Imv = 0;
181
182
183 DenseMtx_colCopy (mtxG, Ik-1, vecR, 0);
184
185 /*
186 ------------------------------
187 ZMLBiCGSTABR Iteration start
188 ------------------------------
189 */
190
191 MARKTIME(t1) ;
192
193
194 while ( Iter <= itermax ){
195
196 Iter++;
197
198 FrontMtx_solveOneColumn(Precond, vecGt, 0, mtxG, Ik-1,
199 Precond->manager, cpus, msglvl, msgFile) ;
200
201 switch ( symmetryflag ) {
202 case SPOOLES_SYMMETRIC :
203 InpMtx_sym_gmmm(mtxA, zero, vecT, one, vecGt) ;
204 break ;
205 case SPOOLES_HERMITIAN :
206 InpMtx_herm_gmmm(mtxA, zero, vecT, one, vecGt) ;
207 break ;
208 case SPOOLES_NONSYMMETRIC :
209 InpMtx_nonsym_gmmm(mtxA, zero, vecT, one, vecGt) ;
210 break ;
211 default :
212 fprintf(msgFile, "\n ZTFQMR Matrix type wrong");
213 fprintf(msgFile, "\n Fatal error");
214 goto end;
215 }
216 DenseMtx_colCopy (mtxW, Ik-1, vecT, 0);
217 Imv++;
218
219 /* c[Ik] = *DenseMtx_colDotProduct (mtxQ, 0, mtxG, Ik-1); */
220 DenseMtx_colDotProduct (mtxQ, 0, mtxG, Ik-1, Rtmp);
221 c[2*Ik] = Rtmp[0];
222 c[2*Ik+1]= Rtmp[1];
223
224 if (zabs(Rtmp) == 0){
225 fprintf(msgFile, "\n\n Fatal Error, \n"
226 " ZMLBiCGSTABR Breakdown, c[k] = 0 !!") ;
227 return_flag = -1;
228 goto end;
229 };
230
231 /* Alpha =( *DenseMtx_colDotProduct (mtxQ,0, vecR, 0))/c[Ik] ; */
232
233 DenseMtx_colDotProduct (mtxQ,0, vecR, 0, Rtmp);
234 /* Alpha = Rtmp/c[Ik]; */
235 zdiv(Rtmp, c+(2*Ik), Alpha);
236 DenseMtx_colCopy (mtxU, 0, vecR, 0);
237 /* Rtmp = -Alpha; */
238 zsub(zero, Alpha, Rtmp);
239 DenseMtx_colGenAxpy (one, mtxU, 0, Rtmp, mtxW, Ik-1);
240
241
242 FrontMtx_solveOneColumn(Precond, vecT1, 0, mtxU, 0,
243 Precond->manager, cpus, msglvl, msgFile) ;
244
245 DenseMtx_colCopy (mtxU, 1, vecT1, 0);
246 switch ( symmetryflag ) {
247 case SPOOLES_SYMMETRIC :
248 InpMtx_sym_gmmm(mtxA, zero, vecT, one, vecT1) ;
249 break ;
250 case SPOOLES_HERMITIAN :
251 InpMtx_herm_gmmm(mtxA, zero, vecT, one, vecT1) ;
252 break ;
253 case SPOOLES_NONSYMMETRIC :
254 InpMtx_nonsym_gmmm(mtxA, zero, vecT, one, vecT1) ;
255 break ;
256 default :
257 fprintf(msgFile, "\n ZTFQMR Matrix type wrong");
258 fprintf(msgFile, "\n Fatal error");
259 goto end;
260 }
261
262 Imv++;
263
264 /* Rho = *DenseMtx_colDotProduct (vecT, 0, vecT, 0); */
265 DenseMtx_colDotProduct (vecT, 0, vecT, 0, Rho);
266 if (zabs(Rho) == 0){
267 fprintf(msgFile, "\n\n Fatal Error, \n"
268 " ZMLBiCGSTABR Breakdown, Rho = 0 !!") ;
269 return_flag = -1;
270 goto end;
271 };
272
273 /* Rho = -(*DenseMtx_colDotProduct (mtxU, 0, vecT, 0))/Rho; */
274 DenseMtx_colDotProduct (mtxU, 0, vecT, 0, Rtmp);
275 /* Rho = -Rtmp/Rho; */
276 zdiv(Rtmp, Rho, Ttmp);
277 zsub(zero, Ttmp, Rho);
278
279 DenseMtx_colGenAxpy(one, vecX, 0, Alpha, vecGt, 0);
280 /* Rtmp = -Rho; */
281 zsub(zero, Rho, Rtmp);
282 DenseMtx_colGenAxpy (one, vecX, 0, Rtmp, mtxU, 1) ;
283 DenseMtx_colCopy (vecR, 0, mtxU, 0);
284 DenseMtx_colGenAxpy (one, vecR, 0, Rho, vecT, 0) ;
285 /* Iter++; */
286
287 /*
288 ----------------
289 Convergence Test
290 ---------------
291 */
292 Rtmp[0] = DenseMtx_twoNormOfColumn ( vecR, 0);
293 ratio = Rtmp[0]/Init_norm;
294 fprintf(msgFile, "\n\n At iteration %d"
295 " the convergence ratio is %12.4e"
296 "\n Residual norm is %6.2ee",
297 Imv, ratio, Rtmp[0]) ;
298 fflush(msgFile) ;
299 if ( Rtmp[0] <= error_trol ) {
300 fprintf(msgFile, "\n # iterations = %d", Imv) ;
301 return_flag = Imv;
302 goto end;
303 };
304
305 for (ii = 1; ii < Ik+1; ii++ ){
306 if (Iter > itermax ){
307 fprintf(msgFile, "\n # iterations = %d", Imv) ;
308 fprintf(msgFile, "\n\n ZMLBiCGSTABR did not Converge !") ;
309 return_flag = Imv;
310 goto end;
311 };
312
313
314 DenseMtx_colCopy (vecZD,0, mtxU, 0);
315 DenseMtx_colCopy (vecZG,0, vecR, 0);
316 DenseMtx_zero(vecZW) ;
317
318 if ( Iter > 1 ){
319 for ( is = ii ; is < Ik ; is++ ){
320 /* Beta =-( *(DenseMtx_colDotProduct(mtxQ, is, vecZD, 0)))/c[is]; */
321 DenseMtx_colDotProduct(mtxQ, is, vecZD, 0, Rtmp);
322 /* Beta = -Rtmp/c[is]; */
323 zdiv(Rtmp, c+(2*is), Ttmp);
324 zsub(zero, Ttmp, Beta);
325 DenseMtx_colGenAxpy (one, vecZD, 0, Beta, mtxD, is-1);
326 DenseMtx_colGenAxpy (one, vecZG, 0, Beta, mtxG, is-1);
327 DenseMtx_colGenAxpy (one, vecZW, 0, Beta, mtxW, is-1);
328 };
329 };
330
331 /* Beta = Rho * c[Ik]; */
332 zmul(Rho, c+(2*Ik), Beta);
333
334 if (zabs(Beta) == 0){
335 fprintf(msgFile, "\n\n Fatal Error, \n"
336 " ZMLBiCGSTABR Breakdown, Beta = 0 !!") ;
337 return_flag = -1;
338 goto end;
339 };
340 /*
341 Beta = - Q(:,1)'* (r + Rho* zw )/ Beta;
342 */
343 DenseMtx_colCopy (vecT, 0, vecR, 0);
344 DenseMtx_colGenAxpy (one, vecT, 0, Rho, vecZW, 0);
345 /* Beta = - (*DenseMtx_colDotProduct(mtxQ, 0, vecT, 0))/Beta; */
346 DenseMtx_colDotProduct(mtxQ, 0, vecT, 0, Rtmp);
347 /* Beta = - Rtmp/Beta; */
348 zdiv(Rtmp, Beta, Ttmp);
349 zsub(zero, Ttmp, Beta);
350 /*
351 zg = zg + beta*G(:,k);
352 zw = rho*(zw + beta*W(:,k));
353 zd = r + zw;
354 */
355 DenseMtx_colGenAxpy (one, vecZG, 0, Beta, mtxG, Ik-1);
356 /* Rtmp = Rho*Beta; */
357 zmul(Rho, Beta, Rtmp);
358 DenseMtx_colGenAxpy (Rho, vecZW, 0, Rtmp, mtxW, Ik-1);
359 DenseMtx_colCopy (vecZD, 0, vecR, 0);
360 DenseMtx_colGenAxpy (one, vecZD, 0, one, vecZW, 0);
361 /*
362
363 for s = 1:i-1
364 beta = -Q(:,s+1)'*zd/c(s);
365 zd = zd + beta*D(:,s);
366 zg = zg + beta*G(:,s);
367 end
368 */
369 for ( is = 1; is < ii - 1; is ++){
370 /* Beta = -(*DenseMtx_colDotProduct(mtxQ, is, vecZD, 0))/c[is] ; */
371 DenseMtx_colDotProduct(mtxQ, is, vecZD, 0, Rtmp);
372 /* Beta = - Rtmp/c[is]; */
373 zdiv(Rtmp, c+(2*is), Ttmp);
374 zsub(zero, Ttmp, Beta);
375 DenseMtx_colGenAxpy (one, vecZD, 0, Beta, mtxD, is-1);
376 DenseMtx_colGenAxpy (one, vecZG, 0, Beta, mtxG, is-1);
377 };
378
379 /*
380 D(:,i) = zd - u;
381 G(:,i) = zg + zw;
382 */
383 DenseMtx_colCopy (mtxD, ii-1, vecZD, 0);
384 DenseMtx_colGenAxpy (one, mtxD, ii-1, minusone, mtxU, 0);
385 DenseMtx_colCopy (mtxG, ii-1, vecZG, 0);
386 DenseMtx_colGenAxpy (one, mtxG, ii-1, one, vecZW, 0);
387
388 /*
389 if i < k
390 c(i) = Q(:,i+1)'*D(:,i);
391 */
392 if ( ii < Ik ){
393 /* c[ii] = *DenseMtx_colDotProduct(mtxQ, ii, mtxD, ii-1); */
394 DenseMtx_colDotProduct(mtxQ, ii, mtxD, ii-1, Rtmp);
395 c[2*ii] = Rtmp[0];
396 c[2*ii+1] = Rtmp[1];
397 /*
398 If breakdown ?
399 */
400 if (zabs(Rtmp) == 0){
401 fprintf(msgFile, "\n\n Fatal Error, \n"
402 " ZMLBiCGSTABR Breakdown, c[ii] = 0 !!") ;
403 return_flag = -1;
404 goto end;
405 };
406 /*
407 alpha = Q(:,i+1)'*u/c(i);
408 u = u - alpha*D(:,i);
409 g_tld = U\(L\G(:,i));
410 */
411 /* Alpha = (*DenseMtx_colDotProduct(mtxQ, ii, mtxU, 0))/c[ii]; */
412 DenseMtx_colDotProduct(mtxQ, ii, mtxU, 0, Rtmp);
413 /* Alpha = Rtmp/c[ii]; */
414 zdiv(Rtmp, c+(2*ii), Alpha);
415
416 /* Rtmp = -Alpha; */
417 zsub(zero, Alpha, Rtmp);
418 DenseMtx_colGenAxpy (one, mtxU, 0, Rtmp, mtxD, ii-1);
419 /*
420 DenseMtx_colCopy (vecT, 0, mtxG, ii-1);
421 FrontMtx_solve(Precond, vecGt, vecT, Precond->manager,
422 cpus, msglvl, msgFile) ;
423 */
424 FrontMtx_solveOneColumn(Precond, vecGt, 0, mtxG, ii-1,
425 Precond->manager, cpus, msglvl, msgFile) ;
426
427 /*
428 x = x + rho*alpha*g_tld;
429 W(:,i) = A*g_tld;
430 r = r - rho*alpha*W(:,i);
431 */
432 /* Rtmp = Rho * Alpha; */
433 zmul(Rho, Alpha, Rtmp);
434 DenseMtx_colGenAxpy (one, vecX, 0, Rtmp, vecGt, 0);
435 switch ( symmetryflag ) {
436 case SPOOLES_SYMMETRIC :
437 InpMtx_sym_gmmm(mtxA, zero, vecT, one, vecGt) ;
438 break ;
439 case SPOOLES_HERMITIAN :
440 InpMtx_herm_gmmm(mtxA, zero, vecT, one, vecGt) ;
441 break ;
442 case SPOOLES_NONSYMMETRIC :
443 InpMtx_nonsym_gmmm(mtxA, zero, vecT, one, vecGt) ;
444 break ;
445 default :
446 fprintf(msgFile, "\n ZTFQMR Matrix type wrong");
447 fprintf(msgFile, "\n Fatal error");
448 goto end;
449 }
450 Imv++;
451 DenseMtx_colCopy (mtxW, ii-1, vecT, 0);
452 /* Rtmp = -Rtmp; */
453 zsub(zero, Rtmp, Rtmp);
454 DenseMtx_colGenAxpy (one, vecR, 0, Rtmp, mtxW, ii-1);
455
456 /*
457 ----------------
458 Convergence Test
459 ---------------
460 */
461 Rtmp[0] = DenseMtx_twoNormOfColumn ( vecR, 0);
462 if ( Rtmp[0] <= error_trol ) {
463 fprintf(msgFile, "\n # iterations = %d", Imv) ;
464 return_flag = Imv;
465 goto end;
466 };
467 };
468
469 };
470
471 }
472 /* End of while loop */
473 MARKTIME(t2) ;
474 fprintf(msgFile, "\n CPU : Total iteration time is : %8.3f ", t2 - t1) ;
475 fprintf(msgFile, "\n # iterations = %d", Imv) ;
476 fprintf(msgFile, "\n\n ZMLBiCGSTABR did not Converge !") ;
477 DenseMtx_colCopy(mtxX, 0, vecX, 0);
478
479 /*
480
481 ------------------------
482 free the working storage
483 ------------------------
484 */
485 end:
486 MARKTIME(t2) ;
487 fprintf(msgFile, "\n CPU : Total iteration time is : %8.3f ", t2 - t1) ;
488 DenseMtx_colCopy(mtxX, 0, vecX, 0);
489 switch ( symmetryflag ) {
490 case SPOOLES_SYMMETRIC :
491 InpMtx_sym_gmmm(mtxA, zero, vecT, one, mtxX) ;
492 break ;
493 case SPOOLES_HERMITIAN :
494 InpMtx_herm_gmmm(mtxA, zero, vecT, one, mtxX) ;
495 break ;
496 case SPOOLES_NONSYMMETRIC :
497 InpMtx_nonsym_gmmm(mtxA, zero, vecT, one, mtxX) ;
498 break ;
499 default :
500 fprintf(msgFile, "\n ZTFQMR Matrix type wrong");
501 fprintf(msgFile, "\n Fatal error");
502 goto end;
503 }
504
505 DenseMtx_sub(vecT, mtxB) ;
506
507 Rtmp[0] = DenseMtx_twoNormOfColumn(vecT, 0);
508
509
510 fprintf(msgFile, "\n ZMLBiCGSTABR True Residual norm: %6.2e ",
511 Rtmp[0]) ;
512 fprintf(msgFile, "\n\n after ZMLBiCGSTABR") ;
513 DenseMtx_free(mtxD) ;
514 DenseMtx_free(mtxG) ;
515 DenseMtx_free(mtxU) ;
516 DenseMtx_free(mtxW) ;
517 DenseMtx_free(vecGt) ;
518 DenseMtx_free(vecR) ;
519 DenseMtx_free(vecT) ;
520 DenseMtx_free(vecT1) ;
521 DenseMtx_free(vecX) ;
522 DenseMtx_free(vecZD) ;
523 DenseMtx_free(vecZG) ;
524 DenseMtx_free(vecZW) ;
525
526
527 fprintf(msgFile, "\n") ;
528
529 return(return_flag) ; }
530
531 /*--------------------------------------------------------------------*/
532