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