1 #include "BSprivate.h"
2
3 /*@ BSpar_isolve - Solve a symmetric indefinite system of equations
4 using symmlq preconditioned by one of several
5 preconditioners.
6
7 Input Parameters:
8 . A - a sparse matrix
9 . fact_A - the incomplete factored version of A, if any (NULL if not exist)
10 . comm_A - the communication structure for A
11 . B - a sparse matrix
12 . comm_B - the communication structure for B
13 . in_rhs - the rhs
14 . shift - the shift to multiply B by
15 . pre_option - the preconditioner to use
16 $ PRE_ICC: incomplete Cholesky factorization of A
17 $ PRE_ILU: incomplete LU factorization of A
18 $ PRE_SSOR: Successive over relaxation (using just A)
19 $ PRE_BJACOBI: Block Jacobi (using just A)
20 . residual - the final computed residual
21 . procinfo - the usual processor stuff
22
23 Output Parameters:
24 . out_x - the solution vector
25
26 Returns:
27 The number of iterations.
28
29 Notes:
30 The system solved is (A-shift*B)out_x = in_rhs.
31 The preconditioners must be computed prior to calling BSpar_isolve.
32 For more information on the preconditioners, see the manual.
33
34 $ The following are now specified in the context:
35 $ err_tol - the tolerance to which to solve the problem
36 $ stop if the estimated norm of the residual divided by
37 $ the norm of the rhs is less than err_tol
38 $ max_iter - the maximum number of iterations to take
39 $ guess - if TRUE, then initialize out_x to 0, otherwise
40 $ the program assumes that out_x contains an initial
41 $ guess
42
43 @*/
BSpar_isolve(BSpar_mat * A,BSpar_mat * fact_A,BScomm * comm_A,BSpar_mat * B,BScomm * comm_B,FLOAT * in_rhs,FLOAT * out_x,FLOAT shift,FLOAT * residual,BSprocinfo * procinfo)44 int BSpar_isolve(BSpar_mat *A, BSpar_mat *fact_A, BScomm *comm_A,
45 BSpar_mat *B, BScomm *comm_B, FLOAT *in_rhs, FLOAT *out_x,
46 FLOAT shift, FLOAT *residual, BSprocinfo *procinfo)
47 {
48 /* temporary space to hold an incoming vector */
49 int i;
50 int cur_step, cur_phase;
51 int done, conv;
52 FLOAT *resid, *z, *v, *w, *r2;
53 FLOAT cg_beta, cg_alpha;
54 FLOAT obeta;
55 FLOAT gbar, dbar;
56 FLOAT rhs1, rhs2;
57 FLOAT snprod;
58 FLOAT beta1;
59 FLOAT bstep;
60 FLOAT ynorm2;
61 FLOAT tnorm;
62 FLOAT bestnm;
63 int cgpt;
64 FLOAT bnorm;
65 FLOAT *x;
66 FLOAT *rhs;
67 FLOAT *temp1, *temp2;
68 FLOAT tval;
69 int pre_option, max_iter, guess;
70 FLOAT err_tol;
71
72 if ((procinfo->print) && (PSISROOT(procinfo))) {
73 BSctx_print(procinfo);
74 }
75
76 pre_option = procinfo->preconditioner;
77 err_tol = procinfo->residual_tol;
78 max_iter = procinfo->max_iterations;
79 guess = procinfo->guess;
80
81 /* reorder the rhs */
82 MY_MALLOCN(rhs,(FLOAT *),sizeof(FLOAT)*A->num_rows,1);
83 BSperm_dvec(in_rhs,rhs,A->perm); CHKERRN(0);
84
85 /* allocate space for x */
86 MY_MALLOCN(x,(FLOAT *),sizeof(FLOAT)*A->num_rows,1);
87
88 /* allocate space for symmlq vectors */
89 MY_MALLOCN(resid,(FLOAT *),sizeof(FLOAT)*A->num_rows,1);
90 MY_MALLOCN(v,(FLOAT *),sizeof(FLOAT)*A->num_rows,1);
91 MY_MALLOCN(r2,(FLOAT *),sizeof(FLOAT)*A->num_rows,1);
92 MY_MALLOCN(w,(FLOAT *),sizeof(FLOAT)*A->num_rows,1);
93 MY_MALLOCN(z,(FLOAT *),sizeof(FLOAT)*A->num_rows,1);
94 MY_MALLOCN(temp1,(FLOAT *),sizeof(FLOAT)*A->num_rows,1);
95 MY_MALLOCN(temp2,(FLOAT *),sizeof(FLOAT)*A->num_rows,1);
96
97 /* form the initial guess */
98 if (guess) {
99 for (i=0;i<A->num_rows;i++) {
100 x[i] = 0.0;
101 }
102 } else {
103 BSperm_dvec(out_x,x,A->perm);
104 }
105
106 /* scale the rhs */
107 for (i=0;i<A->num_rows;i++) {
108 tval = sqrt(A->scale_diag[i]);
109 rhs[i] /= tval;
110 x[i] *= tval;
111 }
112
113 /* get the norm of B */
114 bnorm = BSpar_ip(A->num_rows,rhs,rhs,procinfo); CHKERRN(0);
115 bnorm = sqrt(bnorm);
116 done = FALSE;
117 cur_step = 0;
118 cur_phase = 0;
119 if (bnorm == 0.0) done = TRUE;
120 conv = FALSE;
121 while ((!done) && (cur_step < max_iter)) {
122 switch(BSpar_symmlq(A->num_rows,rhs,x,resid,z,v,w,r2,&cg_beta,&cg_alpha,
123 &obeta,&gbar,&dbar,&rhs1,&rhs2,
124 &bstep,&snprod,&beta1,&tnorm,&ynorm2,&bestnm,&cgpt,conv,
125 &cur_step,&cur_phase,procinfo)) {
126 case SYMM_MATVECXR: {
127 BStri_mult(A,comm_A,B,comm_B,x,resid,temp1,temp2,shift,1,
128 procinfo); CHKERRN(0);
129 break;
130 }
131 case SYMM_MATVECVZ: {
132 /* check for convergence first */
133 if ((cur_step > 0) && (bestnm/bnorm < err_tol)) conv = TRUE;
134 if (!conv) {
135 BStri_mult(A,comm_A,B,comm_B,v,z,temp1,temp2,shift,1,
136 procinfo); CHKERRN(0);
137 }
138 break;
139 }
140 case SYMM_MSOLVE1: {
141 for (i=0;i<A->num_rows;i++) {
142 z[i] = resid[i];
143 }
144 BStri_solve(A,fact_A,comm_A,z,pre_option,1,procinfo); CHKERRN(0);
145 break;
146 }
147 case SYMM_MSOLVE2: {
148 for (i=0;i<A->num_rows;i++) {
149 z[i] = r2[i];
150 }
151 BStri_solve(A,fact_A,comm_A,z,pre_option,1,procinfo); CHKERRN(0);
152 break;
153 }
154 case SYMM_MSOLVEB: {
155 for (i=0;i<A->num_rows;i++) {
156 z[i] = rhs[i];
157 }
158 BStri_solve(A,fact_A,comm_A,z,pre_option,1,procinfo); CHKERRN(0);
159 break;
160 }
161 case SYMM_DONE: {
162 done = TRUE;
163 break;
164 }
165 default: {
166 return(-cur_step);
167 }
168 }
169 CHKERRN(0);
170 }
171
172 BStri_mult(A,comm_A,B,comm_B,x,resid,temp1,temp2,shift,1,procinfo);
173 CHKERRN(0);
174 for (i=0;i<A->num_rows;i++) {
175 resid[i] = rhs[i] - resid[i];
176 }
177 (*residual) = BSpar_ip(A->num_rows,resid,resid,procinfo); CHKERRN(0);
178 if (bnorm != 0.0) {
179 (*residual) = sqrt((*residual))/bnorm;
180 } else {
181 (*residual) = sqrt((*residual));
182 }
183
184 MY_FREE(temp1);
185 MY_FREE(temp2);
186 MY_FREE(z);
187 MY_FREE(v);
188 MY_FREE(w);
189 MY_FREE(r2);
190 MY_FREE(resid);
191
192 /* Rescale X */
193 for (i=0;i<A->num_rows;i++) {
194 x[i] /= sqrt(fabs(A->scale_diag[i]));
195 }
196
197 /* reorder the solution vector */
198 BSiperm_dvec(x,out_x,A->perm); CHKERRN(0);
199 MY_FREE(rhs);
200 MY_FREE(x);
201
202 /* return the number of iterations */
203 return(cur_step);
204 }
205