1 /* pcgl.c */
2
3 #include "../Iter.h"
4
5 /*--------------------------------------------------------------------*/
6 /*
7 ---------------------------------------------------------------------
8 purpose -- to solve a symmetric nonnegative matrix equation
9
10 Ax=b
11 using left preconditioned conjugate gradient method
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
27 created -- Oct. 27, 1998 Wei-Pai Tang
28 ---------------------------------------------------------------------
29 */
30
31 int
pcgl(int n_matrixSize,int type,int symmetryflag,InpMtx * mtxA,FrontMtx * Precond,DenseMtx * mtxX,DenseMtx * mtxB,int itermax,double convergetol,int msglvl,FILE * msgFile)32 pcgl (
33 int n_matrixSize,
34 int type,
35 int symmetryflag,
36 InpMtx *mtxA,
37 FrontMtx *Precond,
38 DenseMtx *mtxX,
39 DenseMtx *mtxB,
40 int itermax,
41 double convergetol,
42 int msglvl,
43 FILE *msgFile
44 )
45 {
46 Chv *chv, *rootchv ;
47 ChvManager *chvmanager ;
48 DenseMtx *mtxZ ;
49 DenseMtx *vecP, *vecR, *vecQ ;
50 DenseMtx *vecX, *vecZ ;
51 double Alpha, Beta, Rho, Rho0, Init_norm, ratio, Res_norm, Rtmp ;
52 double t1, t2, cpus[9] ;
53 double one[2] = {1.0, 0.0}, zero[2] = {0.0, 0.0} ;
54 double Tiny = 0.1e-28;
55 int Iter, neqns;
56 int stats[6] ;
57
58 if (symmetryflag != SPOOLES_SYMMETRIC){
59 fprintf(msgFile, "\n\n Fatal Error, \n"
60 " Matrix is not symmetric in PCGL !!") ;
61 return (-1);
62 };
63
64 neqns = n_matrixSize;
65
66 /*
67 --------------------
68 init the vectors in bicgstab
69 --------------------
70 */
71 vecP = DenseMtx_new() ;
72 DenseMtx_init(vecP, type, 0, 0, neqns, 1, 1, neqns) ;
73
74 vecR = DenseMtx_new() ;
75 DenseMtx_init(vecR, type, 0, 0, neqns, 1, 1, neqns) ;
76
77 vecX = DenseMtx_new() ;
78 DenseMtx_init(vecX, type, 0, 0, neqns, 1, 1, neqns) ;
79
80 vecQ = DenseMtx_new() ;
81 DenseMtx_init(vecQ, type, 0, 0, neqns, 1, 1, neqns) ;
82
83 vecZ = DenseMtx_new() ;
84 DenseMtx_init(vecZ, type, 0, 0, neqns, 1, 1, neqns) ;
85
86
87 /*
88 --------------------------
89 Initialize the iterations
90 --------------------------
91 */
92 FrontMtx_solve(Precond, vecR, mtxB, Precond->manager,
93 cpus, msglvl, msgFile) ;
94
95 Init_norm = DenseMtx_twoNormOfColumn(vecR,0);
96 if ( Init_norm == 0.0 ){
97 Init_norm = 1.0; };
98 ratio = 1.0;
99 DenseMtx_zero(vecX) ;
100
101 /* DenseMtx_copy(vecR, mtxB); */
102
103
104 MARKTIME(t1) ;
105 Iter = 0;
106
107 /*
108 ------------------------------
109 Main Loop of the iterations
110 ------------------------------
111 */
112
113 while ( ratio > convergetol && Iter <= itermax )
114 {
115 Iter++;
116 /* */
117 DenseMtx_colDotProduct(vecR, 0, vecR, 0, &Rho);
118 /* */
119 if ( Iter == 1 ) {
120 DenseMtx_colCopy(vecP, 0, vecR, 0);
121 } else {
122 Beta = Rho /(Rho0 + Tiny);
123 DenseMtx_colGenAxpy(&Beta, vecP, 0, one, vecR, 0);
124 };
125
126 InpMtx_sym_gmmm(mtxA, zero, vecZ, one, vecP) ;
127 FrontMtx_solve(Precond, vecQ, vecZ, Precond->manager,
128 cpus, msglvl, msgFile) ;
129
130 DenseMtx_colDotProduct(vecP, 0, vecQ, 0, &Rtmp);
131 Alpha = Rho/(Rtmp+Tiny);
132 DenseMtx_colGenAxpy(one, vecX, 0, &Alpha, vecP, 0);
133 Rtmp=-Alpha;
134 DenseMtx_colGenAxpy(one, vecR, 0, &Rtmp, vecQ, 0);
135 Rho0 = Rho;
136
137 /* */
138 Res_norm = DenseMtx_twoNormOfColumn(vecR,0);
139 ratio = Res_norm/Init_norm;
140 fprintf(msgFile, "\n\n At iteration %d"
141 " the convergence ratio is %12.4e",
142 Iter, ratio) ;
143 }
144 /* End of while loop */
145 MARKTIME(t2) ;
146 fprintf(msgFile, "\n CPU : Converges in time: %8.3f ", t2 - t1) ;
147 fprintf(msgFile, "\n # iterations = %d", Iter) ;
148
149 fprintf(msgFile, "\n\n after PCGL") ;
150 DenseMtx_colCopy(mtxX, 0, vecX, 0);
151
152 /*
153
154 ------------------------
155 free the working storage
156 ------------------------
157 */
158
159
160 end:
161
162 DenseMtx_free(vecP) ;
163 DenseMtx_free(vecR) ;
164 DenseMtx_free(vecX) ;
165 DenseMtx_free(vecQ) ;
166 DenseMtx_free(vecZ) ;
167
168
169 fprintf(msgFile, "\n") ;
170
171 return(1) ; }
172
173 /*--------------------------------------------------------------------*/
174