Subversion Repositories Scribus

Rev

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