1 /* zbicgstabr.c */
2
3 #include "../Iter.h"
4
5 /*--------------------------------------------------------------------*/
6 /*
7 ---------------------------------------------------------------------
8 purpose -- to solve a complex matrix equation
9
10 Ax=b
11 using right preconditioned BiCGSTABR
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
zbicgstabr(int n_matrixSize,int type,int symmetryflag,InpMtx * mtxA,FrontMtx * Precond,DenseMtx * mtxX,DenseMtx * mtxB,int itermax,double convergetol,int msglvl,FILE * msgFile)33 zbicgstabr (
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 --------------------
67 init the vectors in bicgstab
68 --------------------
69 */
70 vecP = DenseMtx_new() ;
71 DenseMtx_init(vecP, type, 0, 0, neqns, 1, 1, neqns) ;
72
73 vecR = DenseMtx_new() ;
74 DenseMtx_init(vecR, type, 0, 0, neqns, 1, 1, neqns) ;
75
76 vecR0 = DenseMtx_new() ;
77 DenseMtx_init(vecR0, type, 0, 0, neqns, 1, 1, neqns) ;
78
79 vecT = DenseMtx_new() ;
80 DenseMtx_init(vecT, type, 0, 0, neqns, 1, 1, neqns) ;
81
82 vecV = DenseMtx_new() ;
83 DenseMtx_init(vecV, type, 0, 0, neqns, 1, 1, neqns) ;
84
85 vecW = DenseMtx_new() ;
86 DenseMtx_init(vecW, type, 0, 0, neqns, 1, 1, neqns) ;
87
88 vecX = DenseMtx_new() ;
89 DenseMtx_init(vecX, type, 0, 0, neqns, 1, 1, neqns) ;
90
91 vecY = DenseMtx_new() ;
92 DenseMtx_init(vecY, type, 0, 0, neqns, 1, 1, neqns) ;
93
94 vecZ = DenseMtx_new() ;
95 DenseMtx_init(vecZ, type, 0, 0, neqns, 1, 1, neqns) ;
96
97
98
99 /*
100 --------------------------
101 Initialize the iterations
102 --------------------------
103 */
104 Init_norm = DenseMtx_twoNormOfColumn(mtxB, 0);
105
106 ratio = 1.0;
107 DenseMtx_zero(vecX) ;
108
109
110 DenseMtx_colCopy (vecR0, 0, mtxB, 0);
111 DenseMtx_colCopy (vecR, 0, vecR0, 0);
112 Iter = 0;
113 Imv = 0;
114
115 fprintf(msgFile, "\n ZBiCGSTABR Initial norml: %6.2ef ", Init_norm ) ;
116 fprintf(msgFile, "\n ZBiCGSTABR Conveg. Control: %12.8f ", convergetol ) ;
117
118 /*
119
120 */
121 Rho[0] = 1.0;
122 Alpha[0] = 1.0;
123 Omega[0] = 1.0;
124 Rho[1] = 0.0;
125 Alpha[1] = 0.0;
126 Omega[1] = 0.0;
127
128 DenseMtx_zero(vecV) ;
129 DenseMtx_zero(vecP) ;
130 /*
131 ------------------------------
132
133 ------------------------------
134 */
135
136 MARKTIME(t1) ;
137
138 /*
139 -----------------
140 factor the matrix
141 -----------------
142 */
143 while ( ratio > convergetol && Iter <= itermax )
144 {
145 Iter++;
146 /* Rho_new = DenseMtx_dot(vecR0, vecR); */
147 DenseMtx_colDotProduct (vecR0, 0, vecR,0, Rho_new);
148 if ( zabs(Rho_new) == 0.0 || zabs(Omega)== 0.0 ){
149 fprintf(stderr, "\n breakdown in ZBiCGSTABR !! "
150 "\n Fatal error \n");
151 exit(-1) ;
152 }
153
154 /* Beta = (Rho_new / (Rho+Tiny)) * (Alpha / (Omega+Tiny)); */
155 zdiv(Rho_new, Rho, Beta);
156 zmul(Beta, Alpha, Rtmp);
157 zdiv(Rtmp, Omega, Beta);
158
159 Rho[0] = Rho_new[0];
160 Rho[1] = Rho_new[1];
161 /* DenseMtx_axpy(vecP, vecV, -Omega); */
162 /* DenseMtx_aypx(vecP, vecR, Beta); */
163 zsub(zero, Omega, Rtmp);
164 DenseMtx_colGenAxpy (one, vecP, 0, Rtmp, vecV, 0 );
165 DenseMtx_colGenAxpy (Beta, vecP, 0, one, vecR, 0 );
166 /* */
167 FrontMtx_solve(Precond, vecY, vecP, Precond->manager,
168 cpus, msglvl, msgFile) ;
169 /* */
170 /* InpMtx_nonsym_gmmm(mtxA, zero, vecV, one, vecY) ; */
171 switch ( symmetryflag ) {
172 case SPOOLES_SYMMETRIC :
173 InpMtx_sym_gmmm(mtxA, zero, vecV, one, vecY) ;
174 break ;
175 case SPOOLES_HERMITIAN :
176 InpMtx_herm_gmmm(mtxA, zero, vecV, one, vecY) ;
177 break ;
178 case SPOOLES_NONSYMMETRIC :
179 InpMtx_nonsym_gmmm(mtxA, zero, vecV, one, vecY) ;
180 break ;
181 default :
182 fprintf(msgFile, "\n BiCGSTABR Matrix type wrong");
183 fprintf(msgFile, "\n Fatal error");
184 goto end;
185 }
186 Imv++;
187
188 /* Alpha = Rho / (DenseMtx_dot(vecR0, vecV)+Tiny); */
189
190 DenseMtx_colDotProduct (vecR0, 0, vecV,0, Rtmp);
191 if ( zabs(Rtmp) == 0.0 ){
192 fprintf(stderr, "\n breakdown in ZBiCGSTABR !! "
193 "\n Fatal error \n");
194 exit(-1) ;
195 }
196 zdiv(Rho, Rtmp, Alpha);
197
198 /* */
199 /* DenseMtx_axpy(vecR, vecV, -Alpha); */
200 zsub(zero, Alpha, Rtmp);
201 DenseMtx_colGenAxpy (one, vecR, 0, Rtmp, vecV, 0 );
202
203 /* */
204 FrontMtx_solve(Precond, vecZ, vecR, Precond->manager,
205 cpus, msglvl, msgFile) ;
206
207 /* InpMtx_nonsym_gmmm(mtxA, zero, vecT, one, vecZ) ; */
208 switch ( symmetryflag ) {
209 case SPOOLES_SYMMETRIC :
210 InpMtx_sym_gmmm(mtxA, zero, vecT, one, vecZ) ;
211 break ;
212 case SPOOLES_HERMITIAN :
213 InpMtx_herm_gmmm(mtxA, zero, vecT, one, vecZ) ;
214 break ;
215 case SPOOLES_NONSYMMETRIC :
216 InpMtx_nonsym_gmmm(mtxA, zero, vecT, one, vecZ) ;
217 break ;
218 default :
219 fprintf(msgFile, "\n BiCGSTABR Matrix type wrong");
220 fprintf(msgFile, "\n Fatal error");
221 goto end;
222 }
223 Imv++;
224
225 /* Omega = DenseMtx_dot(vecT, vecR)/(DenseMtx_dot(vecT, vecT)+Tiny); */
226
227 DenseMtx_colDotProduct (vecT, 0, vecT,0, Ttmp);
228 if ( zabs(Ttmp) == 0.0 ){
229 fprintf(stderr, "\n breakdown in ZBiCGSTABR !! "
230 "\n Fatal error \n");
231 exit(-1) ;
232 };
233 DenseMtx_colDotProduct (vecT, 0, vecR, 0, Rtmp);
234 zdiv(Rtmp, Ttmp, Omega);
235
236
237 DenseMtx_colGenAxpy (one, vecX, 0, Alpha, vecY, 0);
238 DenseMtx_colGenAxpy (one, vecX, 0, Omega, vecZ, 0);
239
240 /* DenseMtx_axpy(vecR, vecT, -Omega); */
241 zsub(zero, Omega, Rtmp);
242 DenseMtx_colGenAxpy (one, vecR, 0, Rtmp, vecT, 0);
243
244 Res_norm = DenseMtx_twoNormOfColumn (vecR, 0);
245 ratio = Res_norm/Init_norm;
246 fprintf(msgFile, "\n\n At iteration %d"
247 " the convergence ratio is %12.4e"
248 "\n Residual norm is %6.2e",
249 Imv, ratio, Res_norm) ;
250 }
251 /* End of while loop */
252 MARKTIME(t2) ;
253 fprintf(msgFile, "\n CPU : Converges in time: %8.3f ", t2 - t1) ;
254 fprintf(msgFile, "\n # iterations = %d", Imv) ;
255
256 fprintf(msgFile, "\n\n after ZBICGSTABR") ;
257 /* DenseMtx_copy(mtxX, vecX); */
258 DenseMtx_colCopy (mtxX, 0, vecX, 0);
259 /*
260
261 ------------------------
262 free the working storage
263 ------------------------
264 */
265 end:
266 DenseMtx_free(vecP) ;
267 DenseMtx_free(vecR) ;
268 DenseMtx_free(vecR0) ;
269 DenseMtx_free(vecT) ;
270 DenseMtx_free(vecV) ;
271 DenseMtx_free(vecW) ;
272 DenseMtx_free(vecX) ;
273 DenseMtx_free(vecY) ;
274 DenseMtx_free(vecZ) ;
275
276 fprintf(msgFile, "\n") ;
277
278 return(1) ; }
279
280 /*--------------------------------------------------------------------*/
281