Teuchos Package Browser (Single Doxygen Collection)  Version of the Day
numerics/example/hilbert/cxx_main.cpp
Go to the documentation of this file.
2 // *****************************************************************************
3 // Teuchos: Common Tools Package
4 //
5 // Copyright 2004 NTESS and the Teuchos contributors.
7 // *****************************************************************************
9
10 // Teuchos Example: Hilbert
11
12 // This example showcases the usage of BLAS generics with an arbitrary precision
13 // library -- ARPREC.
14
15 // Hilbert matrices are classical examples of ill-conditioned matrices. Cholesky
16 // factorization fails on higher-order Hilbert matrices because they lose their
17 // positive definiteness when represented with floating-point numbers. We have
18 // attempted to alleviate this problem with arbitrary precision.
19
20 // The example program compares two datatypes, scalar type 1 and scalar type 2,
21 // which can be customized below using #defines. Default types are mp_real (from
22 // ARPREC) and double. The mp_real datatype must be initialized with a maximum
23 // precision value, also customizable below. (Default is 32.)
24
25 // For a given size n, an n-by-n Hilbert matrix H and a n-by-1 std::vector b are
26 // constructed such that, if Hx* = b, the true solution x* is a one-std::vector.
27 // Cholesky factorization is attempted on H; if it fails, no further tests are
28 // attempted for that datatype. If it is successful, the approximate solution x~
29 // is computed with a pair of BLAS TRSM (triangular solve) calls. Then, the
30 // two-norm of (x* - x~) is computed with BLAS AXPY (std::vector update) and BLAS
31 // NRM2. The program output is of the form:
32
33 // [size of Hilbert matrix]: [two-norm of (x* - x~)]
34
35 // Tests for scalar type 2 are performed before scalar type 1 because scalar
36 // type 2 fails at Cholesky factorization for much lower values of n if the
37 // mp_real precision is sufficiently large.
38
39 // Timing analysis still remains to be done for this example, which should be
40 // easily accomplished with the timing mechanisms native to Teuchos.
41
43 #include "Teuchos_ConfigDefs.hpp"
44 #include "Teuchos_BLAS.hpp"
45 #include "Teuchos_Version.hpp"
46 #include <typeinfo>
47
48 #ifdef HAVE_TEUCHOS_ARPREC
49 #include <arprec/mp_real.h>
50 #endif
51
52 #ifdef HAVE_TEUCHOS_GNU_MP
53 #include "gmp.h"
54 #include "gmpxx.h"
55 #endif
56
57
58 using namespace Teuchos;
59
60 #ifdef HAVE_TEUCHOS_ARPREC
61 #define SType1 mp_real
62 #elif defined(HAVE_TEUCHOS_GNU_MP)
63 #define SType1 mpf_class
64 #else
65 #define SType1 double
66 #endif
67 #define SType2 double
68 #define OType int
69
70 template<typename TYPE>
71 void ConstructHilbertMatrix(TYPE*, int);
72
73 template<typename TYPE>
74 void ConstructHilbertSumVector(TYPE*, int);
75
76 template<typename TYPE1, typename TYPE2>
77 void ConvertHilbertMatrix(TYPE1*, TYPE2*, int);
78
79 template<typename TYPE1, typename TYPE2>
80 void ConvertHilbertSumVector(TYPE1*, TYPE2*, int);
81
82 #ifdef HAVE_TEUCHOS_ARPREC
83 template<>
84 void ConvertHilbertMatrix(mp_real*, double*, int);
85
86 template<>
87 void ConvertHilbertSumVector(mp_real*, double*, int);
88 #endif
89
90 #ifdef HAVE_TEUCHOS_GNU_MP
91 template<>
92 void ConvertHilbertMatrix(mpf_class*, double*, int);
93
94 template<>
95 void ConvertHilbertSumVector(mpf_class*, double*, int);
96 #endif
97
98 template<typename TYPE>
99 int Cholesky(TYPE*, int);
100
101 template<typename TYPE>
102 int Solve(int, TYPE*, TYPE*, TYPE*);
103
104 template<typename TYPE>
105 void PrintArrayAsVector(TYPE*, int);
106
107 template<typename TYPE>
108 void PrintArrayAsMatrix(TYPE*, int, int);
109
110 #ifdef HAVE_TEUCHOS_ARPREC
111 template<>
112 void PrintArrayAsVector(mp_real*, int);
113
114 template<>
115 void PrintArrayAsMatrix(mp_real*, int, int);
116 #endif
117
118 int main(int argc, char *argv[]) {
119
120  std::cout << Teuchos::Teuchos_Version() << std::endl << std::endl;
121  //
122  // Create command line processor.
123  //
124  Teuchos::CommandLineProcessor hilbertCLP(true, false);
125  //
126  // Set option for precision and verbosity
127  int precision = 32;
128  hilbertCLP.setOption("precision", &precision, "Arbitrary precision");
129  bool verbose = false;
130  hilbertCLP.setOption("verbose", "quiet", &verbose, "Verbosity of example");
131  //
132  // Parse command line.
133  hilbertCLP.parse( argc, argv );
134
135 #ifdef HAVE_TEUCHOS_ARPREC
136  mp::mp_init( precision );
137 #endif
138
139 #ifdef HAVE_TEUCHOS_GNU_MP
140  mpf_set_default_prec( precision );
141  std::cout<< "The precision of the GNU MP variable is (in bits) : "<< mpf_get_default_prec() << std::endl;
142 #endif
143  //
144  // Keep track of valid datatypes
145  //
146  int compSType1 = 1; // Perform cholesky factorization of matrices of SType1
147  int compSType2 = 1; // Perform cholesky factorization of matrices of SType2
148  int convSType1 = 1; // Perform cholesky factorization of matrices of SType1 (that were converted from SType2)
149  int convSType2 = 1; // Perform cholesky factorization of matrices of SType2 (that were converted from SType1)
150
151  int n = 2; // Initial dimension of hilbert matrix.
152  //
153  // Error in solution.
154  //
155  SType1 result1, result2_1;
156  SType2 result2, result1_2;
157  //
158  // Create pointers to necessary matrices/vectors.
159  //
160  SType1 *H1=0, *b1=0;
161  SType2 *H2=0, *b2=0;
162  //
163  while ( compSType1>0 || compSType2>0 || convSType1>0 || convSType2>0 ) {
164
165  if (compSType1 > 0) {
166  H1 = new SType1[ n*n ];
167  b1 = new SType1[ n ];
168  //
169  // Construct problem.
170  //
171  ConstructHilbertMatrix< SType1 >(H1, n);
172  ConstructHilbertSumVector< SType1 >(b1, n);
173  //
174  // Try to solve it.
175  //
176  compSType1 = Solve(n, H1, b1, &result1);
177  if (compSType1 < 0 && verbose)
178  std::cout << typeid( result1 ).name() << " -- Cholesky factorization failed (negative diagonal) at row "<<-compSType1<< std::endl;
179  //
180  // Clean up always;
181  delete [] H1; H1 = 0;
182  delete [] b1; b1 = 0;
183  }
184  if (compSType2 > 0) {
185  H2 = new SType2[ n*n ];
186  b2 = new SType2[ n ];
187  //
188  // Construct problem.
189  //
190  ConstructHilbertMatrix< SType2 >(H2, n);
191  ConstructHilbertSumVector< SType2 >(b2, n);
192  //
193  // Try to solve it.
194  //
195  compSType2 = Solve(n, H2, b2, &result2);
196  if (compSType2 < 0 && verbose)
197  std::cout << typeid( result2 ).name() << " -- Cholesky factorization failed (negative diagonal) at row "<<-compSType2<< std::endl;
198  //
199  // Clean up always.
200  delete [] H2; H2 = 0;
201  delete [] b2; b2 = 0;
202  }
203  if (convSType2 > 0) {
204  //
205  // Create and construct the problem in higher precision
206  //
207  if (!H1) H1 = new SType1[ n*n ];
208  if (!b1) b1 = new SType1[ n ];
209  ConstructHilbertMatrix( H1, n );
210  ConstructHilbertSumVector( b1, n );
211  //
212  if (!H2) H2 = new SType2[ n*n ];
213  if (!b2) b2 = new SType2[ n ];
214  //
215  // Convert the problem from SType1 to SType2 ( which should be of lesser precision )
216  //
217  ConvertHilbertMatrix(H1, H2, n);
218  ConvertHilbertSumVector(b1, b2, n);
219  //
220  // Try to solve it.
221  //
222  convSType2 = Solve(n, H2, b2, &result1_2);
223  if (convSType2 < 0 && verbose)
224  std::cout << typeid( result1_2 ).name() << " (converted) -- Cholesky factorization failed (negative diagonal) at row "<<-convSType2<< std::endl;
225  //
226  // Clean up
227  //
228  delete [] H2; H2 = 0;
229  delete [] b2; b2 = 0;
230  delete [] H1; H1 = 0;
231  delete [] b1; b1 = 0;
232  }
233  if (convSType1 > 0) {
234  //
235  // Create and construct the problem in lower precision
236  //
237  if (!H2) H2 = new SType2[ n*n ];
238  if (!b2) b2 = new SType2[ n ];
239  ConstructHilbertMatrix(H2, n);
241  //
242  if (!H1) H1 = new SType1[ n*n ];
243  if (!b1) b1 = new SType1[ n ];
244  //
245  // Convert the problem from SType2 to SType1 ( which should be of higher precision )
246  //
247  ConvertHilbertMatrix(H2, H1, n);
248  ConvertHilbertSumVector(b2, b1, n);
249  //
250  // Try to solve it.
251  //
252  convSType1 = Solve(n, H1, b1, &result2_1);
253  if (convSType1 < 0 && verbose)
254  std::cout << typeid( result2_1 ).name() << " (converted) -- Cholesky factorization failed (negative diagonal) at row "<<-convSType1<< std::endl;
255  //
256  // Clean up
257  //
258  delete [] H1; H1 = 0;
259  delete [] b1; b1 = 0;
260  delete [] H2; H2 = 0;
261  delete [] b2; b2 = 0;
262  }
263  if (verbose && (compSType1>0 || compSType2>0 || convSType1>0 || convSType2>0) ) {
264  std::cout << "***************************************************" << std::endl;
265  std::cout << "Dimension of Hilbert Matrix : "<< n << std::endl;
266  std::cout << "***************************************************" << std::endl;
267  std::cout << "Datatype : Absolute error || x_hat - x ||"<< std::endl;
268  std::cout << "---------------------------------------------------" << std::endl;
269  }
270  if (compSType1>0 && verbose)
271  std::cout << typeid( result1 ).name() << "\t : "<< result1 << std::endl;
272
273  if (convSType1>0 && verbose)
274  std::cout << typeid( result2_1 ).name() <<"(converted) : "<< result2_1 << std::endl;
275
276  if (convSType2>0 && verbose)
277  std::cout << typeid( result1_2 ).name() <<"(converted) : "<< result2_1 << std::endl;
278
279  if (compSType2>0 && verbose)
280  std::cout << typeid( result2 ).name() << "\t : "<< result2 << std::endl;
281  //
282  // Increment counter.
283  //
284  n++;
285  }
286
287 #ifdef HAVE_TEUCHOS_ARPREC
288  mp::mp_finalize();
289 #endif
290
291  return 0;
292 }
293
294 template<typename TYPE>
295 void ConstructHilbertMatrix(TYPE* A, int n) {
296  TYPE scal1 = ScalarTraits<TYPE>::one();
297  for(int i = 0; i < n; i++) {
298  for(int j = 0; j < n; j++) {
299  A[i + (j * n)] = (scal1 / (i + j + scal1));
300  }
301  }
302 }
303
304 template<typename TYPE>
305 void ConstructHilbertSumVector(TYPE* x, int n) {
306  TYPE scal0 = ScalarTraits<TYPE>::zero();
307  TYPE scal1 = ScalarTraits<TYPE>::one();
308  for(int i = 0; i < n; i++) {
309  x[i] = scal0;
310  for(int j = 0; j < n; j++) {
311  x[i] += (scal1 / (i + j + scal1));
312  }
313  }
314 }
315
316 template<typename TYPE1, typename TYPE2>
317 void ConvertHilbertMatrix(TYPE1* A, TYPE2* B, int n) {
318  for(int i = 0; i < n; i++) {
319  for(int j = 0; j < n; j++) {
320  B[i + (j * n)] = A[i + (j * n)];
321  }
322  }
323 }
324
325 template<typename TYPE1, typename TYPE2>
326 void ConvertHilbertSumVector(TYPE1* x, TYPE2* y, int n) {
327  for(int i = 0; i < n; i++) {
328  y[i] = x[i];
329  }
330 }
331
332 #ifdef HAVE_TEUCHOS_ARPREC
333 template<>
334 void ConvertHilbertMatrix(mp_real* A, double* B, int n) {
335  for(int i = 0; i < n; i++) {
336  for(int j = 0; j < n; j++) {
337  B[i + (j * n)] = dble( A[i + (j * n)] );
338  }
339  }
340 }
341
342 template<>
343 void ConvertHilbertSumVector(mp_real* x, double* y, int n) {
344  for(int i = 0; i < n; i++) {
345  y[i] = dble( x[i] );
346  }
347 }
348 #endif
349
350 #ifdef HAVE_TEUCHOS_GNU_MP
351 template<>
352 void ConvertHilbertMatrix(mpf_class* A, double* B, int n) {
353  for(int i = 0; i < n; i++) {
354  for(int j = 0; j < n; j++) {
355  B[i + (j * n)] = A[i + (j * n)].get_d();
356  }
357  }
358 }
359
360 template<>
361 void ConvertHilbertSumVector(mpf_class* x, double* y, int n) {
362  for(int i = 0; i < n; i++) {
363  y[i] = x[i].get_d();
364  }
365 }
366 #endif
367
368 template<typename TYPE>
369 int Cholesky(TYPE* A, int n) {
370  TYPE scal0 = ScalarTraits<TYPE>::zero();
371  for(int k = 0; k < n; k++) {
372  for(int j = k + 1; j < n; j++) {
373  TYPE alpha = A[k + (j * n)] / A[k + (k * n)];
374  for(int i = j; i < n; i++) {
375  A[j + (i * n)] -= (alpha * A[k + (i * n)]);
376  }
377  }
378  if(A[k + (k * n)] <= scal0) {
379  return -k;
380  }
381  TYPE beta = ScalarTraits<TYPE>::squareroot(A[k + (k * n)]);
382  for(int i = k; i < n; i++) {
383  A[k + (i * n)] /= beta;
384  }
385  }
386  return 1;
387 }
388
389 template<typename TYPE>
390 int Solve(int n, TYPE* H, TYPE* b, TYPE* err) {
391  TYPE scal0 = ScalarTraits<TYPE>::zero();
392  TYPE scal1 = ScalarTraits<TYPE>::one();
393  TYPE scalNeg1 = scal0 - scal1;
394  BLAS<int, TYPE> blasObj;
395  TYPE* x = new TYPE[n];
396  for(int i = 0; i < n; i++) {
397  x[i] = scal1;
398  }
399  int choleskySuccessful = Cholesky(H, n);
400  if(choleskySuccessful > 0) {
403  blasObj.AXPY(n, scalNeg1, x, 1, b, 1);
404  *err = blasObj.NRM2(n, b, 1);
405  }
406  delete[] x;
407  return choleskySuccessful;
408 }
409
410 template<typename TYPE>
411 void PrintArrayAsVector(TYPE* x, int n) {
412  std::cout << "[";
413  for(int i = 0; i < n; i++) {
414  std::cout << " " << x[i];
415  }
416  std::cout << " ]" << std::endl;
417 }
418
419 template<typename TYPE>
420 void PrintArrayAsMatrix(TYPE* a, int m, int n) {
421  std::cout << "[";
422  for(int i = 0; i < m; i++) {
423  if(i != 0) {
424  std::cout << " ";
425  }
426  std::cout << "[";
427  for(int j = 0; j < n; j++) {
428  std::cout << " " << a[i + (j * m)];
429  }
430  std::cout << " ]";
431  if(i != (m - 1)) {
432  std::cout << std::endl;
433  }
434  }
435  std::cout << "]" << std::endl;
436 }
437
438 #ifdef HAVE_TEUCHOS_ARPREC
439 template<>
440 void PrintArrayAsVector(mp_real* x, int n) {
441  std::cout << "[ ";
442  for(int i = 0; i < n; i++) {
443  if(i != 0) {
444  std::cout << " ";
445  }
446  std::cout << x[i];
447  }
448  std::cout << "]" << std::endl;
449 }
450
451 template<>
452 void PrintArrayAsMatrix(mp_real* 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  if(j != 0) {
461  std::cout << " ";
462  }
463  std::cout << " " << a[i + (j * m)];
464  }
465  std::cout << " ]";
466  if(i != (m - 1)) {
467  std::cout << std::endl;
468  }
469  }
470  std::cout << "]" << std::endl;
471 }
472 #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.