27 using std::setprecision;
38 const std::vector< double >& y )
43 z.resize( nrows, 0.0 );
45 for (
int i = 0; i < nrows; i++)
53 const std::vector< double >& y )
58 z.resize( nrows, 0.0 );
60 for (
int i = 0; i < nrows; i++)
72 z.resize( nrows, 0.0 );
74 for (
int i = 0; i < nrows; i++ )
80 std::vector< vector< double > >
81 operator + (
const std::vector< std::vector< double > >&A,
82 const std::vector< std::vector< double > >&B )
85 int ncols = A[0].size();
87 vector< vector< double > > C( nrows );
88 for(
int i = 0; i < nrows; i++ )
89 C[i].resize( ncols, 0.0 );
91 for (
int i = 0; i < nrows; i++)
92 for (
int j = 0; j < ncols; j++)
93 C[i][j] = A[i][j] + B[i][j];
98 std::vector< vector< double > >
99 operator - (
const std::vector< std::vector< double > >&A,
100 const std::vector< std::vector< double > >&B )
102 int nrows = A.size();
103 int ncols = A[0].size();
105 vector< vector< double > > C( nrows );
106 for(
int i = 0; i < nrows; i++ )
107 C[i].resize( ncols, 0.0 );
109 for (
int i = 0; i < nrows; i++)
110 for (
int j = 0; j < ncols; j++)
111 C[i][j] = A[i][j] - B[i][j];
117 std::vector< double >
120 int nrows = x.size();
123 y.resize( nrows, 0.0 );
125 for (
int i = 0; i < nrows; i++)
131 std::vector< double >
134 int nrows = x.size();
137 assert( abs( a ) > DBL_EPSILON );
139 y.resize( nrows, 0.0 );
141 for (
int i = 0; i < nrows; i++)
147 std::vector< std::vector< double > >
148 operator * (
double a,
const std::vector< std::vector< double > >&A )
150 int nrows = A.size();
151 int ncols = A[0].size();
153 vector< vector< double > > B( nrows );
154 for(
int i = 0; i < nrows; i++ )
155 B[i].resize( ncols, 0.0 );
157 for (
int i = 0; i < nrows; i++)
158 for (
int j = 0; j < ncols; j++)
159 B[i][j] = a * A[i][j];
164 std::vector< std::vector< double > >
165 operator / (
const std::vector< std::vector< double > >& A,
double a )
167 int nrows = A.size();
168 int ncols = A[0].size();
170 assert( abs( a ) > DBL_EPSILON );
172 vector< vector< double > > B( nrows );
173 for(
int i = 0; i < nrows; i++ )
174 B[i].resize( ncols, 0.0 );
176 for (
int i = 0; i < nrows; i++)
177 for (
int j = 0; j < ncols; j++)
183 std::vector< double >
185 const std::vector< double >& x )
187 int nrows = A.size();
188 int ncols = A[0].size();
191 y.resize( nrows, 0.0 );
193 for (
int i = 0; i < nrows; i++)
194 for (
int j = 0; j < ncols; j++)
195 y[i] += A[i][j] * x[j];
200 std::vector< double >
202 const std::vector< std::vector< double > >& A )
204 int nrows = A.size();
205 int ncols = A[0].size();
208 y.resize( ncols, 0.0 );
210 for (
int j = 0; j < ncols; j++)
211 for (
int i = 0; i < nrows; i++)
212 y[j] += A[i][j] * x[i];
217 std::vector< vector< double > >
219 const std::vector< std::vector< double > >&B )
225 vector< vector< double > > C( m );
226 for(
int i = 0; i < m; i++ )
227 C[i].resize( n, 0.0 );
229 for (
int i = 0; i < m; i++)
230 for (
int j = 0; j < n; j++)
231 for (
int k = 0; k < r; k++)
232 C[i][j] += A[i][k] * B[k][j];
239 const std::vector< double > & b )
243 for (
unsigned int i = 0; i < a.size(); i++ )
250 vector< vector< double > >
252 const std::vector< double >& b )
254 vector< vector< double> > C( a.size() );
255 for(
unsigned int i = 0; i < a.size(); i++ )
256 C[i].resize( b.size() );
258 for(
unsigned int i = 0; i < a.size(); i++ )
259 for(
unsigned int j = 0; j < b.size(); j++ )
260 C[i][j] = a[i] * b[j];
268 const std::vector< double > x )
271 unsigned int n = A.size();
273 assert ( x.size() == n );
275 for (
unsigned int i = 0; i < n; i++ )
276 for (
unsigned int j = 0; j < n; j++ )
277 prod += x[i] * A[i][j] * x[j];
284 norm (
const std::vector< double > & a )
293 unsigned int n = A.size();
295 for (
unsigned int i = 0; i < n ; ++i )
296 for (
unsigned int j = 0; j <= i ; ++j)
298 double sum = A[i][j];
302 for (
unsigned int k = 0; k < j; ++k )
303 sum -= A[i][k] * A[j][k];
307 if (sum < 0)
return EXIT_FAILURE;
311 if ( fabs(sum) < DBL_EPSILON )
return EXIT_FAILURE;
316 A[i][j] = sum / A[j][j];
325 std::vector< double > & x,
326 const std::vector< double > & b )
328 unsigned int n = L.size();
333 for (
unsigned int i = 0; i < n ; i++)
336 for (
unsigned int j = 0; j < i ; ++j)
337 sum -= x[j] * L[i][j];
338 x[i] = sum / L[i][i];
342 for (
int i = n - 1; i >= 0; i-- )
345 for (
unsigned int j = i + 1; j < n ; j++)
346 sum -= x[j] * L[j][i];
347 x[i] = sum / L[i][i];
356 std::vector < std::vector < double > > & Ainv )
358 unsigned int n = A.size();
359 vector< double > x, b;
360 vector< vector< double > > L;
369 for (
unsigned int i = 0; i < n; i++ )
372 L[i].resize ( n, 0.0 );
375 Ainv[i].resize ( n, 0.0 );
378 for (
unsigned int i = 0; i < n; i++ )
379 for (
unsigned int j = 0; j < n; j++ )
386 for (
unsigned int j = 0; j < n; j++ )
389 for (
unsigned int i = 0; i < n; i++ )
390 b[i] = ( i == j ) ? 1 : 0;
396 for (
unsigned int i = 0; i < n; i++ )
404 eye ( std::vector < std::vector < double > >& I,
int n )
407 for(
int i = 0; i < n; i++ )
410 I[i].resize( n, 0.0 );
413 for(
int i = 0; i < n; i++ )
420 write (
const std::vector < double > & a )
422 unsigned int n = a.size();
424 for (
unsigned int i = 0; i < n; ++i )
425 cout << setprecision( 15 ) << a[i] << endl;
433 write (
const std::vector < std::vector < double > > & A )
435 unsigned int n = A.size();
437 for (
unsigned int i = 0; i < n; ++i )
439 for (
unsigned int j = 0; j < n; ++j )
440 cout << setw( 8 ) << setprecision( 4 ) << A[i][j];
455 for(
int i = 0; i < m; i++ )
458 A[i].resize( n, 0.0 );