Teuchos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
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
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
90using 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
102template<typename TYPE>
103void ConstructHilbertMatrix(TYPE*, int);
104
105template<typename TYPE>
106void ConstructHilbertSumVector(TYPE*, int);
107
108template<typename TYPE1, typename TYPE2>
109void ConvertHilbertMatrix(TYPE1*, TYPE2*, int);
110
111template<typename TYPE1, typename TYPE2>
112void ConvertHilbertSumVector(TYPE1*, TYPE2*, int);
113
114#ifdef HAVE_TEUCHOS_ARPREC
115template<>
116void ConvertHilbertMatrix(mp_real*, double*, int);
117
118template<>
119void ConvertHilbertSumVector(mp_real*, double*, int);
120#endif
121
122#ifdef HAVE_TEUCHOS_GNU_MP
123template<>
124void ConvertHilbertMatrix(mpf_class*, double*, int);
125
126template<>
127void ConvertHilbertSumVector(mpf_class*, double*, int);
128#endif
129
130template<typename TYPE>
131int Cholesky(TYPE*, int);
132
133template<typename TYPE>
134int Solve(int, TYPE*, TYPE*, TYPE*);
135
136template<typename TYPE>
137void PrintArrayAsVector(TYPE*, int);
138
139template<typename TYPE>
140void PrintArrayAsMatrix(TYPE*, int, int);
141
142#ifdef HAVE_TEUCHOS_ARPREC
143template<>
144void PrintArrayAsVector(mp_real*, int);
145
146template<>
147void PrintArrayAsMatrix(mp_real*, int, int);
148#endif
149
150int 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 );
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 ];
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
326template<typename TYPE>
327void 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
336template<typename TYPE>
337void 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
348template<typename TYPE1, typename TYPE2>
349void 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
357template<typename TYPE1, typename TYPE2>
358void 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
365template<>
366void 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
374template<>
375void 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
383template<>
384void 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
392template<>
393void 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
400template<typename TYPE>
401int 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
421template<typename TYPE>
422int 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
442template<typename TYPE>
443void 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
451template<typename TYPE>
452void 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
471template<>
472void 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
483template<>
484void 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
Templated interface class to BLAS routines.
Basic command line parser for input from (argc,argv[])
Teuchos header file which uses auto-configuration information to include necessary C++ headers.
Templated BLAS wrapper.
Class that helps parse command line input arguments from (argc,argv[]) and set options.
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.
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 <- y+alpha*x.
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,...
ScalarTraits< ScalarType >::magnitudeType NRM2(const OrdinalType &n, const ScalarType *x, const OrdinalType &incx) const
Compute the 2-norm of the vector x.
int main()
Definition: evilMain.cpp:75
Definition: PackageA.cpp:3
Definition: PackageB.cpp:3
std::string Teuchos_Version()
void ConstructHilbertMatrix(TYPE *, int)
void PrintArrayAsVector(TYPE *, int)
int Solve(int, TYPE *, TYPE *, TYPE *)
int Cholesky(TYPE *, int)
void PrintArrayAsMatrix(TYPE *, int, int)
void ConvertHilbertSumVector(TYPE1 *, TYPE2 *, int)
void ConvertHilbertMatrix(TYPE1 *, TYPE2 *, int)
void ConstructHilbertSumVector(TYPE *, int)
static T one()
Returns representation of one for this scalar type.
static T zero()
Returns representation of zero for this scalar type.
static T squareroot(T x)
Returns a number of magnitudeType that is the square root of this scalar type x.