1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
|
#include <iostream>
#include <iomanip>
#include <string>
#include <vector>
#include <algorithm>
#include <numeric>
#include <cmath>
#include <cassert>
using namespace std;
const double SMALL = 1.0E-30; // used to stop divide-by-zero
const double NEARZERO = 1.0E-10; // interpretation of "zero" for printing purposes
using vec = vector<double>; // vector
using matrix = vector<vec>; // matrix (=collection of (row) vectors)
//======================================================================
void print( const string &title, const matrix &A )
{
if ( title != "" ) cout << title << '\n';
for ( auto &row : A )
{
for ( auto x : row ) cout << setw( 15 ) << ( abs( x ) < NEARZERO ? 0.0 : x );
cout << '\n';
}
}
//======================================================================
void print( const string &title, const vec &A )
{
if ( title != "" ) cout << title << '\n';
for ( auto x : A ) cout << setw( 15 ) << ( abs( x ) < NEARZERO ? 0.0 : x );
cout << '\n';
}
//======================================================================
matrix matmul( const matrix &A, const matrix &B ) // Matrix times matrix
{
int rowsA = A.size(), colsA = A[0].size();
int rowsB = B.size(), colsB = B[0].size();
assert( colsA == rowsB );
matrix C( rowsA, vec( colsB, 0.0 ) );
for ( int i = 0; i < rowsA; i++ )
{
for ( int j = 0; j < colsB; j++ )
{
// scalar product of ith row of A and jth column of B
for ( int k = 0; k < colsA; k++ ) C[i][j] += A[i][k] * B[k][j];
}
}
return C;
}
//======================================================================
vec matmul( const matrix &A, const vec &V ) // Matrix times vector
{
int rowsA = A.size(), colsA = A[0].size();
int rowsV = V.size();
assert( colsA == rowsV );
vec C( rowsA, 0.0 );
for ( int i = 0; i < rowsA; i++ )
{
for ( int k = 0; k < colsA; k++ ) C[i] += A[i][k] * V[k];
}
return C;
}
//======================================================================
double innerProduct( const vec &U, const vec &V ) // Inner product of U and V
{
return inner_product( U.begin(), U.end(), V.begin(), 0.0 );
}
//======================================================================
double norm( const vec &V ) // L2 norm of V
{
return sqrt( innerProduct( V, V ) );
}
//======================================================================
bool conjugateGradient( const matrix &A, const vec &B, vec &X )
{
double TOLERANCE = 1.0e-10;
int n = A.size();
X = vec( n, 0.0 );
vec R = B, P = B;
for ( int k = 0; k < n; k++ )
{
vec Rold = R; // Store previous residual
vec AP = matmul( A, P );
double alpha = innerProduct( R, R ) / max( innerProduct( P, AP ), SMALL );
for ( int i = 0; i < n; i++ )
{
X[i] += alpha * P[i]; // Next estimate of solution
R[i] -= alpha * AP[i]; // Residual
}
if ( norm( R ) < TOLERANCE ) break; // Convergence test
double beta = innerProduct( R, R ) / max( innerProduct( Rold, Rold ), SMALL );
for ( int i = 0; i < n; i++ ) P[i] = R[i] + beta * P[i]; // Next gradient
}
return true;
}
//======================================================================
int main()
{
matrix A = { { 7, 3, 1 },
{ 3, 7, 0 },
{ 1, 0, 10 } };
vec B = { 1, -5, 2 };
vec X;
if ( conjugateGradient( A, B, X ) )
{
print( "X:" , X );
cout << "\n\nCheck:\n";
print( "B:" , B );
print( "AX:", matmul( A, X ) );
}
else
{
cout << "Unable to solve\n";
}
}
|
