Teuchos Package Browser (Single Doxygen Collection) Version of the Day
numerics/example/hilbert/cxx_main.cpp
Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 //
00004 //                    Teuchos: Common Tools Package
00005 //                 Copyright (2004) Sandia Corporation
00006 //
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 //
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00038 //
00039 // ***********************************************************************
00040 // @HEADER
00041 
00042 // Teuchos Example: Hilbert
00043   
00044 // This example showcases the usage of BLAS generics with an arbitrary precision
00045 // library -- ARPREC.
00046   
00047 // Hilbert matrices are classical examples of ill-conditioned matrices. Cholesky
00048 // factorization fails on higher-order Hilbert matrices because they lose their
00049 // positive definiteness when represented with floating-point numbers. We have
00050 // attempted to alleviate this problem with arbitrary precision.
00051   
00052 // The example program compares two datatypes, scalar type 1 and scalar type 2,
00053 // which can be customized below using #defines. Default types are mp_real (from
00054 // ARPREC) and double. The mp_real datatype must be initialized with a maximum 
00055 // precision value, also customizable below. (Default is 32.)
00056   
00057 // For a given size n, an n-by-n Hilbert matrix H and a n-by-1 std::vector b are 
00058 // constructed such that, if Hx* = b, the true solution x* is a one-std::vector.
00059 // Cholesky factorization is attempted on H; if it fails, no further tests are
00060 // attempted for that datatype. If it is successful, the approximate solution x~
00061 // is computed with a pair of BLAS TRSM (triangular solve) calls. Then, the
00062 // two-norm of (x* - x~) is computed with BLAS AXPY (std::vector update) and BLAS
00063 // NRM2. The program output is of the form:
00064 
00065 //     [size of Hilbert matrix]: [two-norm of (x* - x~)]
00066   
00067 // Tests for scalar type 2 are performed before scalar type 1 because scalar
00068 // type 2 fails at Cholesky factorization for much lower values of n if the
00069 // mp_real precision is sufficiently large.
00070   
00071 // Timing analysis still remains to be done for this example, which should be
00072 // easily accomplished with the timing mechanisms native to Teuchos.
00073 
00074 #include "Teuchos_CommandLineProcessor.hpp"
00075 #include "Teuchos_ConfigDefs.hpp"
00076 #include "Teuchos_BLAS.hpp"
00077 #include "Teuchos_Version.hpp"
00078 #include <typeinfo>
00079 
00080 #ifdef HAVE_TEUCHOS_ARPREC
00081 #include <arprec/mp_real.h>
00082 #endif
00083 
00084 #ifdef HAVE_TEUCHOS_GNU_MP
00085 #include "gmp.h"
00086 #include "gmpxx.h"
00087 #endif
00088 
00089 
00090 using namespace Teuchos;
00091 
00092 #ifdef HAVE_TEUCHOS_ARPREC
00093 #define SType1     mp_real
00094 #elif defined(HAVE_TEUCHOS_GNU_MP)
00095 #define SType1     mpf_class
00096 #else
00097 #define SType1     double
00098 #endif
00099 #define SType2     double
00100 #define OType    int
00101 
00102 template<typename TYPE>
00103 void ConstructHilbertMatrix(TYPE*, int);
00104 
00105 template<typename TYPE>
00106 void ConstructHilbertSumVector(TYPE*, int);
00107 
00108 template<typename TYPE1, typename TYPE2>
00109 void ConvertHilbertMatrix(TYPE1*, TYPE2*, int);
00110 
00111 template<typename TYPE1, typename TYPE2>
00112 void ConvertHilbertSumVector(TYPE1*, TYPE2*, int);
00113 
00114 #ifdef HAVE_TEUCHOS_ARPREC
00115 template<>
00116 void ConvertHilbertMatrix(mp_real*, double*, int);
00117 
00118 template<>
00119 void ConvertHilbertSumVector(mp_real*, double*, int);
00120 #endif
00121 
00122 #ifdef HAVE_TEUCHOS_GNU_MP
00123 template<>
00124 void ConvertHilbertMatrix(mpf_class*, double*, int);
00125 
00126 template<>
00127 void ConvertHilbertSumVector(mpf_class*, double*, int);
00128 #endif
00129 
00130 template<typename TYPE>
00131 int Cholesky(TYPE*, int);
00132 
00133 template<typename TYPE>
00134 int Solve(int, TYPE*, TYPE*, TYPE*);
00135 
00136 template<typename TYPE>
00137 void PrintArrayAsVector(TYPE*, int);
00138 
00139 template<typename TYPE>
00140 void PrintArrayAsMatrix(TYPE*, int, int);
00141 
00142 #ifdef HAVE_TEUCHOS_ARPREC
00143 template<>
00144 void PrintArrayAsVector(mp_real*, int);
00145 
00146 template<>
00147 void PrintArrayAsMatrix(mp_real*, int, int);
00148 #endif
00149 
00150 int main(int argc, char *argv[]) {
00151 
00152   std::cout << Teuchos::Teuchos_Version() << std::endl << std::endl;
00153   //
00154   // Create command line processor. 
00155   //
00156   Teuchos::CommandLineProcessor hilbertCLP(true, false);
00157   //
00158   // Set option for precision and verbosity  
00159   int precision = 32;
00160   hilbertCLP.setOption("precision", &precision, "Arbitrary precision");
00161   bool verbose = false;
00162   hilbertCLP.setOption("verbose", "quiet", &verbose, "Verbosity of example");
00163   //
00164   // Parse command line.
00165   hilbertCLP.parse( argc, argv );
00166 
00167 #ifdef HAVE_TEUCHOS_ARPREC
00168   mp::mp_init( precision );
00169 #endif
00170 
00171 #ifdef HAVE_TEUCHOS_GNU_MP
00172   mpf_set_default_prec( precision );
00173   std::cout<< "The precision of the GNU MP variable is (in bits) : "<< mpf_get_default_prec() << std::endl;
00174 #endif
00175   //
00176   // Keep track of valid datatypes
00177   //
00178   int compSType1 = 1;  // Perform cholesky factorization of matrices of SType1
00179   int compSType2 = 1;  // Perform cholesky factorization of matrices of SType2
00180   int convSType1 = 1;  // Perform cholesky factorization of matrices of SType1 (that were converted from SType2)
00181   int convSType2 = 1;  // Perform cholesky factorization of matrices of SType2 (that were converted from SType1)
00182 
00183   int n = 2;  // Initial dimension of hilbert matrix.
00184   //
00185   // Error in solution.
00186   //
00187   SType1 result1, result2_1;
00188   SType2 result2, result1_2;
00189   //
00190   // Create pointers to necessary matrices/vectors.
00191   //
00192   SType1 *H1=0, *b1=0;
00193   SType2 *H2=0, *b2=0;
00194   //
00195   while ( compSType1>0 || compSType2>0 || convSType1>0 || convSType2>0 ) {
00196     
00197     if (compSType1 > 0) {
00198       H1 = new SType1[ n*n ];
00199       b1 = new SType1[ n ];
00200       //
00201       // Construct problem.
00202       //
00203       ConstructHilbertMatrix< SType1 >(H1, n);
00204       ConstructHilbertSumVector< SType1 >(b1, n);
00205       //
00206       // Try to solve it.
00207       //
00208       compSType1 = Solve(n, H1, b1, &result1);
00209       if (compSType1 < 0 && verbose) 
00210   std::cout << typeid( result1 ).name() << " -- Cholesky factorization failed (negative diagonal) at row "<<-compSType1<< std::endl;
00211       //
00212       // Clean up always;
00213       delete [] H1; H1 = 0;
00214       delete [] b1; b1 = 0;
00215     }
00216     if (compSType2 > 0) {
00217       H2 = new SType2[ n*n ];
00218       b2 = new SType2[ n ];
00219       //
00220       // Construct problem.
00221       //
00222       ConstructHilbertMatrix< SType2 >(H2, n);
00223       ConstructHilbertSumVector< SType2 >(b2, n);
00224       //
00225       // Try to solve it.
00226       //
00227       compSType2 = Solve(n, H2, b2, &result2);
00228       if (compSType2 < 0 && verbose) 
00229   std::cout << typeid( result2 ).name() << " -- Cholesky factorization failed (negative diagonal) at row "<<-compSType2<< std::endl;
00230       //
00231       // Clean up always.
00232       delete [] H2; H2 = 0;
00233       delete [] b2; b2 = 0;      
00234     }
00235     if (convSType2 > 0) {
00236       //
00237       // Create and construct the problem in higher precision
00238       //
00239       if (!H1) H1 = new SType1[ n*n ];
00240       if (!b1) b1 = new SType1[ n ];
00241       ConstructHilbertMatrix( H1, n );
00242       ConstructHilbertSumVector( b1, n );
00243       //
00244       if (!H2) H2 = new SType2[ n*n ];
00245       if (!b2) b2 = new SType2[ n ];
00246       //
00247       // Convert the problem from SType1 to SType2 ( which should be of lesser precision )
00248       //
00249       ConvertHilbertMatrix(H1, H2, n);
00250       ConvertHilbertSumVector(b1, b2, n);
00251       //
00252       // Try to solve it.
00253       //
00254       convSType2 = Solve(n, H2, b2, &result1_2);
00255       if (convSType2 < 0 && verbose) 
00256   std::cout << typeid( result1_2 ).name() << " (converted) -- Cholesky factorization failed (negative diagonal) at row "<<-convSType2<< std::endl;
00257       //
00258       // Clean up
00259       //
00260       delete [] H2; H2 = 0;
00261       delete [] b2; b2 = 0;
00262       delete [] H1; H1 = 0;
00263       delete [] b1; b1 = 0;
00264     }
00265     if (convSType1 > 0) {
00266       //
00267       // Create and construct the problem in lower precision
00268       //
00269       if (!H2) H2 = new SType2[ n*n ];
00270       if (!b2) b2 = new SType2[ n ];
00271       ConstructHilbertMatrix(H2, n);
00272       ConstructHilbertSumVector(b2, n);
00273       //
00274       if (!H1) H1 = new SType1[ n*n ];
00275       if (!b1) b1 = new SType1[ n ];
00276       //
00277       // Convert the problem from SType2 to SType1 ( which should be of higher precision )
00278       //
00279       ConvertHilbertMatrix(H2, H1, n);
00280       ConvertHilbertSumVector(b2, b1, n);
00281       //
00282       // Try to solve it.
00283       //
00284       convSType1 = Solve(n, H1, b1, &result2_1);
00285       if (convSType1 < 0 && verbose) 
00286   std::cout << typeid( result2_1 ).name() << " (converted) -- Cholesky factorization failed (negative diagonal) at row "<<-convSType1<< std::endl;
00287       //
00288       // Clean up
00289       //
00290       delete [] H1; H1 = 0;
00291       delete [] b1; b1 = 0;
00292       delete [] H2; H2 = 0;
00293       delete [] b2; b2 = 0;
00294     }
00295     if (verbose && (compSType1>0 || compSType2>0 || convSType1>0 || convSType2>0) ) {
00296       std::cout << "***************************************************" << std::endl;
00297       std::cout << "Dimension of Hilbert Matrix : "<< n << std::endl;
00298       std::cout << "***************************************************" << std::endl;
00299       std::cout << "Datatype : Absolute error || x_hat - x ||"<< std::endl;
00300       std::cout << "---------------------------------------------------" << std::endl;
00301     }    
00302     if (compSType1>0 && verbose)
00303       std::cout << typeid( result1 ).name() << "\t : "<< result1 << std::endl;
00304     
00305     if (convSType1>0 && verbose)
00306       std::cout << typeid( result2_1 ).name() <<"(converted) : "<< result2_1 << std::endl;
00307 
00308     if (convSType2>0 && verbose)
00309       std::cout << typeid( result1_2 ).name() <<"(converted) : "<< result2_1 << std::endl;
00310 
00311     if (compSType2>0 && verbose) 
00312       std::cout << typeid( result2 ).name() << "\t : "<< result2 << std::endl;
00313     //
00314     // Increment counter.
00315     //
00316     n++;
00317   }
00318 
00319 #ifdef HAVE_TEUCHOS_ARPREC
00320   mp::mp_finalize();
00321 #endif
00322 
00323   return 0;
00324 }
00325 
00326 template<typename TYPE>
00327 void ConstructHilbertMatrix(TYPE* A, int n) {
00328   TYPE scal1 = ScalarTraits<TYPE>::one();
00329   for(int i = 0; i < n; i++) {
00330     for(int j = 0; j < n; j++) {
00331       A[i + (j * n)] = (scal1 / (i + j + scal1));
00332     }
00333   }
00334 }
00335 
00336 template<typename TYPE>
00337 void ConstructHilbertSumVector(TYPE* x, int n) {
00338   TYPE scal0 = ScalarTraits<TYPE>::zero();
00339   TYPE scal1 = ScalarTraits<TYPE>::one();
00340   for(int i = 0; i < n; i++) {
00341     x[i] = scal0;
00342     for(int j = 0; j < n; j++) {
00343       x[i] += (scal1 / (i + j + scal1));
00344     }
00345   }  
00346 }
00347 
00348 template<typename TYPE1, typename TYPE2>
00349 void ConvertHilbertMatrix(TYPE1* A, TYPE2* B, int n) {
00350   for(int i = 0; i < n; i++) {
00351     for(int j = 0; j < n; j++) {
00352       B[i + (j * n)] = A[i + (j * n)];
00353     }
00354   }
00355 }
00356 
00357 template<typename TYPE1, typename TYPE2>
00358 void ConvertHilbertSumVector(TYPE1* x, TYPE2* y, int n) {
00359   for(int i = 0; i < n; i++) {
00360     y[i] = x[i];
00361   }  
00362 }
00363 
00364 #ifdef HAVE_TEUCHOS_ARPREC
00365 template<>
00366 void ConvertHilbertMatrix(mp_real* A, double* B, int n) {
00367   for(int i = 0; i < n; i++) {
00368     for(int j = 0; j < n; j++) {
00369       B[i + (j * n)] = dble( A[i + (j * n)] );
00370     }
00371   }
00372 }
00373 
00374 template<>
00375 void ConvertHilbertSumVector(mp_real* x, double* y, int n) {
00376   for(int i = 0; i < n; i++) {
00377     y[i] = dble( x[i] );
00378   }  
00379 }
00380 #endif
00381 
00382 #ifdef HAVE_TEUCHOS_GNU_MP
00383 template<>
00384 void ConvertHilbertMatrix(mpf_class* A, double* B, int n) {
00385   for(int i = 0; i < n; i++) {
00386     for(int j = 0; j < n; j++) {
00387       B[i + (j * n)] = A[i + (j * n)].get_d();
00388     }
00389   }
00390 }
00391 
00392 template<>
00393 void ConvertHilbertSumVector(mpf_class* x, double* y, int n) {
00394   for(int i = 0; i < n; i++) {
00395     y[i] = x[i].get_d();
00396   }  
00397 }
00398 #endif
00399 
00400 template<typename TYPE>
00401 int Cholesky(TYPE* A, int n) {
00402   TYPE scal0 = ScalarTraits<TYPE>::zero();
00403   for(int k = 0; k < n; k++) {
00404     for(int j = k + 1; j < n; j++) {
00405       TYPE alpha = A[k + (j * n)] / A[k + (k * n)];
00406       for(int i = j; i < n; i++) {
00407   A[j + (i * n)] -= (alpha * A[k + (i * n)]);
00408       }
00409     }
00410     if(A[k + (k * n)] <= scal0) {
00411       return -k;
00412     }
00413     TYPE beta = ScalarTraits<TYPE>::squareroot(A[k + (k * n)]);
00414     for(int i = k; i < n; i++) {
00415       A[k + (i * n)] /= beta;
00416     }
00417   }
00418   return 1;
00419 }
00420 
00421 template<typename TYPE>
00422 int Solve(int n, TYPE* H, TYPE* b, TYPE* err) {
00423   TYPE scal0 = ScalarTraits<TYPE>::zero();
00424   TYPE scal1 = ScalarTraits<TYPE>::one();
00425   TYPE scalNeg1 = scal0 - scal1;
00426   BLAS<int, TYPE> blasObj;
00427   TYPE* x = new TYPE[n];
00428   for(int i = 0; i < n; i++) {
00429     x[i] = scal1;
00430   }
00431   int choleskySuccessful = Cholesky(H, n);
00432   if(choleskySuccessful > 0) {
00433     blasObj.TRSM(Teuchos::LEFT_SIDE, Teuchos::UPPER_TRI, Teuchos::TRANS, Teuchos::NON_UNIT_DIAG, n, 1, scal1, H, n, b, n);
00434     blasObj.TRSM(Teuchos::LEFT_SIDE, Teuchos::UPPER_TRI, Teuchos::NO_TRANS, Teuchos::NON_UNIT_DIAG, n, 1, scal1, H, n, b, n);
00435     blasObj.AXPY(n, scalNeg1, x, 1, b, 1);
00436     *err = blasObj.NRM2(n, b, 1);
00437   }  
00438   delete[] x;
00439   return choleskySuccessful;
00440 }
00441 
00442 template<typename TYPE>
00443 void PrintArrayAsVector(TYPE* x, int n) {
00444   std::cout << "[";
00445   for(int i = 0; i < n; i++) {
00446     std::cout << " " << x[i];
00447   }
00448   std::cout << " ]" << std::endl;
00449 }
00450 
00451 template<typename TYPE>
00452 void PrintArrayAsMatrix(TYPE* a, int m, int n) {
00453   std::cout << "[";
00454   for(int i = 0; i < m; i++) {
00455     if(i != 0) {
00456       std::cout << " ";
00457     }
00458     std::cout << "[";
00459     for(int j = 0; j < n; j++) {
00460       std::cout << " " << a[i + (j * m)];
00461     }
00462     std::cout << " ]";
00463     if(i != (m - 1)) {
00464       std::cout << std::endl;
00465     }
00466   }
00467   std::cout << "]" << std::endl;
00468 }
00469 
00470 #ifdef HAVE_TEUCHOS_ARPREC
00471 template<>
00472 void PrintArrayAsVector(mp_real* x, int n) {
00473   std::cout << "[ ";
00474   for(int i = 0; i < n; i++) {
00475     if(i != 0) {
00476       std::cout << "  ";
00477     }
00478     std::cout << x[i];
00479   }
00480   std::cout << "]" << std::endl;
00481 }
00482 
00483 template<>
00484 void PrintArrayAsMatrix(mp_real* a, int m, int n) {
00485   std::cout << "[";
00486   for(int i = 0; i < m; i++) {
00487     if(i != 0) {
00488       std::cout << " ";
00489     }
00490     std::cout << "[";
00491     for(int j = 0; j < n; j++) {
00492       if(j != 0) {
00493   std::cout << "  ";
00494       }
00495       std::cout << " " << a[i + (j * m)];
00496     }
00497     std::cout << " ]";
00498     if(i != (m - 1)) {
00499       std::cout << std::endl;
00500     }
00501   }
00502   std::cout << "]" << std::endl; 
00503 }
00504 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines