1 #include "BSprivate.h"
2
3 /*@ BSpar_sym_solve - Solve a symmetric positive definite system of equations
4 using conjugate gradients preconditioned by one of several
5 preconditioners. The rhs can be a block of vectors. The
6 user should not call this function directly, but BSpar_solve().
7
8 Input Parameters:
9 . BS - the number of vectors in the RHS
10 . A - a sparse matrix
11 . fact_A - the incomplete factored version of A, if any
12 . comm_A - the communication structure for A
13 . in_rhs - the contiguous block of vectors forming the rhs
14 . pre_option - the preconditioner to use
15 $ PRE_ICC: incomplete Cholesky factorization
16 $ PRE_ILU: incomplete LU factorization
17 $ PRE_SSOR: Successive over relaxation
18 $ PRE_BJACOBI: Block Jacobi
19 . err_tol - the tolerance to which to solve the problem
20 stop if the estimated norm of the residual divided by
21 the norm of the rhs is less than err_tol
22 . max_iter - the maximum number of iterations to take
23 . residual - the final computed residual
24 . guess - if TRUE, then initialize out_x to 0, otherwise the program
25 assumes that out_x contains an initial guess
26 . procinfo - the usual processor stuff
27
28 Output Parameters:
29 . out_x - the contiguous block of vectors containing the solution
30
31 Returns:
32 The number of iterations or a negative number indicating the number
33 of iterations prior to finding that the matrix (or preconditioner)
34 is not positive definite.
35
36 Notes:
37 The preconditioners must be computed prior to calling BSpar_isolve.
38 For more information on the preconditioners, see the manual.
39 @*/
BSpar_sym_solve(int BS,BSpar_mat * A,BSpar_mat * fact_A,BScomm * comm_A,FLOAT * in_rhs,FLOAT * out_x,int pre_option,FLOAT err_tol,int max_iter,FLOAT * residual,int guess,BSprocinfo * procinfo)40 int BSpar_sym_solve(int BS, BSpar_mat *A, BSpar_mat *fact_A, BScomm *comm_A,
41 FLOAT *in_rhs, FLOAT *out_x, int pre_option, FLOAT err_tol, int max_iter,
42 FLOAT *residual, int guess, BSprocinfo *procinfo)
43 {
44 int i, j;
45 int cur_step, cur_phase;
46 int done;
47 FLOAT *resid, *z, *p;
48 FLOAT *cg_beta, *cg_alpha;
49 FLOAT *bnorm;
50 FLOAT *x;
51 FLOAT *rhs;
52 FLOAT tval;
53 FLOAT *t_x, *t_rhs;
54 int n;
55
56 if(!A->symmetric) {
57 MY_SETERRCN(PAR_SOLVE_ERROR,"Trying to solve nonsymmetric system with CG\n");
58 }
59
60 n = A->num_rows;
61
62 /* reorder the rhs */
63 MY_MALLOCN(rhs,(FLOAT *),sizeof(FLOAT)*n*BS,1);
64 for (i=0;i<BS;i++) {
65 BSperm_dvec(&(in_rhs[i*n]),&(rhs[i*n]),A->perm); CHKERRN(0);
66 }
67
68 /* allocate space for x */
69 MY_MALLOCN(x,(FLOAT *),sizeof(FLOAT)*n*BS,2);
70
71 /* allocate space for cg vectors */
72 MY_MALLOCN(resid,(FLOAT *),sizeof(FLOAT)*n*BS,3);
73 MY_MALLOCN(p,(FLOAT *),sizeof(FLOAT)*n*BS,4);
74 MY_MALLOCN(z,(FLOAT *),sizeof(FLOAT)*n*BS,5);
75 MY_MALLOCN(cg_alpha,(FLOAT *),sizeof(FLOAT)*BS,6);
76 MY_MALLOCN(cg_beta,(FLOAT *),sizeof(FLOAT)*BS,7);
77 MY_MALLOCN(bnorm,(FLOAT *),sizeof(FLOAT)*BS,8);
78
79 /* form the initial guess */
80 if (guess) {
81 for (i=0;i<n*BS;i++) {
82 x[i] = 0.0;
83 }
84 } else {
85 for (i=0;i<BS;i++) {
86 BSperm_dvec(&(out_x[i*n]),&(x[i*n]),A->perm); CHKERRN(0);
87 }
88 }
89
90 /* scale the rhs and x */
91 if(A->scale_diag!=NULL) {
92 for (i=0;i<BS;i++) {
93 t_rhs = &(rhs[i*n]);
94 t_x = &(x[i*n]);
95 for (j=0;j<n;j++) {
96 tval = sqrt(A->scale_diag[j]);
97 t_rhs[j] /= tval;
98 t_x[j] *= tval;
99 }
100 }
101 }
102
103 /* get the norm of B */
104 BSpar_bip(n,rhs,rhs,BS,bnorm,procinfo); CHKERRN(0);
105 done = TRUE;
106 for (i=0;i<BS;i++) {
107 bnorm[i] = sqrt(bnorm[i]);
108 if (bnorm[i] != 0.0) done = FALSE;
109 }
110 cur_step = 0;
111 cur_phase = 0;
112 while ((!done) && (cur_step < max_iter)) {
113 switch(BSpar_bcg(n,rhs,x,resid,z,p,cg_beta,cg_alpha,
114 &cur_step,&cur_phase,BS,procinfo)) {
115 case CG_MATVECR: {
116 BStri_mult(A,comm_A,NULL,NULL,x,resid,NULL,NULL,0.0,BS,
117 procinfo); CHKERRN(0);
118 break;
119 }
120 case CG_MATVECZ: {
121 BStri_mult(A,comm_A,NULL,NULL,p,z,NULL,NULL,0.0,BS,
122 procinfo); CHKERRN(0);
123 break;
124 }
125 case CG_MSOLVE: {
126 BSpar_bip(n,resid,resid,BS,residual,procinfo); CHKERRN(0);
127 for (i=0;i<BS;i++) {
128 if (bnorm[i] != 0.0) {
129 residual[i] = sqrt(residual[i])/bnorm[i];
130 } else {
131 residual[i] = sqrt(residual[i]);
132 }
133 }
134 done = TRUE;
135 for (i=0;i<BS;i++) {
136 if (residual[i] >= err_tol) done = FALSE;
137 }
138 if (!done) {
139 for (i=0;i<n*BS;i++) {
140 z[i] = resid[i];
141 }
142 BStri_solve(A,fact_A,comm_A,z,pre_option,BS,procinfo);
143 CHKERRN(0);
144 }
145 break;
146 }
147 default: {
148 return(-cur_step);
149 }
150 }
151 CHKERRN(0);
152 }
153
154 BStri_mult(A,comm_A,NULL,NULL,x,resid,NULL,NULL,0.0,BS,procinfo); CHKERRN(0);
155 for (i=0;i<n*BS;i++) {
156 resid[i] = rhs[i]-resid[i];
157 }
158 BSpar_bip(n,resid,resid,BS,residual,procinfo); CHKERRN(0);
159 for (i=0;i<BS;i++) {
160 if (bnorm[i] != 0.0) {
161 residual[i] = sqrt(residual[i])/bnorm[i];
162 } else {
163 residual[i] = sqrt(residual[i]);
164 }
165 }
166
167 MY_FREE(z);
168 MY_FREE(p);
169 MY_FREE(resid);
170 MY_FREE(cg_alpha);
171 MY_FREE(cg_beta);
172 MY_FREE(bnorm);
173
174 /* Rescale X */
175 if(A->scale_diag!=NULL) {
176 for (i=0;i<BS;i++) {
177 t_x = &(x[i*n]);
178 for (j=0;j<n;j++) {
179 t_x[j] /= sqrt(A->scale_diag[j]);
180 }
181 }
182 }
183
184 /* reorder the solution vector */
185 for (i=0;i<BS;i++) {
186 BSiperm_dvec(&(x[i*n]),&(out_x[i*n]),A->perm); CHKERRN(0);
187 }
188
189 MY_FREE(rhs);
190 MY_FREE(x);
191
192 /* return the number of iterations */
193 return(cur_step);
194 }
195