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
25int conjugate_gradient
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
82cleanup0 :
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
98int conjugate_gradient_f
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 }
156cleanup1:
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
167int
168conjugate_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
235cleanup2 :
236 free(r);
237 free(p);
238 free(Ap);
239 free(Ax);
240 return rv;
241}
242