Rev 15347 | Details | Compare with Previous | Last modification | View Log | RSS feed
Rev | Author | Line No. | Line |
---|---|---|---|
12065 | fschmid | 1 | /* |
2 | * conjugate_gradient.cpp |
||
3 | * |
||
4 | * Copyright 2006 Nathan Hurst <njh@mail.csse.monash.edu.au> |
||
5 | * |
||
6 | * This library is free software; you can redistribute it and/or |
||
7 | * modify it either under the terms of the GNU Lesser General Public |
||
8 | * License version 2.1 as published by the Free Software Foundation |
||
9 | * (the "LGPL") or, at your option, under the terms of the Mozilla |
||
10 | * Public License Version 1.1 (the "MPL"). If you do not alter this |
||
11 | * notice, a recipient may use your version of this file under either |
||
12 | * the MPL or the LGPL. |
||
13 | * |
||
14 | * You should have received a copy of the LGPL along with this library |
||
15 | * in the file COPYING-LGPL-2.1; if not, write to the Free Software |
||
18122 | mrdocs | 16 | * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA |
12065 | fschmid | 17 | * You should have received a copy of the MPL along with this library |
18 | * in the file COPYING-MPL-1.1 |
||
19 | * |
||
20 | * The contents of this file are subject to the Mozilla Public License |
||
21 | * Version 1.1 (the "License"); you may not use this file except in |
||
22 | * compliance with the License. You may obtain a copy of the License at |
||
23 | * http://www.mozilla.org/MPL/ |
||
24 | * |
||
25 | * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY |
||
26 | * OF ANY KIND, either express or implied. See the LGPL or the MPL for |
||
27 | * the specific language governing rights and limitations. |
||
28 | * |
||
29 | */ |
||
30 | |||
31 | #include <math.h> |
||
32 | #include <stdlib.h> |
||
33 | #include <valarray> |
||
34 | #include <cassert> |
||
15347 | fschmid | 35 | #include "conjugate_gradient.h" |
12065 | fschmid | 36 | |
37 | /* lifted wholely from wikipedia. */ |
||
38 | |||
39 | using std::valarray; |
||
40 | |||
41 | static void |
||
42 | matrix_times_vector(valarray<double> const &matrix, /* m * n */ |
||
43 | valarray<double> const &vec, /* n */ |
||
44 | valarray<double> &result) /* m */ |
||
45 | { |
||
46 | unsigned n = vec.size(); |
||
47 | unsigned m = result.size(); |
||
48 | assert(m*n == matrix.size()); |
||
49 | for (unsigned i = 0; i < m; i++) { |
||
50 | double res = 0; |
||
51 | for (unsigned j = 0; j < n; j++) |
||
15347 | fschmid | 52 | res += matrix[i*m+j] * vec[j]; |
12065 | fschmid | 53 | result[i] = res; |
54 | } |
||
55 | } |
||
15347 | fschmid | 56 | /* |
12065 | fschmid | 57 | static double Linfty(valarray<double> const &vec) { |
58 | return std::max(vec.max(), -vec.min()); |
||
59 | } |
||
15347 | fschmid | 60 | */ |
12065 | fschmid | 61 | double |
62 | inner(valarray<double> const &x, |
||
63 | valarray<double> const &y) { |
||
64 | double total = 0; |
||
65 | for(unsigned i = 0; i < x.size(); i++) |
||
66 | total += x[i]*y[i]; |
||
67 | return total;// (x*y).sum(); <- this is more concise, but ineff |
||
68 | } |
||
69 | |||
70 | void |
||
71 | conjugate_gradient(double **A, |
||
72 | double *x, |
||
73 | double *b, |
||
74 | unsigned n, |
||
75 | double tol, |
||
76 | int max_iterations, |
||
77 | bool ortho1) { |
||
78 | valarray<double> vA(n*n); |
||
79 | valarray<double> vx(n); |
||
80 | valarray<double> vb(n); |
||
81 | for(unsigned i=0;i<n;i++) { |
||
82 | vx[i]=x[i]; |
||
83 | vb[i]=b[i]; |
||
84 | for(unsigned j=0;j<n;j++) { |
||
85 | vA[i*n+j]=A[i][j]; |
||
86 | } |
||
87 | } |
||
88 | conjugate_gradient(vA,vx,vb,n,tol,max_iterations,ortho1); |
||
89 | for(unsigned i=0;i<n;i++) { |
||
90 | x[i]=vx[i]; |
||
91 | } |
||
92 | } |
||
93 | void |
||
94 | conjugate_gradient(valarray<double> const &A, |
||
95 | valarray<double> &x, |
||
96 | valarray<double> const &b, |
||
97 | unsigned n, double tol, |
||
15347 | fschmid | 98 | unsigned max_iterations, bool ortho1) { |
12065 | fschmid | 99 | valarray<double> Ap(n), p(n), r(n); |
100 | matrix_times_vector(A,x,Ap); |
||
101 | r=b-Ap; |
||
102 | double r_r = inner(r,r); |
||
103 | unsigned k = 0; |
||
104 | tol *= tol; |
||
105 | while(k < max_iterations && r_r > tol) { |
||
106 | k++; |
||
107 | double r_r_new = r_r; |
||
108 | if(k == 1) |
||
109 | p = r; |
||
110 | else { |
||
111 | r_r_new = inner(r,r); |
||
112 | p = r + (r_r_new/r_r)*p; |
||
113 | } |
||
114 | matrix_times_vector(A, p, Ap); |
||
115 | double alpha_k = r_r_new / inner(p, Ap); |
||
116 | x += alpha_k*p; |
||
117 | r -= alpha_k*Ap; |
||
118 | r_r = r_r_new; |
||
119 | } |
||
120 | //printf("njh: %d iters, Linfty = %g L2 = %g\n", k, |
||
121 | //std::max(-r.min(), r.max()), sqrt(r_r)); |
||
122 | // x is solution |
||
123 | } |
||
124 | /* |
||
125 | Local Variables: |
||
126 | mode:c++ |
||
127 | c-file-style:"stroustrup" |
||
128 | c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +)) |
||
129 | indent-tabs-mode:nil |
||
130 | fill-column:99 |
||
131 | End: |
||
132 | */ |
||
15347 | fschmid | 133 | // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=4:softtabstop=4 |