1 /* pcgr.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 right 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 return value -- final residual norm
27
28 created -- Oct. 27, 1998 Wei-Pai Tang
29 ---------------------------------------------------------------------
30 */
31
32 int
pcgr(int n_matrixSize,int type,int symmetryflag,InpMtx * mtxA,FrontMtx * Precond,DenseMtx * mtxX,DenseMtx * mtxB,int itermax,double convergetol,int msglvl,FILE * msgFile)33 pcgr (
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, *vecQ ;
51 DenseMtx *vecX, *vecZ ;
52 double Alpha, Beta, Rho, Rho0, Init_norm, ratio, Res_norm, Rtmp ;
53 double t1, t2, cpus[9] ;
54 double one[2] = {1.0, 0.0} ;
55 double Tiny = 0.1e-28;
56 int Iter, neqns;
57 int stats[6] ;
58
59 if (symmetryflag != SPOOLES_SYMMETRIC){
60 fprintf(msgFile, "\n\n Fatal Error, \n"
61 " Matrix is not symmetric in PCGR !!") ;
62 return(-1);
63 };
64
65 neqns = n_matrixSize;
66
67 /*
68 --------------------
69 init the vectors in bicgstab
70 --------------------
71 */
72 vecP = DenseMtx_new() ;
73 DenseMtx_init(vecP, type, 0, 0, neqns, 1, 1, neqns) ;
74
75 vecR = DenseMtx_new() ;
76 DenseMtx_init(vecR, type, 0, 0, neqns, 1, 1, neqns) ;
77
78 vecX = DenseMtx_new() ;
79 DenseMtx_init(vecX, type, 0, 0, neqns, 1, 1, neqns) ;
80
81 vecQ = DenseMtx_new() ;
82 DenseMtx_init(vecQ, type, 0, 0, neqns, 1, 1, neqns) ;
83
84 vecZ = DenseMtx_new() ;
85 DenseMtx_init(vecZ, type, 0, 0, neqns, 1, 1, neqns) ;
86
87
88 /*
89 --------------------------
90 Initialize the iterations
91 --------------------------
92 */
93 Init_norm = DenseMtx_twoNormOfColumn(mtxB,0);
94 if ( Init_norm == 0.0 ){
95 Init_norm = 1.0; };
96 ratio = 1.0;
97 DenseMtx_zero(vecX) ;
98
99 DenseMtx_colCopy(vecR, 0, mtxB, 0);
100
101
102 MARKTIME(t1) ;
103 Iter = 0;
104
105 /*
106 ------------------------------
107 Main Loop of the iterations
108 ------------------------------
109 */
110
111 while ( ratio > convergetol && Iter <= itermax )
112 {
113 Iter++;
114 /* */
115 FrontMtx_solve(Precond, vecZ, vecR, Precond->manager,
116 cpus, msglvl, msgFile) ;
117 DenseMtx_colDotProduct(vecR, 0, vecZ, 0, &Rho);
118 /* */
119 if ( Iter == 1 ) {
120 DenseMtx_colCopy(vecP, 0, vecZ, 0);
121 } else {
122 Beta = Rho /(Rho0 + Tiny);
123 DenseMtx_colGenAxpy(&Beta, vecP, 0, one, vecZ, 0);
124 };
125
126 DenseMtx_zero(vecQ);
127 InpMtx_sym_mmm(mtxA, vecQ, one, vecP) ;
128 DenseMtx_colDotProduct(vecP, 0, vecQ, 0, &Rtmp);
129 Alpha = Rho/(Rtmp+Tiny);
130 DenseMtx_colGenAxpy(one, vecX, 0, &Alpha, vecP, 0);
131 Rtmp=-Alpha;
132 DenseMtx_colGenAxpy(one, vecR, 0, &Rtmp, vecQ, 0);
133 Rho0 = Rho;
134
135 /* */
136 Res_norm = DenseMtx_twoNormOfColumn(vecR,0);
137 ratio = Res_norm/Init_norm;
138 fprintf(msgFile, "\n\n At iteration %d"
139 " the convergence ratio is %12.4e",
140 Iter, ratio) ;
141 }
142 /* End of while loop */
143 MARKTIME(t2) ;
144 fprintf(msgFile, "\n CPU : Converges in time: %8.3f ", t2 - t1) ;
145 fprintf(msgFile, "\n # iterations = %d", Iter) ;
146
147 fprintf(msgFile, "\n\n after PCGR") ;
148 DenseMtx_colCopy(mtxX, 0, vecX, 0);
149
150 /*
151
152 ------------------------
153 free the working storage
154 ------------------------
155 */
156
157 DenseMtx_free(vecP) ;
158 DenseMtx_free(vecR) ;
159 DenseMtx_free(vecX) ;
160 DenseMtx_free(vecQ) ;
161 DenseMtx_free(vecZ) ;
162
163
164 fprintf(msgFile, "\n") ;
165
166 return(1) ; }
167
168 /*--------------------------------------------------------------------*/
169