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 | |
25 | int 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 | |
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 | |
98 | int 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 | } |
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 | } |
242 | |