Graphviz  2.41.20170921.2350
sparse_solve.c
Go to the documentation of this file.
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 
30 };
31 
32 
34  FREE(o->data);
35 }
36 
39  SparseMatrix A = d->A;
40  real alpha = d->alpha;
41  real xsum = 0.;
42  int m = A->m, i;
43 
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 
57  Operator o;
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;
65  return o;
66 }
67 
68 
72  return y;
73 }
74 
76  Operator o;
77 
78  o = GNEW(struct Operator_struct);
79  o->data = (void*) A;
81  return o;
82 }
83 
84 
86  if (o) FREE(o);
87 }
88 
89 
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 
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 
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 
124 
125  return o;
126 }
127 
128 
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 
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 
153 
154  return o;
155 }
156 
158  if (o->data) FREE(o->data);
159  if (o) FREE(o);
160 }
161 
162 static 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 
227 real 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 
249 real* 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);
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 
296 real 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);
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 
#define ABS(a)
Definition: arith.h:45
real * vector_saxpy(int n, real *x, real *y, real beta)
Definition: general.c:108
#define FREE
Definition: PriorityQueue.c:23
#define assert(x)
Definition: cghdr.h:47
struct SparseMatrix_struct * SparseMatrix
Definition: SparseMatrix.h:40
real * vector_saxpy2(int n, real *x, real *y, real beta)
Definition: general.c:115
Operator Operator_diag_precon_new(SparseMatrix A)
Definition: sparse_solve.c:129
real * vector_subtract_to(int n, real *x, real *y)
Definition: general.c:88
void SparseMatrix_multiply_vector(SparseMatrix A, real *v, real **res, int transposed)
int conjugate_gradient(vtx_data *A, double *x, double *b, int n, double tol, int max_iterations)
Definition: conjgrad.c:26
void Operator_diag_precon_delete(Operator o)
Definition: sparse_solve.c:157
real * Operator_diag_precon_apply(Operator o, real *x, real *y)
Definition: sparse_solve.c:90
real * Operator_uniform_stress_matmul_apply(Operator o, real *x, real *y)
Definition: sparse_solve.c:37
Operator Operator_matmul_new(SparseMatrix A)
Definition: sparse_solve.c:75
#define MALLOC
Definition: PriorityQueue.c:21
real vector_product(int n, real *x, real *y)
Definition: general.c:101
Operator Operator_uniform_stress_diag_precon_new(SparseMatrix A, real alpha)
Definition: sparse_solve.c:100
real * jacobi(SparseMatrix A, int dim, real *x0, real *rhs, int maxit, int *flag)
Definition: sparse_solve.c:249
#define MEMCPY
Definition: PriorityQueue.c:24
real cg(Operator Ax, Operator precond, int n, int dim, real *x0, real *rhs, real tol, int maxit, int *flag)
Definition: sparse_solve.c:227
#define GNEW(t)
Definition: memory.h:37
real SparseMatrix_solve(SparseMatrix A, int dim, real *x0, real *rhs, real tol, int maxit, int method, int *flag)
Definition: sparse_solve.c:296
EXTERN unsigned char Verbose
Definition: globals.h:64
for(;;)
Definition: grammar.c:1846
Operator Operator_uniform_stress_matmul(SparseMatrix A, real alpha)
Definition: sparse_solve.c:56
real *(* Operator_apply)(Operator o, real *in, real *out)
Definition: sparse_solve.h:26
#define alpha
Definition: shapes.c:3902
void Operator_uniform_stress_matmul_delete(Operator o)
Definition: sparse_solve.c:33
real * Operator_matmul_apply(Operator o, real *x, real *y)
Definition: sparse_solve.c:69
#define N_GNEW(n, t)
Definition: agxbuf.c:20
#define FALSE
Definition: cgraph.h:35
void Operator_matmul_delete(Operator o)
Definition: sparse_solve.c:85
#define real
Definition: general.h:34