MechSys: linalg/tst/tlinsolver.cpp Source File

00001 /*************************************************************************************
00002  * MechSys - A C++ library to simulate (Continuum) Mechanical Systems                *
00003  * Copyright (C) 2005 Dorival de Moraes Pedroso <dorival.pedroso at gmail.com>       *
00004  * Copyright (C) 2005 Raul Dario Durand Farfan  <raul.durand at gmail.com>           *
00005  *                                                                                   *
00006  * This file is part of MechSys.                                                     *
00007  *                                                                                   *
00008  * MechSys is free software; you can redistribute it and/or modify it under the      *
00009  * terms of the GNU General Public License as published by the Free Software         *
00010  * Foundation; either version 2 of the License, or (at your option) any later        *
00011  * version.                                                                          *
00012  *                                                                                   *
00013  * MechSys is distributed in the hope that it will be useful, but WITHOUT ANY        *
00014  * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A   *
00015  * PARTICULAR PURPOSE. See the GNU General Public License for more details.          *
00016  *                                                                                   *
00017  * You should have received a copy of the GNU General Public License along with      *
00018  * MechSys; if not, write to the Free Software Foundation, Inc., 51 Franklin Street, *
00019  * Fifth Floor, Boston, MA 02110-1301, USA                                           *
00020  *************************************************************************************/
00021 
00022 // STL
00023 #include <iostream>
00024 #include <cmath>
00025 #include <iomanip>
00026 
00027 // MPI (must be before #define REAL double)
00028 #ifdef USE_SCALAPACK
00029   #include <mpi.h>
00030 #endif
00031 
00032 // MechSys
00033 #include "linalg/lawrap.h"
00034 #include "linalg/matrix.h"
00035 #include "linalg/vector.h"
00036 #include "util/exception.h"
00037 
00038 using std::cout;
00039 using std::endl;
00040 using LinAlg::Matrix;
00041 using LinAlg::Vector;
00042 
00043 int SolveAndOutput(Matrix<REAL> & A, Vector<REAL> & Y, Vector<REAL> const & Ycorr)
00044 {
00045     // Tolerance
00046     const double MAXERR = 1.0e-8;
00047 
00048     // Size
00049     int sz = A.Rows();
00050 
00051     // Solve
00052     LinAlg::Gesv(A,Y); // Y <- inv(A)*Y
00053 
00054 #ifdef USE_SCALAPACK
00055     int rank = MPI::COMM_WORLD.Get_rank(); // This process
00056     if (rank==0)
00057 #endif
00058     {
00059         // Output
00060         cout << sz << " x " << sz << endl;
00061         cout << "Y     = " << Y;
00062         cout << "Ycorr = " << Ycorr << endl; // Ycorrect
00063 
00064         // Check
00065         bool ok = true;
00066         for (int i=0; i<sz; ++i) if (fabs(Y(i)-Ycorr(i))>MAXERR) ok=false;
00067         if (ok) cout << "..ok.." << endl;
00068         else    cout << ".. FAILED .... FAILED .... FAILED .... FAILED .... FAILED .... FAILED .." << endl;
00069         cout << "\n\n" << endl;
00070         return (ok ? 0 : 1);
00071     }
00072     return 0;
00073 }
00074 
00075 int main(int argc, char **argv) try
00076 {
00077 
00078 #ifdef USE_SCALAPACK
00079     // Initialize MPI
00080     MPI::Init(argc, argv);
00081 #endif
00082 
00083     int num_errors=0;
00084 
00085     {
00086         Matrix<REAL> A(9,9," 19   3   1   12   1   16   1   3  11 \
00087                             -19   3   1   12   1   16   1   3  11 \
00088                             -19  -3   1   12   1   16   1   3  11 \
00089                             -19  -3  -1   12   1   16   1   3  11 \
00090                             -19  -3  -1  -12   1   16   1   3  11 \
00091                             -19  -3  -1  -12  -1   16   1   3  11 \
00092                             -19  -3  -1  -12  -1  -16   1   3  11 \
00093                             -19   3  -1  -12  -1  -16  -1   3  11 \
00094                             -19  -3  -1  -12  -1  -16  -1  -3  11");
00095         Vector<REAL> Y    (9,"0 0 1  0 0 0  0 0 0");
00096         Vector<REAL> Ycorr(9,"0  -0.16666666666666666 0.5 0  0  0  -0.5 0.16666666666666666 0");
00097         num_errors += SolveAndOutput(A,Y,Ycorr);
00098     }
00099 
00100     {
00101         Matrix<REAL> A(5,5,"2  3  0  0  0 \
00102                             3  0  4  0  6 \
00103                             0 -1 -3  2  0 \
00104                             0  0  1  0  0 \
00105                             0  4  2  0  1");
00106         Vector<REAL> Y    (5,"8 45 -3 3 19");
00107         Vector<REAL> Ycorr(5, "1 2 3 4 5");
00108         num_errors += SolveAndOutput(A,Y,Ycorr);
00109     }
00110 
00111     {
00112         Matrix<REAL> A(5,5,"5  4  3  2  1 \
00113                             4  4  3  2  1 \
00114                             0  3  3  2  1 \
00115                             0  0  2  2  1 \
00116                             0  0  0  1  1");
00117         Vector<REAL> Y    (5,"1 2 3 4 5");
00118         Vector<REAL> Ycorr(5, "-1 3 -10  19 -14");
00119         num_errors += SolveAndOutput(A,Y,Ycorr);
00120     }
00121 
00122     {
00123         Matrix<REAL> A(6,6,"6  5  4  3  2  1 \
00124                             5  5  4  3  2  1 \
00125                             0  4  4  3  2  1 \
00126                             0  0  3  3  2  1 \
00127                             0  0  0  2  2  1 \
00128                             0  0  0  0  1  1");
00129         Vector<REAL> Y    (6,"1 2 3 4 5 6");
00130         Vector<REAL> Ycorr(6,"-1  4  -17  50 -101  107");
00131         num_errors += SolveAndOutput(A,Y,Ycorr);
00132     }
00133 
00134     {
00135         Matrix<REAL> A(7,7,"7  6  5  4  3  2  1 \
00136                             6  6  5  4  3  2  1 \
00137                             0  5  5  4  3  2  1 \
00138                             0  0  4  4  3  2  1 \
00139                             0  0  0  3  3  2  1 \
00140                             0  0  0  0  2  2  1 \
00141                             0  0  0  0  0  1  1"); 
00142         Vector<REAL> Y    (7,"1 2 3 4 5 6 7");
00143         Vector<REAL> Ycorr(7,"-1    5  -26  103 -310  619 -612");
00144         num_errors += SolveAndOutput(A,Y,Ycorr);
00145     }
00146 
00147     {
00148         Matrix<REAL> A(10,10,"10  9  8  7  6  5  4  3  2  1 \
00149                                9  9  8  7  6  5  4  3  2  1 \
00150                                0  8  8  7  6  5  4  3  2  1 \
00151                                0  0  7  7  6  5  4  3  2  1 \
00152                                0  0  0  6  6  5  4  3  2  1 \
00153                                0  0  0  0  5  5  4  3  2  1 \
00154                                0  0  0  0  0  4  4  3  2  1 \
00155                                0  0  0  0  0  0  3  3  2  1 \
00156                                0  0  0  0  0  0  0  2  2  1 \
00157                                0  0  0  0  0  0  0  0  1  1");
00158         Vector<REAL> Y    (10,"1 2 3 4 5 6 7 8 9 10");
00159         Vector<REAL> Ycorr(10,"-1       8     -65     454   -2725   13624  -54497  163490 -326981 326991");
00160         num_errors += SolveAndOutput(A,Y,Ycorr);
00161     }
00162 
00163     {
00164         Matrix<REAL> A(20,20,"0                    0                    0                   3.5675854203801425  5.586999952339394    0                   9.124672011124929    9.145908826954795   1.2416508676904114  6.38300935148842    8.372597363621663   0                    4.8937055932552855   2.122155506388934   1.9249910701053075  8.10075353261826    0                    0                   0                   0                  \
00165                               0                    8.034018687286288    9.138640099458767   0                   7.167512496608924    0                   0                    0.5224156775260869  9.151075423706791   0.5527127803943133  0                   0                    0                    0                   2.272425659398155   0                   0.09821762809327228  0                   1.0346169405287153  0                  \
00166                               3.933323734304129    0                    6.045976642762936   0.6280320819039942  0                    0                   0.6140203425194934   6.39310615769223    0                   4.037854158235598   1.9017865519581822  3.602733303541066    0                    2.686808823172111   5.9502664007032156  0                   4.002550260487107    4.627875905089286   1.4056849404831984  0                  \
00167                               0                    9.660110167917344    0                   5.075437326736201   8.229344625037534    0                   8.583627552501572    1.7807198014928283  0                   1.6191894089316983  3.180176627991914   0.5028035436711342   0                    0                   0                   0                   0                    1.1694508123497538  0                   4.973393581488593  \
00168                               0                    0                    0                   0                   8.955387440376533    0                   0                    0                   0                   6.737033231793585   0                   0                    0.8406995803220718   0                   0                   4.3216109960992855  0                    3.4671405554185952  0                   0                  \
00169                               1.385736985109698    0                    0                   0                   0                    0                   5.756620211148549    8.43296076919799    5.689564996693641   0                   0                   0                    0                    0                   0                   0                   0                    0                   0                   3.4900697970019823 \
00170                               2.0789969292331802   5.974867554471115    4.338542348294354   8.046348353368327   0.08751824121803531  0                   0.43480370689359615  0                   0                   0                   0                   0.21197661565830694  5.878441983814603    0                   0                   0                   0                    6.788854452000153   6.062662267233091   6.728162823511978  \
00171                               0.5374390217351777   8.771369977500436    0                   1.8018373020167933  0                    5.246951774894724   0                    0                   4.200189050235592   0                   5.824521414137147   7.705584595537237    8.316324494690875    0                   0                   0                   0                    0                   1.277254092094513   7.574874413062603  \
00172                               0                    3.9584190446893897   0                   5.013657729990476   0                    6.385530108566956   0                    0                   0                   0                   8.277569130450104   6.899156203121713    3.0655986254352285   3.302314983284859   0                   5.656808001497019   0                    3.562520664465325   0.920620182277091   0                  \
00173                               2.9652891618821853   9.605098735961636    5.969015261347524   8.486840233026586   6.570461155812925    2.7689233965970916  4.891827503661292    6.103225568342411   0                   1.192249182768621   9.755511421595996   0                    1.315292257739018    0                   0                   5.035028083029985   0                    0                   7.71746739264559    0                  \
00174                               3.236407921267915    0                    0                   0                   0.1393753667108022   7.857758401045775   0.11713860216390648  1.4631471861831646  8.147193514176983   0                   0                   8.32428646975454     0                    0                   0                   9.360562742016489   0                    0                   3.0181093147037608  0.1813351619034509 \
00175                               0                    0                    7.7372331399768415  0                   4.300505563006509    0                   0                    0                   0                   8.967138515878984   0                   0                    0                    0                   5.516559493051853   0                   0                    0                   0                   7.5308398601016755 \
00176                               0                    0                    0                   0                   0.7260314156986958   9.641906392874098   8.647227787313017    0                   4.843826767243984   7.162499107286214   0.7161359833354453  0                    5.055131381762084    0.8666467595989003  0                   1.0712736864402794  0                    4.6832902662542     6.741092859798396   0                  \
00177                               6.723683426451984    1.7347989745713865   2.81269142489325    3.161818033414278   3.789589775557017    0.2973896849528046  9.636707282777255    0                   7.243793402664997   0                   9.141605456332742   0                    0                    1.7563111778436968  0                   0                   0                    0                   0                   9.107769992791393  \
00178                               7.660177348481911    9.757878945492967    0                   0                   5.454970391258724    0                   1.480182691130142    0                   3.980512003770281   5.912202052303711   5.275390362376423   9.402134883591792    1.5325609264828544   0                   0                   0.9688270177141955  8.961998004688787    9.682075344044918   2.3331578962435398  0                  \
00179                               0                    0.7420651996198635   5.499413067550481   0                   0                    0                   8.985616883119826    0                   0                   0                   1.632453790716828   6.233979849609247    0                    6.179287061888416   2.8706643104880083  0                   0                    0                   0                   0                  \
00180                               0                    0                    8.522563946980364   6.962833350821894   6.390638362572334    7.7503519153029155  6.6588096153608      0                   3.148540643401213   0                   0                   0                    2.7875828171548225   0                   0                   0                   6.355005714351777    4.733037768787803   0.683004053735865   0                  \
00181                               0                    2.654866147252987    0                   0.2777425437946157  0                    4.871086319563457   0                    0                   0.2486323604990659  6.157788721206641   9.262154654959115   0                    0.02717699005374219  0                   0                   6.033192161053535   0                    4.740589637450165   6.601074976942364   0.6917756252120999 \
00182                               0.22916959815091897  3.9489130728041766   0                   3.816313079140624   7.486380981073855    0                   2.5056964723817443   0                   4.839742316423194   0.309707325646863   0                   1.3273283407419922   0                    7.707075561098562   8.049161871827437   0                   7.586628280796931    3.544465739650361   3.133452745731952   9.449910658959729  \
00183                               0                    0.44951693816252636  0                   2.696042434612286   9.721936217783874    0                   4.732090227606945    9.77128819033402    0                   0                   3.7389230532001907  0                    7.470756543006561    0                   0                   0.3021467296140301  4.713176433137849    7.67146049488812    0                   3.053629908626001  ");
00184         Vector<REAL> Y    (20,"1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20");
00185         Vector<REAL> Ycorr(20,"-5.83298035975012485466 -5.02528090315008135747 -0.65339477675371748777 \
00186                                 3.28342871409988967812  2.17233316234833395697 -3.02582636229451429344 \
00187                                 0.49178772642741530596 -0.61316205636004328383  2.81474223668284828648 \
00188                                 0.26549722703766331922  3.03791749623312412609  6.32611698147693246597 \
00189                                -2.48386676081233570557 -4.50044353516062134446  0.80704125950832916736 \
00190                                -3.97526809879191134200 -0.57278944457015512626  0.87247668887712592767 \
00191                                 5.23101362560551041980  0.11691861137762919742");
00192         num_errors += SolveAndOutput(A,Y,Ycorr);
00193     }
00194     
00195     {
00196         Matrix<REAL> A(30,30,"30  29  28  27  26  25  24  23  22  21  20  19  18  17  16  15  14  13  12  11  10  9  8  7  6  5  4  3  2  1 \
00197                               29  29  28  27  26  25  24  23  22  21  20  19  18  17  16  15  14  13  12  11  10  9  8  7  6  5  4  3  2  1 \
00198                                0  28  28  27  26  25  24  23  22  21  20  19  18  17  16  15  14  13  12  11  10  9  8  7  6  5  4  3  2  1 \
00199                                0   0  27  27  26  25  24  23  22  21  20  19  18  17  16  15  14  13  12  11  10  9  8  7  6  5  4  3  2  1 \
00200                                0   0   0  26  26  25  24  23  22  21  20  19  18  17  16  15  14  13  12  11  10  9  8  7  6  5  4  3  2  1 \
00201                                0   0   0   0  25  25  24  23  22  21  20  19  18  17  16  15  14  13  12  11  10  9  8  7  6  5  4  3  2  1 \
00202                                0   0   0   0   0  24  24  23  22  21  20  19  18  17  16  15  14  13  12  11  10  9  8  7  6  5  4  3  2  1 \
00203                                0   0   0   0   0   0  23  23  22  21  20  19  18  17  16  15  14  13  12  11  10  9  8  7  6  5  4  3  2  1 \
00204                                0   0   0   0   0   0   0  22  22  21  20  19  18  17  16  15  14  13  12  11  10  9  8  7  6  5  4  3  2  1 \
00205                                0   0   0   0   0   0   0   0  21  21  20  19  18  17  16  15  14  13  12  11  10  9  8  7  6  5  4  3  2  1 \
00206                                0   0   0   0   0   0   0   0   0  20  20  19  18  17  16  15  14  13  12  11  10  9  8  7  6  5  4  3  2  1 \
00207                                0   0   0   0   0   0   0   0   0   0  19  19  18  17  16  15  14  13  12  11  10  9  8  7  6  5  4  3  2  1 \
00208                                0   0   0   0   0   0   0   0   0   0   0  18  18  17  16  15  14  13  12  11  10  9  8  7  6  5  4  3  2  1 \
00209                                0   0   0   0   0   0   0   0   0   0   0   0  17  17  16  15  14  13  12  11  10  9  8  7  6  5  4  3  2  1 \
00210                                0   0   0   0   0   0   0   0   0   0   0   0   0  16  16  15  14  13  12  11  10  9  8  7  6  5  4  3  2  1 \
00211                                0   0   0   0   0   0   0   0   0   0   0   0   0   0  15  15  14  13  12  11  10  9  8  7  6  5  4  3  2  1 \
00212                                0   0   0   0   0   0   0   0   0   0   0   0   0   0   0  14  14  13  12  11  10  9  8  7  6  5  4  3  2  1 \
00213                                0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0  13  13  12  11  10  9  8  7  6  5  4  3  2  1 \
00214                                0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0  12  12  11  10  9  8  7  6  5  4  3  2  1 \
00215                                0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0  11  11  10  9  8  7  6  5  4  3  2  1 \
00216                                0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0  10  10  9  8  7  6  5  4  3  2  1 \
00217                                0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   9  9  8  7  6  5  4  3  2  1 \
00218                                0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0  8  8  7  6  5  4  3  2  1 \
00219                                0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0  0  7  7  6  5  4  3  2  1 \
00220                                0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0  0  0  6  6  5  4  3  2  1 \
00221                                0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0  0  0  0  5  5  4  3  2  1 \
00222                                0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0  0  0  0  0  4  4  3  2  1 \
00223                                0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0  0  0  0  0  0  3  3  2  1 \
00224                                0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0  0  0  0  0  0  0  2  2  1 \
00225                                0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0  0  0  0  0  0  0  0  1  1");
00226         Vector<REAL> Y(30,"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");
00227     }
00228 
00229     cout << "Number of errors = " << num_errors << endl;
00230 
00231 #ifdef USE_SCALAPACK
00232     // Finalize MPI
00233     MPI::Finalize();
00234 #endif
00235 
00236     return 0;
00237 }
00238 catch (Exception * e) //{{{ 
00239 {
00240     e->Cout();
00241     if (e->IsFatal()) {delete e; exit(1);}
00242     delete e;
00243 }
00244 catch (char const * m)
00245 {
00246     std::cout << "Fatal: " << m << std::endl;
00247     exit (1);
00248 }
00249 catch (...)
00250 {
00251     std::cout << "Some exception (...) ocurred\n";
00252 } //}}} 
00253 
00254 // vim:fdm=marker