Graphviz  2.39.20141222.0545
conjgrad.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 
15 #include "matrix_ops.h"
16 #include "conjgrad.h"
17 /* #include <math.h> */
18 #include <stdlib.h>
19 
20 
21 /*************************
22 ** C.G. method - SPARSE *
23 *************************/
24 
26  (vtx_data * A, double *x, double *b, int n, double tol,
27  int max_iterations) {
28  /* Solves Ax=b using Conjugate-Gradients method */
29  /* 'x' and 'b' are orthogonalized against 1 */
30 
31  int i, rv = 0;
32 
33  double alpha, beta, r_r, r_r_new, p_Ap;
34  double *r = N_GNEW(n, double);
35  double *p = N_GNEW(n, double);
36  double *Ap = N_GNEW(n, double);
37  double *Ax = N_GNEW(n, double);
38  double *alphap = N_GNEW(n, double);
39 
40  double *orth_b = N_GNEW(n, double);
41  copy_vector(n, b, orth_b);
42  orthog1(n, orth_b);
43  orthog1(n, x);
44  right_mult_with_vector(A, n, x, Ax);
45  vectors_subtraction(n, orth_b, Ax, r);
46  copy_vector(n, r, p);
47  r_r = vectors_inner_product(n, r, r);
48 
49  for (i = 0; i < max_iterations && max_abs(n, r) > tol; i++) {
50  right_mult_with_vector(A, n, p, Ap);
51  p_Ap = vectors_inner_product(n, p, Ap);
52  if (p_Ap == 0)
53  break; /*exit(1); */
54  alpha = r_r / p_Ap;
55 
56  /* derive new x: */
57  vectors_scalar_mult(n, p, alpha, alphap);
58  vectors_addition(n, x, alphap, x);
59 
60  /* compute values for next iteration: */
61  if (i < max_iterations - 1) { /* not last iteration */
62  vectors_scalar_mult(n, Ap, alpha, Ap);
63  vectors_subtraction(n, r, Ap, r); /* fast computation of r, the residual */
64 
65  /* Alternaive accurate, but slow, computation of the residual - r */
66  /* right_mult_with_vector(A, n, x, Ax); */
67  /* vectors_subtraction(n,b,Ax,r); */
68 
69  r_r_new = vectors_inner_product(n, r, r);
70  if (r_r == 0) {
71  agerr (AGERR, "conjugate_gradient: unexpected length 0 vector\n");
72  rv = 1;
73  goto cleanup0;
74  }
75  beta = r_r_new / r_r;
76  r_r = r_r_new;
77  vectors_scalar_mult(n, p, beta, p);
78  vectors_addition(n, r, p, p);
79  }
80  }
81 
82 cleanup0 :
83  free(r);
84  free(p);
85  free(Ap);
86  free(Ax);
87  free(alphap);
88  free(orth_b);
89 
90  return rv;
91 }
92 
93 
94 /****************************
95 ** C.G. method - DENSE *
96 ****************************/
97 
99  (float **A, double *x, double *b, int n, double tol,
100  int max_iterations, boolean ortho1) {
101  /* Solves Ax=b using Conjugate-Gradients method */
102  /* 'x' and 'b' are orthogonalized against 1 if 'ortho1=true' */
103 
104  int i, rv = 0;
105 
106  double alpha, beta, r_r, r_r_new, p_Ap;
107  double *r = N_GNEW(n, double);
108  double *p = N_GNEW(n, double);
109  double *Ap = N_GNEW(n, double);
110  double *Ax = N_GNEW(n, double);
111  double *alphap = N_GNEW(n, double);
112 
113  double *orth_b = N_GNEW(n, double);
114  copy_vector(n, b, orth_b);
115  if (ortho1) {
116  orthog1(n, orth_b);
117  orthog1(n, x);
118  }
119  right_mult_with_vector_f(A, n, x, Ax);
120  vectors_subtraction(n, orth_b, Ax, r);
121  copy_vector(n, r, p);
122  r_r = vectors_inner_product(n, r, r);
123 
124  for (i = 0; i < max_iterations && max_abs(n, r) > tol; i++) {
125  right_mult_with_vector_f(A, n, p, Ap);
126  p_Ap = vectors_inner_product(n, p, Ap);
127  if (p_Ap == 0)
128  break; /*exit(1); */
129  alpha = r_r / p_Ap;
130 
131  /* derive new x: */
132  vectors_scalar_mult(n, p, alpha, alphap);
133  vectors_addition(n, x, alphap, x);
134 
135  /* compute values for next iteration: */
136  if (i < max_iterations - 1) { /* not last iteration */
137  vectors_scalar_mult(n, Ap, alpha, Ap);
138  vectors_subtraction(n, r, Ap, r); /* fast computation of r, the residual */
139 
140  /* Alternaive accurate, but slow, computation of the residual - r */
141  /* right_mult_with_vector(A, n, x, Ax); */
142  /* vectors_subtraction(n,b,Ax,r); */
143 
144  r_r_new = vectors_inner_product(n, r, r);
145  if (r_r == 0) {
146  rv = 1;
147  agerr (AGERR, "conjugate_gradient: unexpected length 0 vector\n");
148  goto cleanup1;
149  }
150  beta = r_r_new / r_r;
151  r_r = r_r_new;
152  vectors_scalar_mult(n, p, beta, p);
153  vectors_addition(n, r, p, p);
154  }
155  }
156 cleanup1:
157  free(r);
158  free(p);
159  free(Ap);
160  free(Ax);
161  free(alphap);
162  free(orth_b);
163 
164  return rv;
165 }
166 
167 int
168 conjugate_gradient_mkernel(float *A, float *x, float *b, int n,
169  double tol, int max_iterations)
170 {
171  /* Solves Ax=b using Conjugate-Gradients method */
172  /* A is a packed symmetric matrix */
173  /* matrux A is "packed" (only upper triangular portion exists, row-major); */
174 
175  int i, rv = 0;
176 
177  double alpha, beta, r_r, r_r_new, p_Ap;
178  float *r = N_NEW(n, float);
179  float *p = N_NEW(n, float);
180  float *Ap = N_NEW(n, float);
181  float *Ax = N_NEW(n, float);
182 
183  /* centering x and b */
184  orthog1f(n, x);
185  orthog1f(n, b);
186 
187  right_mult_with_vector_ff(A, n, x, Ax);
188  /* centering Ax */
189  orthog1f(n, Ax);
190 
191 
192  vectors_substractionf(n, b, Ax, r);
193  copy_vectorf(n, r, p);
194 
195  r_r = vectors_inner_productf(n, r, r);
196 
197  for (i = 0; i < max_iterations && max_absf(n, r) > tol; i++) {
198  orthog1f(n, p);
199  orthog1f(n, x);
200  orthog1f(n, r);
201 
202  right_mult_with_vector_ff(A, n, p, Ap);
203  /* centering Ap */
204  orthog1f(n, Ap);
205 
206  p_Ap = vectors_inner_productf(n, p, Ap);
207  if (p_Ap == 0)
208  break;
209  alpha = r_r / p_Ap;
210 
211  /* derive new x: */
212  vectors_mult_additionf(n, x, (float) alpha, p);
213 
214  /* compute values for next iteration: */
215  if (i < max_iterations - 1) { /* not last iteration */
216  vectors_mult_additionf(n, r, (float) -alpha, Ap);
217 
218 
219  r_r_new = vectors_inner_productf(n, r, r);
220 
221  if (r_r == 0) {
222  rv = 1;
223  agerr (AGERR, "conjugate_gradient: unexpected length 0 vector\n");
224  goto cleanup2;
225  }
226  beta = r_r_new / r_r;
227  r_r = r_r_new;
228 
229  vectors_scalar_multf(n, p, (float) beta, p);
230 
231  vectors_additionf(n, r, p, p);
232  }
233  }
234 
235 cleanup2 :
236  free(r);
237  free(p);
238  free(Ap);
239  free(Ax);
240  return rv;
241 }
void vectors_additionf(int n, float *vector1, float *vector2, float *result)
Definition: matrix_ops.c:629
Definition: cgraph.h:389
void orthog1(int n, double *vec)
Definition: matrix_ops.c:319
#define N_NEW(n, t)
Definition: memory.h:36
void copy_vectorf(int n, float *source, float *dest)
Definition: matrix_ops.c:657
void vectors_addition(int n, double *vector1, double *vector2, double *result)
Definition: matrix_ops.c:399
void copy_vector(int n, double *source, double *dest)
Definition: matrix_ops.c:431
int agerr(agerrlevel_t level, const char *fmt,...)
Definition: agerror.c:142
void vectors_mult_additionf(int n, float *vector1, float alpha, float *vector2)
Definition: matrix_ops.c:639
void right_mult_with_vector_f(float **matrix, int n, double *vector, double *result)
Definition: matrix_ops.c:371
int conjugate_gradient(vtx_data *A, double *x, double *b, int n, double tol, int max_iterations)
Definition: conjgrad.c:26
void vectors_scalar_multf(int n, float *vector, float alpha, float *result)
Definition: matrix_ops.c:648
double vectors_inner_productf(int n, float *vector1, float *vector2)
Definition: matrix_ops.c:665
void free()
int i
Definition: gvdevice.c:448
double vectors_inner_product(int n, double *vector1, double *vector2)
Definition: matrix_ops.c:439
#define alpha
Definition: shapes.c:3864
void orthog1f(int n, float *vec)
Definition: matrix_ops.c:544
void vectors_substractionf(int n, float *vector1, float *vector2, float *result)
Definition: matrix_ops.c:619
int conjugate_gradient_f(float **A, double *x, double *b, int n, double tol, int max_iterations, boolean ortho1)
Definition: conjgrad.c:99
void vectors_subtraction(int n, double *vector1, double *vector2, double *result)
Definition: matrix_ops.c:388
void right_mult_with_vector_ff(float *packed_matrix, int n, float *vector, float *result)
Definition: matrix_ops.c:594
void vectors_scalar_mult(int n, double *vector, double alpha, double *result)
Definition: matrix_ops.c:422
int conjugate_gradient_mkernel(float *A, float *x, float *b, int n, double tol, int max_iterations)
Definition: conjgrad.c:168
#define N_GNEW(n, t)
Definition: agxbuf.c:20
void right_mult_with_vector(vtx_data *matrix, int n, double *vector, double *result)
Definition: matrix_ops.c:354