1/* $Id$ $Revision$ */
2/* vim:set shiftwidth=4 ts=8: */
3
4/*************************************************************************
5 * Copyright (c) 2011 AT&T Intellectual Property
6 * All rights reserved. This program and the accompanying materials
7 * are made available under the terms of the Eclipse Public License v1.0
8 * which accompanies this distribution, and is available at
9 * http://www.eclipse.org/legal/epl-v10.html
10 *
11 * Contributors: See CVS logs. Details at http://www.graphviz.org/
12 *************************************************************************/
13
14#include <assert.h>
15#include <string.h>
16#include "sparse_solve.h"
17#include "sfdpinternal.h"
18#include "memory.h"
19#include "logic.h"
20#include "math.h"
21#include "arith.h"
22#include "types.h"
23#include "globals.h"
24
25/* #define DEBUG_PRINT */
26
27struct uniform_stress_matmul_data{
28 real alpha;
29 SparseMatrix A;
30};
31
32
33void Operator_uniform_stress_matmul_delete(Operator o){
34 FREE(o->data);
35}
36
37real *Operator_uniform_stress_matmul_apply(Operator o, real *x, real *y){
38 struct uniform_stress_matmul_data *d = (struct uniform_stress_matmul_data*) (o->data);
39 SparseMatrix A = d->A;
40 real alpha = d->alpha;
41 real xsum = 0.;
42 int m = A->m, i;
43
44 SparseMatrix_multiply_vector(A, x, &y, FALSE);
45
46 /* alpha*V*x */
47 for (i = 0; i < m; i++) xsum += x[i];
48
49 for (i = 0; i < m; i++) y[i] += alpha*(m*x[i] - xsum);
50
51 return y;
52}
53
54
55
56Operator Operator_uniform_stress_matmul(SparseMatrix A, real alpha){
57 Operator o;
58 struct uniform_stress_matmul_data *d;
59
60 o = MALLOC(sizeof(struct Operator_struct));
61 o->data = d = MALLOC(sizeof(struct uniform_stress_matmul_data));
62 d->alpha = alpha;
63 d->A = A;
64 o->Operator_apply = Operator_uniform_stress_matmul_apply;
65 return o;
66}
67
68
69real *Operator_matmul_apply(Operator o, real *x, real *y){
70 SparseMatrix A = (SparseMatrix) o->data;
71 SparseMatrix_multiply_vector(A, x, &y, FALSE);
72 return y;
73}
74
75Operator Operator_matmul_new(SparseMatrix A){
76 Operator o;
77
78 o = GNEW(struct Operator_struct);
79 o->data = (void*) A;
80 o->Operator_apply = Operator_matmul_apply;
81 return o;
82}
83
84
85void Operator_matmul_delete(Operator o){
86 if (o) FREE(o);
87}
88
89
90real* Operator_diag_precon_apply(Operator o, real *x, real *y){
91 int i, m;
92 real *diag = (real*) o->data;
93 m = (int) diag[0];
94 diag++;
95 for (i = 0; i < m; i++) y[i] = x[i]*diag[i];
96 return y;
97}
98
99
100Operator Operator_uniform_stress_diag_precon_new(SparseMatrix A, real alpha){
101 Operator o;
102 real *diag;
103 int i, j, m = A->m, *ia = A->ia, *ja = A->ja;
104 real *a = (real*) A->a;
105
106 assert(A->type == MATRIX_TYPE_REAL);
107
108 assert(a);
109
110 o = MALLOC(sizeof(struct Operator_struct));
111 o->data = MALLOC(sizeof(real)*(m + 1));
112 diag = (real*) o->data;
113
114 diag[0] = m;
115 diag++;
116 for (i = 0; i < m; i++){
117 diag[i] = 1./(m-1);
118 for (j = ia[i]; j < ia[i+1]; j++){
119 if (i == ja[j] && ABS(a[j]) > 0) diag[i] = 1./((m-1)*alpha+a[j]);
120 }
121 }
122
123 o->Operator_apply = Operator_diag_precon_apply;
124
125 return o;
126}
127
128
129Operator Operator_diag_precon_new(SparseMatrix A){
130 Operator o;
131 real *diag;
132 int i, j, m = A->m, *ia = A->ia, *ja = A->ja;
133 real *a = (real*) A->a;
134
135 assert(A->type == MATRIX_TYPE_REAL);
136
137 assert(a);
138
139 o = N_GNEW(1,struct Operator_struct);
140 o->data = N_GNEW((A->m + 1),real);
141 diag = (real*) o->data;
142
143 diag[0] = m;
144 diag++;
145 for (i = 0; i < m; i++){
146 diag[i] = 1.;
147 for (j = ia[i]; j < ia[i+1]; j++){
148 if (i == ja[j] && ABS(a[j]) > 0) diag[i] = 1./a[j];
149 }
150 }
151
152 o->Operator_apply = Operator_diag_precon_apply;
153
154 return o;
155}
156
157void Operator_diag_precon_delete(Operator o){
158 if (o->data) FREE(o->data);
159 if (o) FREE(o);
160}
161
162static real conjugate_gradient(Operator A, Operator precon, int n, real *x, real *rhs, real tol, int maxit, int *flag){
163 real *z, *r, *p, *q, res = 10*tol, alpha;
164 real rho = 1.0e20, rho_old = 1, res0, beta;
165 real* (*Ax)(Operator o, real *in, real *out) = A->Operator_apply;
166 real* (*Minvx)(Operator o, real *in, real *out) = precon->Operator_apply;
167 int iter = 0;
168
169 z = N_GNEW(n,real);
170 r = N_GNEW(n,real);
171 p = N_GNEW(n,real);
172 q = N_GNEW(n,real);
173
174 r = Ax(A, x, r);
175 r = vector_subtract_to(n, rhs, r);
176
177 res0 = res = sqrt(vector_product(n, r, r))/n;
178#ifdef DEBUG_PRINT
179 if (Verbose){
180 fprintf(stderr, "on entry, cg iter = %d of %d, residual = %g, tol = %g\n", iter, maxit, res, tol);
181 }
182#endif
183
184 while ((iter++) < maxit && res > tol*res0){
185 z = Minvx(precon, r, z);
186 rho = vector_product(n, r, z);
187
188 if (iter > 1){
189 beta = rho/rho_old;
190 p = vector_saxpy(n, z, p, beta);
191 } else {
192 MEMCPY(p, z, sizeof(real)*n);
193 }
194
195 q = Ax(A, p, q);
196
197 alpha = rho/vector_product(n, p, q);
198
199 x = vector_saxpy2(n, x, p, alpha);
200 r = vector_saxpy2(n, r, q, -alpha);
201
202 res = sqrt(vector_product(n, r, r))/n;
203
204#ifdef DEBUG_PRINT
205 if (Verbose && 0){
206 fprintf(stderr, " cg iter = %d, residual = %g, relative res = %g\n", iter, res, res/res0);
207 }
208#endif
209
210
211
212 rho_old = rho;
213 }
214 FREE(z); FREE(r); FREE(p); FREE(q);
215#ifdef DEBUG
216 _statistics[0] += iter - 1;
217#endif
218
219#ifdef DEBUG_PRINT
220 if (Verbose){
221 fprintf(stderr, " cg iter = %d, residual = %g, relative res = %g\n", iter, res, res/res0);
222 }
223#endif
224 return res;
225}
226
227real cg(Operator Ax, Operator precond, int n, int dim, real *x0, real *rhs, real tol, int maxit, int *flag){
228 real *x, *b, res = 0;
229 int k, i;
230 x = N_GNEW(n, real);
231 b = N_GNEW(n, real);
232 for (k = 0; k < dim; k++){
233 for (i = 0; i < n; i++) {
234 x[i] = x0[i*dim+k];
235 b[i] = rhs[i*dim+k];
236 }
237
238 res += conjugate_gradient(Ax, precond, n, x, b, tol, maxit, flag);
239 for (i = 0; i < n; i++) {
240 rhs[i*dim+k] = x[i];
241 }
242 }
243 FREE(x);
244 FREE(b);
245 return res;
246
247}
248
249real* jacobi(SparseMatrix A, int dim, real *x0, real *rhs, int maxit, int *flag){
250 /* maxit iteration of jacobi */
251 real *x, *y, *b, sum, diag, *a;
252 int k, i, j, n = A->n, *ia, *ja, iter;
253 x = MALLOC(sizeof(real)*n);
254 y = MALLOC(sizeof(real)*n);
255 b = MALLOC(sizeof(real)*n);
256 assert(A->type = MATRIX_TYPE_REAL);
257 ia = A->ia; ja = A->ja; a = (real*) A->a;
258
259 for (k = 0; k < dim; k++){
260 for (i = 0; i < n; i++) {
261 x[i] = x0[i*dim+k];
262 b[i] = rhs[i*dim+k];
263 }
264
265 for (iter = 0; iter < maxit; iter++){
266 for (i = 0; i < n; i++){
267 sum = 0; diag = 0;
268 for (j = ia[i]; j < ia[i+1]; j++){
269 if (ja[j] != i){
270 sum += a[j]*x[ja[j]];
271 } else {
272 diag = a[j];
273 }
274 }
275 if (sum == 0) fprintf(stderr,"neighb=%d\n",ia[i+1]-ia[i]);
276 assert(diag != 0);
277 y[i] = (b[i] - sum)/diag;
278
279 }
280 MEMCPY(x, y, sizeof(real)*n);
281 }
282
283 for (i = 0; i < n; i++) {
284 rhs[i*dim+k] = x[i];
285 }
286 }
287
288
289 FREE(x);
290 FREE(y);
291 FREE(b);
292 return rhs;
293
294}
295
296real SparseMatrix_solve(SparseMatrix A, int dim, real *x0, real *rhs, real tol, int maxit, int method, int *flag){
297 Operator Ax, precond;
298 int n = A->m;
299 real res = 0;
300 *flag = 0;
301
302 switch (method){
303 case SOLVE_METHOD_CG:
304 Ax = Operator_matmul_new(A);
305 precond = Operator_diag_precon_new(A);
306 res = cg(Ax, precond, n, dim, x0, rhs, tol, maxit, flag);
307 Operator_matmul_delete(Ax);
308 Operator_diag_precon_delete(precond);
309 break;
310 case SOLVE_METHOD_JACOBI:{
311 jacobi(A, dim, x0, rhs, maxit, flag);
312 break;
313 }
314 default:
315 assert(0);
316 break;
317
318 }
319 return res;
320}
321
322