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