Teuchos Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
numerics/example/hilbert/cxx_main.cpp
Go to the documentation of this file.
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Teuchos: Common Tools Package
5 // Copyright (2004) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38 //
39 // ***********************************************************************
40 // @HEADER
41 
42 // Teuchos Example: Hilbert
43 
44 // This example showcases the usage of BLAS generics with an arbitrary precision
45 // library -- ARPREC.
46 
47 // Hilbert matrices are classical examples of ill-conditioned matrices. Cholesky
48 // factorization fails on higher-order Hilbert matrices because they lose their
49 // positive definiteness when represented with floating-point numbers. We have
50 // attempted to alleviate this problem with arbitrary precision.
51 
52 // The example program compares two datatypes, scalar type 1 and scalar type 2,
53 // which can be customized below using #defines. Default types are mp_real (from
54 // ARPREC) and double. The mp_real datatype must be initialized with a maximum
55 // precision value, also customizable below. (Default is 32.)
56 
57 // For a given size n, an n-by-n Hilbert matrix H and a n-by-1 std::vector b are
58 // constructed such that, if Hx* = b, the true solution x* is a one-std::vector.
59 // Cholesky factorization is attempted on H; if it fails, no further tests are
60 // attempted for that datatype. If it is successful, the approximate solution x~
61 // is computed with a pair of BLAS TRSM (triangular solve) calls. Then, the
62 // two-norm of (x* - x~) is computed with BLAS AXPY (std::vector update) and BLAS
63 // NRM2. The program output is of the form:
64 
65 // [size of Hilbert matrix]: [two-norm of (x* - x~)]
66 
67 // Tests for scalar type 2 are performed before scalar type 1 because scalar
68 // type 2 fails at Cholesky factorization for much lower values of n if the
69 // mp_real precision is sufficiently large.
70 
71 // Timing analysis still remains to be done for this example, which should be
72 // easily accomplished with the timing mechanisms native to Teuchos.
73 
75 #include "Teuchos_ConfigDefs.hpp"
76 #include "Teuchos_BLAS.hpp"
77 #include "Teuchos_Version.hpp"
78 #include <typeinfo>
79 
80 #ifdef HAVE_TEUCHOS_ARPREC
81 #include <arprec/mp_real.h>
82 #endif
83 
84 #ifdef HAVE_TEUCHOS_GNU_MP
85 #include "gmp.h"
86 #include "gmpxx.h"
87 #endif
88 
89 
90 using namespace Teuchos;
91 
92 #ifdef HAVE_TEUCHOS_ARPREC
93 #define SType1 mp_real
94 #elif defined(HAVE_TEUCHOS_GNU_MP)
95 #define SType1 mpf_class
96 #else
97 #define SType1 double
98 #endif
99 #define SType2 double
100 #define OType int
101 
102 template<typename TYPE>
103 void ConstructHilbertMatrix(TYPE*, int);
104 
105 template<typename TYPE>
106 void ConstructHilbertSumVector(TYPE*, int);
107 
108 template<typename TYPE1, typename TYPE2>
109 void ConvertHilbertMatrix(TYPE1*, TYPE2*, int);
110 
111 template<typename TYPE1, typename TYPE2>
112 void ConvertHilbertSumVector(TYPE1*, TYPE2*, int);
113 
114 #ifdef HAVE_TEUCHOS_ARPREC
115 template<>
116 void ConvertHilbertMatrix(mp_real*, double*, int);
117 
118 template<>
119 void ConvertHilbertSumVector(mp_real*, double*, int);
120 #endif
121 
122 #ifdef HAVE_TEUCHOS_GNU_MP
123 template<>
124 void ConvertHilbertMatrix(mpf_class*, double*, int);
125 
126 template<>
127 void ConvertHilbertSumVector(mpf_class*, double*, int);
128 #endif
129 
130 template<typename TYPE>
131 int Cholesky(TYPE*, int);
132 
133 template<typename TYPE>
134 int Solve(int, TYPE*, TYPE*, TYPE*);
135 
136 template<typename TYPE>
137 void PrintArrayAsVector(TYPE*, int);
138 
139 template<typename TYPE>
140 void PrintArrayAsMatrix(TYPE*, int, int);
141 
142 #ifdef HAVE_TEUCHOS_ARPREC
143 template<>
144 void PrintArrayAsVector(mp_real*, int);
145 
146 template<>
147 void PrintArrayAsMatrix(mp_real*, int, int);
148 #endif
149 
150 int main(int argc, char *argv[]) {
151 
152  std::cout << Teuchos::Teuchos_Version() << std::endl << std::endl;
153  //
154  // Create command line processor.
155  //
156  Teuchos::CommandLineProcessor hilbertCLP(true, false);
157  //
158  // Set option for precision and verbosity
159  int precision = 32;
160  hilbertCLP.setOption("precision", &precision, "Arbitrary precision");
161  bool verbose = false;
162  hilbertCLP.setOption("verbose", "quiet", &verbose, "Verbosity of example");
163  //
164  // Parse command line.
165  hilbertCLP.parse( argc, argv );
166 
167 #ifdef HAVE_TEUCHOS_ARPREC
168  mp::mp_init( precision );
169 #endif
170 
171 #ifdef HAVE_TEUCHOS_GNU_MP
172  mpf_set_default_prec( precision );
173  std::cout<< "The precision of the GNU MP variable is (in bits) : "<< mpf_get_default_prec() << std::endl;
174 #endif
175  //
176  // Keep track of valid datatypes
177  //
178  int compSType1 = 1; // Perform cholesky factorization of matrices of SType1
179  int compSType2 = 1; // Perform cholesky factorization of matrices of SType2
180  int convSType1 = 1; // Perform cholesky factorization of matrices of SType1 (that were converted from SType2)
181  int convSType2 = 1; // Perform cholesky factorization of matrices of SType2 (that were converted from SType1)
182 
183  int n = 2; // Initial dimension of hilbert matrix.
184  //
185  // Error in solution.
186  //
187  SType1 result1, result2_1;
188  SType2 result2, result1_2;
189  //
190  // Create pointers to necessary matrices/vectors.
191  //
192  SType1 *H1=0, *b1=0;
193  SType2 *H2=0, *b2=0;
194  //
195  while ( compSType1>0 || compSType2>0 || convSType1>0 || convSType2>0 ) {
196 
197  if (compSType1 > 0) {
198  H1 = new SType1[ n*n ];
199  b1 = new SType1[ n ];
200  //
201  // Construct problem.
202  //
203  ConstructHilbertMatrix< SType1 >(H1, n);
204  ConstructHilbertSumVector< SType1 >(b1, n);
205  //
206  // Try to solve it.
207  //
208  compSType1 = Solve(n, H1, b1, &result1);
209  if (compSType1 < 0 && verbose)
210  std::cout << typeid( result1 ).name() << " -- Cholesky factorization failed (negative diagonal) at row "<<-compSType1<< std::endl;
211  //
212  // Clean up always;
213  delete [] H1; H1 = 0;
214  delete [] b1; b1 = 0;
215  }
216  if (compSType2 > 0) {
217  H2 = new SType2[ n*n ];
218  b2 = new SType2[ n ];
219  //
220  // Construct problem.
221  //
222  ConstructHilbertMatrix< SType2 >(H2, n);
223  ConstructHilbertSumVector< SType2 >(b2, n);
224  //
225  // Try to solve it.
226  //
227  compSType2 = Solve(n, H2, b2, &result2);
228  if (compSType2 < 0 && verbose)
229  std::cout << typeid( result2 ).name() << " -- Cholesky factorization failed (negative diagonal) at row "<<-compSType2<< std::endl;
230  //
231  // Clean up always.
232  delete [] H2; H2 = 0;
233  delete [] b2; b2 = 0;
234  }
235  if (convSType2 > 0) {
236  //
237  // Create and construct the problem in higher precision
238  //
239  if (!H1) H1 = new SType1[ n*n ];
240  if (!b1) b1 = new SType1[ n ];
241  ConstructHilbertMatrix( H1, n );
242  ConstructHilbertSumVector( b1, n );
243  //
244  if (!H2) H2 = new SType2[ n*n ];
245  if (!b2) b2 = new SType2[ n ];
246  //
247  // Convert the problem from SType1 to SType2 ( which should be of lesser precision )
248  //
249  ConvertHilbertMatrix(H1, H2, n);
250  ConvertHilbertSumVector(b1, b2, n);
251  //
252  // Try to solve it.
253  //
254  convSType2 = Solve(n, H2, b2, &result1_2);
255  if (convSType2 < 0 && verbose)
256  std::cout << typeid( result1_2 ).name() << " (converted) -- Cholesky factorization failed (negative diagonal) at row "<<-convSType2<< std::endl;
257  //
258  // Clean up
259  //
260  delete [] H2; H2 = 0;
261  delete [] b2; b2 = 0;
262  delete [] H1; H1 = 0;
263  delete [] b1; b1 = 0;
264  }
265  if (convSType1 > 0) {
266  //
267  // Create and construct the problem in lower precision
268  //
269  if (!H2) H2 = new SType2[ n*n ];
270  if (!b2) b2 = new SType2[ n ];
271  ConstructHilbertMatrix(H2, n);
273  //
274  if (!H1) H1 = new SType1[ n*n ];
275  if (!b1) b1 = new SType1[ n ];
276  //
277  // Convert the problem from SType2 to SType1 ( which should be of higher precision )
278  //
279  ConvertHilbertMatrix(H2, H1, n);
280  ConvertHilbertSumVector(b2, b1, n);
281  //
282  // Try to solve it.
283  //
284  convSType1 = Solve(n, H1, b1, &result2_1);
285  if (convSType1 < 0 && verbose)
286  std::cout << typeid( result2_1 ).name() << " (converted) -- Cholesky factorization failed (negative diagonal) at row "<<-convSType1<< std::endl;
287  //
288  // Clean up
289  //
290  delete [] H1; H1 = 0;
291  delete [] b1; b1 = 0;
292  delete [] H2; H2 = 0;
293  delete [] b2; b2 = 0;
294  }
295  if (verbose && (compSType1>0 || compSType2>0 || convSType1>0 || convSType2>0) ) {
296  std::cout << "***************************************************" << std::endl;
297  std::cout << "Dimension of Hilbert Matrix : "<< n << std::endl;
298  std::cout << "***************************************************" << std::endl;
299  std::cout << "Datatype : Absolute error || x_hat - x ||"<< std::endl;
300  std::cout << "---------------------------------------------------" << std::endl;
301  }
302  if (compSType1>0 && verbose)
303  std::cout << typeid( result1 ).name() << "\t : "<< result1 << std::endl;
304 
305  if (convSType1>0 && verbose)
306  std::cout << typeid( result2_1 ).name() <<"(converted) : "<< result2_1 << std::endl;
307 
308  if (convSType2>0 && verbose)
309  std::cout << typeid( result1_2 ).name() <<"(converted) : "<< result2_1 << std::endl;
310 
311  if (compSType2>0 && verbose)
312  std::cout << typeid( result2 ).name() << "\t : "<< result2 << std::endl;
313  //
314  // Increment counter.
315  //
316  n++;
317  }
318 
319 #ifdef HAVE_TEUCHOS_ARPREC
320  mp::mp_finalize();
321 #endif
322 
323  return 0;
324 }
325 
326 template<typename TYPE>
327 void ConstructHilbertMatrix(TYPE* A, int n) {
328  TYPE scal1 = ScalarTraits<TYPE>::one();
329  for(int i = 0; i < n; i++) {
330  for(int j = 0; j < n; j++) {
331  A[i + (j * n)] = (scal1 / (i + j + scal1));
332  }
333  }
334 }
335 
336 template<typename TYPE>
337 void ConstructHilbertSumVector(TYPE* x, int n) {
338  TYPE scal0 = ScalarTraits<TYPE>::zero();
339  TYPE scal1 = ScalarTraits<TYPE>::one();
340  for(int i = 0; i < n; i++) {
341  x[i] = scal0;
342  for(int j = 0; j < n; j++) {
343  x[i] += (scal1 / (i + j + scal1));
344  }
345  }
346 }
347 
348 template<typename TYPE1, typename TYPE2>
349 void ConvertHilbertMatrix(TYPE1* A, TYPE2* B, int n) {
350  for(int i = 0; i < n; i++) {
351  for(int j = 0; j < n; j++) {
352  B[i + (j * n)] = A[i + (j * n)];
353  }
354  }
355 }
356 
357 template<typename TYPE1, typename TYPE2>
358 void ConvertHilbertSumVector(TYPE1* x, TYPE2* y, int n) {
359  for(int i = 0; i < n; i++) {
360  y[i] = x[i];
361  }
362 }
363 
364 #ifdef HAVE_TEUCHOS_ARPREC
365 template<>
366 void ConvertHilbertMatrix(mp_real* A, double* B, int n) {
367  for(int i = 0; i < n; i++) {
368  for(int j = 0; j < n; j++) {
369  B[i + (j * n)] = dble( A[i + (j * n)] );
370  }
371  }
372 }
373 
374 template<>
375 void ConvertHilbertSumVector(mp_real* x, double* y, int n) {
376  for(int i = 0; i < n; i++) {
377  y[i] = dble( x[i] );
378  }
379 }
380 #endif
381 
382 #ifdef HAVE_TEUCHOS_GNU_MP
383 template<>
384 void ConvertHilbertMatrix(mpf_class* A, double* B, int n) {
385  for(int i = 0; i < n; i++) {
386  for(int j = 0; j < n; j++) {
387  B[i + (j * n)] = A[i + (j * n)].get_d();
388  }
389  }
390 }
391 
392 template<>
393 void ConvertHilbertSumVector(mpf_class* x, double* y, int n) {
394  for(int i = 0; i < n; i++) {
395  y[i] = x[i].get_d();
396  }
397 }
398 #endif
399 
400 template<typename TYPE>
401 int Cholesky(TYPE* A, int n) {
402  TYPE scal0 = ScalarTraits<TYPE>::zero();
403  for(int k = 0; k < n; k++) {
404  for(int j = k + 1; j < n; j++) {
405  TYPE alpha = A[k + (j * n)] / A[k + (k * n)];
406  for(int i = j; i < n; i++) {
407  A[j + (i * n)] -= (alpha * A[k + (i * n)]);
408  }
409  }
410  if(A[k + (k * n)] <= scal0) {
411  return -k;
412  }
413  TYPE beta = ScalarTraits<TYPE>::squareroot(A[k + (k * n)]);
414  for(int i = k; i < n; i++) {
415  A[k + (i * n)] /= beta;
416  }
417  }
418  return 1;
419 }
420 
421 template<typename TYPE>
422 int Solve(int n, TYPE* H, TYPE* b, TYPE* err) {
423  TYPE scal0 = ScalarTraits<TYPE>::zero();
424  TYPE scal1 = ScalarTraits<TYPE>::one();
425  TYPE scalNeg1 = scal0 - scal1;
426  BLAS<int, TYPE> blasObj;
427  TYPE* x = new TYPE[n];
428  for(int i = 0; i < n; i++) {
429  x[i] = scal1;
430  }
431  int choleskySuccessful = Cholesky(H, n);
432  if(choleskySuccessful > 0) {
435  blasObj.AXPY(n, scalNeg1, x, 1, b, 1);
436  *err = blasObj.NRM2(n, b, 1);
437  }
438  delete[] x;
439  return choleskySuccessful;
440 }
441 
442 template<typename TYPE>
443 void PrintArrayAsVector(TYPE* x, int n) {
444  std::cout << "[";
445  for(int i = 0; i < n; i++) {
446  std::cout << " " << x[i];
447  }
448  std::cout << " ]" << std::endl;
449 }
450 
451 template<typename TYPE>
452 void PrintArrayAsMatrix(TYPE* a, int m, int n) {
453  std::cout << "[";
454  for(int i = 0; i < m; i++) {
455  if(i != 0) {
456  std::cout << " ";
457  }
458  std::cout << "[";
459  for(int j = 0; j < n; j++) {
460  std::cout << " " << a[i + (j * m)];
461  }
462  std::cout << " ]";
463  if(i != (m - 1)) {
464  std::cout << std::endl;
465  }
466  }
467  std::cout << "]" << std::endl;
468 }
469 
470 #ifdef HAVE_TEUCHOS_ARPREC
471 template<>
472 void PrintArrayAsVector(mp_real* x, int n) {
473  std::cout << "[ ";
474  for(int i = 0; i < n; i++) {
475  if(i != 0) {
476  std::cout << " ";
477  }
478  std::cout << x[i];
479  }
480  std::cout << "]" << std::endl;
481 }
482 
483 template<>
484 void PrintArrayAsMatrix(mp_real* a, int m, int n) {
485  std::cout << "[";
486  for(int i = 0; i < m; i++) {
487  if(i != 0) {
488  std::cout << " ";
489  }
490  std::cout << "[";
491  for(int j = 0; j < n; j++) {
492  if(j != 0) {
493  std::cout << " ";
494  }
495  std::cout << " " << a[i + (j * m)];
496  }
497  std::cout << " ]";
498  if(i != (m - 1)) {
499  std::cout << std::endl;
500  }
501  }
502  std::cout << "]" << std::endl;
503 }
504 #endif
void TRSM(ESide side, EUplo uplo, ETransp transa, EDiag diag, const OrdinalType &m, const OrdinalType &n, const alpha_type alpha, const A_type *A, const OrdinalType &lda, ScalarType *B, const OrdinalType &ldb) const
Solves the matrix equations: op(A)*X=alpha*B or X*op(A)=alpha*B where X and B are m by n matrices...
void AXPY(const OrdinalType &n, const alpha_type alpha, const x_type *x, const OrdinalType &incx, ScalarType *y, const OrdinalType &incy) const
Perform the operation: y &lt;- y+alpha*x.
void ConstructHilbertSumVector(TYPE *, int)
Templated interface class to BLAS routines.
static T squareroot(T x)
Returns a number of magnitudeType that is the square root of this scalar type x.
Teuchos header file which uses auto-configuration information to include necessary C++ headers...
void PrintArrayAsMatrix(TYPE *, int, int)
Templated BLAS wrapper.
ScalarTraits< ScalarType >::magnitudeType NRM2(const OrdinalType &n, const ScalarType *x, const OrdinalType &incx) const
Compute the 2-norm of the vector x.
void PrintArrayAsVector(TYPE *, int)
void ConvertHilbertSumVector(TYPE1 *, TYPE2 *, int)
void setOption(const char option_true[], const char option_false[], bool *option_val, const char documentation[]=NULL)
Set a boolean option.
EParseCommandLineReturn parse(int argc, char *argv[], std::ostream *errout=&std::cerr) const
Parse a command line.
std::string Teuchos_Version()
int main(int argc, char *argv[])
void ConvertHilbertMatrix(TYPE1 *, TYPE2 *, int)
Basic command line parser for input from (argc,argv[])
void ConstructHilbertMatrix(TYPE *, int)
static T zero()
Returns representation of zero for this scalar type.
int Cholesky(TYPE *, int)
int Solve(int, TYPE *, TYPE *, TYPE *)
static T one()
Returns representation of one for this scalar type.
Class that helps parse command line input arguments from (argc,argv[]) and set options.