00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018 #include<iostream>
00019
00020
00021
00022 #include "linear_algebraic_equations.h"
00023
00024
00025 namespace OpenNN
00026 {
00027
00028
00029
00031
00032 LinearAlgebraicEquations::LinearAlgebraicEquations(void)
00033 {
00034
00035 }
00036
00037
00038
00039
00041
00042 LinearAlgebraicEquations::~LinearAlgebraicEquations(void)
00043 {
00044
00045 }
00046
00047
00048
00049
00050
00051
00054
00055 LinearAlgebraicEquations& LinearAlgebraicEquations::operator = (const LinearAlgebraicEquations&)
00056 {
00057 return(*this);
00058 }
00059
00060
00061
00062
00063
00064
00068
00069 bool LinearAlgebraicEquations::operator == (const LinearAlgebraicEquations&) const
00070 {
00071 return(true);
00072 }
00073
00074
00075
00076
00077
00078
00079
00085
00086 void LinearAlgebraicEquations::perform_Gauss_Jordan_elimination(Matrix<double>& a, Matrix<double>& b) const
00087 {
00088 unsigned int icol = 0;
00089 unsigned int irow = 0;
00090
00091 double big = 0.0;
00092 double dum = 0.0;
00093 double pivinv = 0.0;
00094
00095 unsigned int n = a.get_rows_number();
00096 unsigned int m = b.get_columns_number();
00097
00098 Vector<unsigned int> column_index(n, 0);
00099 Vector<unsigned int> row_index(n, 0);
00100 Vector<unsigned int> pivoting_index(n, 0);
00101
00102
00103
00104 for(unsigned int i = 0; i < n; i++)
00105 {
00106 big = 0.0;
00107
00108
00109
00110 for(unsigned int j = 0; j < n; j++)
00111 {
00112 if(pivoting_index[j] != 1)
00113 {
00114 for(unsigned int k = 0; k < n; k++)
00115 {
00116 if(pivoting_index[k] == 0)
00117 {
00118 if(fabs(a[j][k]) >= big)
00119 {
00120 big = fabs(a[j][k]);
00121 irow = j;
00122 icol = k;
00123 }
00124 }
00125 }
00126 }
00127 }
00128
00129 ++(pivoting_index[icol]);
00130
00131 if(irow != icol)
00132 {
00133 for(unsigned int l = 0; l < n; l++)
00134 {
00135 swap(a[irow][l],a[icol][l]);
00136 }
00137
00138 for(unsigned int l = 0; l < m; l++)
00139 {
00140 swap(b[irow][l],b[icol][l]);
00141 }
00142 }
00143
00144 row_index[i] = irow;
00145 column_index[i] = icol;
00146
00147
00148
00149 if (a[icol][icol] == 0.0)
00150 {
00151 std::ostringstream buffer;
00152
00153 buffer << "OpenNN Exception: LinearAlgebraicEquations class. "
00154 << "perform_Gauss_Jordan_elimination(Matrix<double>&, Matrix<double>&) method.\n"
00155 << "Matrix A is singular\n";
00156
00157 throw std::logic_error(buffer.str().c_str());
00158 }
00159
00160 pivinv = 1.0/a[icol][icol];
00161
00162 a[icol][icol] = 1.0;
00163
00164 for(unsigned int l = 0; l < n; l++)
00165 {
00166 a[icol][l] *= pivinv;
00167 }
00168
00169 for(unsigned int l = 0; l < m; l++)
00170 {
00171 b[icol][l] *= pivinv;
00172 }
00173
00174 for(unsigned int ll = 0; ll < n; ll++)
00175 {
00176 if(ll != icol)
00177 {
00178 dum=a[ll][icol];
00179
00180 a[ll][icol] = 0.0;
00181
00182 for(unsigned int l = 0; l < n; l++)
00183 {
00184 a[ll][l] -= a[icol][l]*dum;
00185 }
00186
00187 for(unsigned int l = 0; l < m; l++)
00188 {
00189 b[ll][l] -= b[icol][l]*dum;
00190 }
00191 }
00192 }
00193 }
00194
00195
00196
00197 for(int l = n-1; l >= 0; l--)
00198 {
00199 if(row_index[l] != column_index[l])
00200 {
00201 for(unsigned int k = 0; k < n; k++)
00202 {
00203 swap(a[k][row_index[l]], a[k][column_index[l]]);
00204 }
00205 }
00206 }
00207 }
00208
00209
00210
00211
00216
00217 void LinearAlgebraicEquations::perform_Gauss_Jordan_elimination(Matrix<double>& a, Vector<double>& b) const
00218 {
00219 unsigned int icol = 0;
00220 unsigned int irow = 0;
00221
00222 double big = 0.0;
00223 double dum = 0.0;
00224 double pivinv = 0.0;
00225
00226 unsigned int n = a.get_rows_number();
00227
00228 Vector<unsigned int> column_index(n, 0);
00229 Vector<unsigned int> row_index(n, 0);
00230 Vector<unsigned int> pivoting_index(n, 0);
00231
00232
00233
00234 for(unsigned int i = 0; i < n; i++)
00235 {
00236 big = 0.0;
00237
00238
00239
00240 for(unsigned int j = 0; j < n; j++)
00241 {
00242 if(pivoting_index[j] != 1)
00243 {
00244 for(unsigned int k = 0; k < n; k++)
00245 {
00246 if(pivoting_index[k] == 0)
00247 {
00248 if(fabs(a[j][k]) >= big)
00249 {
00250 big = fabs(a[j][k]);
00251 irow = j;
00252 icol = k;
00253 }
00254 }
00255 }
00256 }
00257 }
00258
00259 ++(pivoting_index[icol]);
00260
00261 if(irow != icol)
00262 {
00263 for(unsigned int l = 0; l < n; l++)
00264 {
00265 swap(a[irow][l],a[icol][l]);
00266 }
00267
00268 swap(b[irow],b[icol]);
00269 }
00270
00271 row_index[i] = irow;
00272 column_index[i] = icol;
00273
00274
00275
00276 if (a[icol][icol] == 0.0)
00277 {
00278 std::ostringstream buffer;
00279
00280 buffer << "OpenNN Exception: LinearAlgebraicEquations class. "
00281 << "perform_Gauss_Jordan_elimination(const Matrix<double>&, const Matrix<double>&) method.\n"
00282 << "Matrix A is singular\n";
00283
00284 throw std::logic_error(buffer.str().c_str());
00285 }
00286
00287 pivinv = 1.0/a[icol][icol];
00288
00289 a[icol][icol] = 1.0;
00290
00291 for(unsigned int l = 0; l < n; l++)
00292 {
00293 a[icol][l] *= pivinv;
00294 }
00295
00296 b[icol] *= pivinv;
00297
00298 for(unsigned int ll = 0; ll < n; ll++)
00299 {
00300 if(ll != icol)
00301 {
00302 dum=a[ll][icol];
00303
00304 a[ll][icol] = 0.0;
00305
00306 for(unsigned int l = 0; l < n; l++)
00307 {
00308 a[ll][l] -= a[icol][l]*dum;
00309 }
00310
00311 b[ll] -= b[icol]*dum;
00312 }
00313 }
00314 }
00315
00316
00317
00318 for(int l = n-1; l >= 0; l--)
00319 {
00320 if(row_index[l] != column_index[l])
00321 {
00322 for(unsigned int k = 0; k < n; k++)
00323 {
00324 swap(a[k][row_index[l]], a[k][column_index[l]]);
00325 }
00326 }
00327 }
00328 }
00329
00330
00331
00332
00336
00337 Vector<double> LinearAlgebraicEquations::calculate_Gauss_Jordan_solution(Matrix<double> a, Vector<double> b) const
00338 {
00339 const unsigned int n = a.get_rows_number();
00340
00341 Matrix<double> b_matrix(n, 1);
00342
00343 b_matrix.set_column(0, b);
00344
00345 perform_Gauss_Jordan_elimination(a, b_matrix);
00346
00347 return(b_matrix.arrange_column(0));
00348 }
00349
00350
00351
00352
00356
00357 void LinearAlgebraicEquations::swap(double& a, double& b) const
00358 {
00359 double temp = a;
00360 a = b;
00361 b = temp;
00362 }
00363
00364 }
00365
00366
00367
00368
00369
00370
00371
00372
00373
00374
00375
00376
00377
00378
00379
00380
00381
00382