numerics/example/hilbert/cxx_main.cpp

Go to the documentation of this file.

// @HEADER
// *****************************************************************************
//                    Teuchos: Common Tools Package
//
// Copyright 2004 NTESS and the Teuchos contributors.
// SPDX-License-Identifier: BSD-3-Clause
// *****************************************************************************
// @HEADER

// Teuchos Example: Hilbert

// This example showcases the usage of BLAS generics with an arbitrary precision
// library -- ARPREC.

// Hilbert matrices are classical examples of ill-conditioned matrices. Cholesky
// factorization fails on higher-order Hilbert matrices because they lose their
// positive definiteness when represented with floating-point numbers. We have
// attempted to alleviate this problem with arbitrary precision.

// The example program compares two datatypes, scalar type 1 and scalar type 2,
// which can be customized below using #defines. Default types are mp_real (from
// ARPREC) and double. The mp_real datatype must be initialized with a maximum
// precision value, also customizable below. (Default is 32.)

// For a given size n, an n-by-n Hilbert matrix H and a n-by-1 std::vector b are
// constructed such that, if Hx* = b, the true solution x* is a one-std::vector.
// Cholesky factorization is attempted on H; if it fails, no further tests are
// attempted for that datatype. If it is successful, the approximate solution x~
// is computed with a pair of BLAS TRSM (triangular solve) calls. Then, the
// two-norm of (x* - x~) is computed with BLAS AXPY (std::vector update) and BLAS
// NRM2. The program output is of the form:

//     [size of Hilbert matrix]: [two-norm of (x* - x~)]

// Tests for scalar type 2 are performed before scalar type 1 because scalar
// type 2 fails at Cholesky factorization for much lower values of n if the
// mp_real precision is sufficiently large.

// Timing analysis still remains to be done for this example, which should be
// easily accomplished with the timing mechanisms native to Teuchos.

#include "Teuchos_CommandLineProcessor.hpp"
#include "Teuchos_ConfigDefs.hpp"
#include "Teuchos_BLAS.hpp"
#include "Teuchos_Version.hpp"
#include <typeinfo>

#ifdef HAVE_TEUCHOS_ARPREC
#include <arprec/mp_real.h>
#endif

#ifdef HAVE_TEUCHOS_GNU_MP
#include "gmp.h"
#include "gmpxx.h"
#endif


using namespace Teuchos;

#ifdef HAVE_TEUCHOS_ARPREC
#define SType1     mp_real
#elif defined(HAVE_TEUCHOS_GNU_MP)
#define SType1     mpf_class
#else
#define SType1     double
#endif
#define SType2     double
#define OType    int

template<typename TYPE>
void ConstructHilbertMatrix(TYPE*, int);

template<typename TYPE>
void ConstructHilbertSumVector(TYPE*, int);

template<typename TYPE1, typename TYPE2>
void ConvertHilbertMatrix(TYPE1*, TYPE2*, int);

template<typename TYPE1, typename TYPE2>
void ConvertHilbertSumVector(TYPE1*, TYPE2*, int);

#ifdef HAVE_TEUCHOS_ARPREC
template<>
void ConvertHilbertMatrix(mp_real*, double*, int);

template<>
void ConvertHilbertSumVector(mp_real*, double*, int);
#endif

#ifdef HAVE_TEUCHOS_GNU_MP
template<>
void ConvertHilbertMatrix(mpf_class*, double*, int);

template<>
void ConvertHilbertSumVector(mpf_class*, double*, int);
#endif

template<typename TYPE>
int Cholesky(TYPE*, int);

template<typename TYPE>
int Solve(int, TYPE*, TYPE*, TYPE*);

template<typename TYPE>
void PrintArrayAsVector(TYPE*, int);

template<typename TYPE>
void PrintArrayAsMatrix(TYPE*, int, int);

#ifdef HAVE_TEUCHOS_ARPREC
template<>
void PrintArrayAsVector(mp_real*, int);

template<>
void PrintArrayAsMatrix(mp_real*, int, int);
#endif

int main(int argc, char *argv[]) {

  std::cout << Teuchos::Teuchos_Version() << std::endl << std::endl;
  //
  // Create command line processor.
  //
  Teuchos::CommandLineProcessor hilbertCLP(true, false);
  //
  // Set option for precision and verbosity
  int precision = 32;
  hilbertCLP.setOption("precision", &precision, "Arbitrary precision");
  bool verbose = false;
  hilbertCLP.setOption("verbose", "quiet", &verbose, "Verbosity of example");
  //
  // Parse command line.
  hilbertCLP.parse( argc, argv );

#ifdef HAVE_TEUCHOS_ARPREC
  mp::mp_init( precision );
#endif

#ifdef HAVE_TEUCHOS_GNU_MP
  mpf_set_default_prec( precision );
  std::cout<< "The precision of the GNU MP variable is (in bits) : "<< mpf_get_default_prec() << std::endl;
#endif
  //
  // Keep track of valid datatypes
  //
  int compSType1 = 1;  // Perform cholesky factorization of matrices of SType1
  int compSType2 = 1;  // Perform cholesky factorization of matrices of SType2
  int convSType1 = 1;  // Perform cholesky factorization of matrices of SType1 (that were converted from SType2)
  int convSType2 = 1;  // Perform cholesky factorization of matrices of SType2 (that were converted from SType1)

  int n = 2;  // Initial dimension of hilbert matrix.
  //
  // Error in solution.
  //
  SType1 result1, result2_1;
  SType2 result2, result1_2;
  //
  // Create pointers to necessary matrices/vectors.
  //
  SType1 *H1=0, *b1=0;
  SType2 *H2=0, *b2=0;
  //
  while ( compSType1>0 || compSType2>0 || convSType1>0 || convSType2>0 ) {

    if (compSType1 > 0) {
      H1 = new SType1[ n*n ];
      b1 = new SType1[ n ];
      //
      // Construct problem.
      //
      ConstructHilbertMatrix< SType1 >(H1, n);
      ConstructHilbertSumVector< SType1 >(b1, n);
      //
      // Try to solve it.
      //
      compSType1 = Solve(n, H1, b1, &result1);
      if (compSType1 < 0 && verbose)
  std::cout << typeid( result1 ).name() << " -- Cholesky factorization failed (negative diagonal) at row "<<-compSType1<< std::endl;
      //
      // Clean up always;
      delete [] H1; H1 = 0;
      delete [] b1; b1 = 0;
    }
    if (compSType2 > 0) {
      H2 = new SType2[ n*n ];
      b2 = new SType2[ n ];
      //
      // Construct problem.
      //
      ConstructHilbertMatrix< SType2 >(H2, n);
      ConstructHilbertSumVector< SType2 >(b2, n);
      //
      // Try to solve it.
      //
      compSType2 = Solve(n, H2, b2, &result2);
      if (compSType2 < 0 && verbose)
  std::cout << typeid( result2 ).name() << " -- Cholesky factorization failed (negative diagonal) at row "<<-compSType2<< std::endl;
      //
      // Clean up always.
      delete [] H2; H2 = 0;
      delete [] b2; b2 = 0;
    }
    if (convSType2 > 0) {
      //
      // Create and construct the problem in higher precision
      //
      if (!H1) H1 = new SType1[ n*n ];
      if (!b1) b1 = new SType1[ n ];
      ConstructHilbertMatrix( H1, n );
      ConstructHilbertSumVector( b1, n );
      //
      if (!H2) H2 = new SType2[ n*n ];
      if (!b2) b2 = new SType2[ n ];
      //
      // Convert the problem from SType1 to SType2 ( which should be of lesser precision )
      //
      ConvertHilbertMatrix(H1, H2, n);
      ConvertHilbertSumVector(b1, b2, n);
      //
      // Try to solve it.
      //
      convSType2 = Solve(n, H2, b2, &result1_2);
      if (convSType2 < 0 && verbose)
  std::cout << typeid( result1_2 ).name() << " (converted) -- Cholesky factorization failed (negative diagonal) at row "<<-convSType2<< std::endl;
      //
      // Clean up
      //
      delete [] H2; H2 = 0;
      delete [] b2; b2 = 0;
      delete [] H1; H1 = 0;
      delete [] b1; b1 = 0;
    }
    if (convSType1 > 0) {
      //
      // Create and construct the problem in lower precision
      //
      if (!H2) H2 = new SType2[ n*n ];
      if (!b2) b2 = new SType2[ n ];
      ConstructHilbertMatrix(H2, n);
      ConstructHilbertSumVector(b2, n);
      //
      if (!H1) H1 = new SType1[ n*n ];
      if (!b1) b1 = new SType1[ n ];
      //
      // Convert the problem from SType2 to SType1 ( which should be of higher precision )
      //
      ConvertHilbertMatrix(H2, H1, n);
      ConvertHilbertSumVector(b2, b1, n);
      //
      // Try to solve it.
      //
      convSType1 = Solve(n, H1, b1, &result2_1);
      if (convSType1 < 0 && verbose)
  std::cout << typeid( result2_1 ).name() << " (converted) -- Cholesky factorization failed (negative diagonal) at row "<<-convSType1<< std::endl;
      //
      // Clean up
      //
      delete [] H1; H1 = 0;
      delete [] b1; b1 = 0;
      delete [] H2; H2 = 0;
      delete [] b2; b2 = 0;
    }
    if (verbose && (compSType1>0 || compSType2>0 || convSType1>0 || convSType2>0) ) {
      std::cout << "***************************************************" << std::endl;
      std::cout << "Dimension of Hilbert Matrix : "<< n << std::endl;
      std::cout << "***************************************************" << std::endl;
      std::cout << "Datatype : Absolute error || x_hat - x ||"<< std::endl;
      std::cout << "---------------------------------------------------" << std::endl;
    }
    if (compSType1>0 && verbose)
      std::cout << typeid( result1 ).name() << "\t : "<< result1 << std::endl;

    if (convSType1>0 && verbose)
      std::cout << typeid( result2_1 ).name() <<"(converted) : "<< result2_1 << std::endl;

    if (convSType2>0 && verbose)
      std::cout << typeid( result1_2 ).name() <<"(converted) : "<< result2_1 << std::endl;

    if (compSType2>0 && verbose)
      std::cout << typeid( result2 ).name() << "\t : "<< result2 << std::endl;
    //
    // Increment counter.
    //
    n++;
  }

#ifdef HAVE_TEUCHOS_ARPREC
  mp::mp_finalize();
#endif

  return 0;
}

template<typename TYPE>
void ConstructHilbertMatrix(TYPE* A, int n) {
  TYPE scal1 = ScalarTraits<TYPE>::one();
  for(int i = 0; i < n; i++) {
    for(int j = 0; j < n; j++) {
      A[i + (j * n)] = (scal1 / (i + j + scal1));
    }
  }
}

template<typename TYPE>
void ConstructHilbertSumVector(TYPE* x, int n) {
  TYPE scal0 = ScalarTraits<TYPE>::zero();
  TYPE scal1 = ScalarTraits<TYPE>::one();
  for(int i = 0; i < n; i++) {
    x[i] = scal0;
    for(int j = 0; j < n; j++) {
      x[i] += (scal1 / (i + j + scal1));
    }
  }
}

template<typename TYPE1, typename TYPE2>
void ConvertHilbertMatrix(TYPE1* A, TYPE2* B, int n) {
  for(int i = 0; i < n; i++) {
    for(int j = 0; j < n; j++) {
      B[i + (j * n)] = A[i + (j * n)];
    }
  }
}

template<typename TYPE1, typename TYPE2>
void ConvertHilbertSumVector(TYPE1* x, TYPE2* y, int n) {
  for(int i = 0; i < n; i++) {
    y[i] = x[i];
  }
}

#ifdef HAVE_TEUCHOS_ARPREC
template<>
void ConvertHilbertMatrix(mp_real* A, double* B, int n) {
  for(int i = 0; i < n; i++) {
    for(int j = 0; j < n; j++) {
      B[i + (j * n)] = dble( A[i + (j * n)] );
    }
  }
}

template<>
void ConvertHilbertSumVector(mp_real* x, double* y, int n) {
  for(int i = 0; i < n; i++) {
    y[i] = dble( x[i] );
  }
}
#endif

#ifdef HAVE_TEUCHOS_GNU_MP
template<>
void ConvertHilbertMatrix(mpf_class* A, double* B, int n) {
  for(int i = 0; i < n; i++) {
    for(int j = 0; j < n; j++) {
      B[i + (j * n)] = A[i + (j * n)].get_d();
    }
  }
}

template<>
void ConvertHilbertSumVector(mpf_class* x, double* y, int n) {
  for(int i = 0; i < n; i++) {
    y[i] = x[i].get_d();
  }
}
#endif

template<typename TYPE>
int Cholesky(TYPE* A, int n) {
  TYPE scal0 = ScalarTraits<TYPE>::zero();
  for(int k = 0; k < n; k++) {
    for(int j = k + 1; j < n; j++) {
      TYPE alpha = A[k + (j * n)] / A[k + (k * n)];
      for(int i = j; i < n; i++) {
  A[j + (i * n)] -= (alpha * A[k + (i * n)]);
      }
    }
    if(A[k + (k * n)] <= scal0) {
      return -k;
    }
    TYPE beta = ScalarTraits<TYPE>::squareroot(A[k + (k * n)]);
    for(int i = k; i < n; i++) {
      A[k + (i * n)] /= beta;
    }
  }
  return 1;
}

template<typename TYPE>
int Solve(int n, TYPE* H, TYPE* b, TYPE* err) {
  TYPE scal0 = ScalarTraits<TYPE>::zero();
  TYPE scal1 = ScalarTraits<TYPE>::one();
  TYPE scalNeg1 = scal0 - scal1;
  BLAS<int, TYPE> blasObj;
  TYPE* x = new TYPE[n];
  for(int i = 0; i < n; i++) {
    x[i] = scal1;
  }
  int choleskySuccessful = Cholesky(H, n);
  if(choleskySuccessful > 0) {
    blasObj.TRSM(Teuchos::LEFT_SIDE, Teuchos::UPPER_TRI, Teuchos::TRANS, Teuchos::NON_UNIT_DIAG, n, 1, scal1, H, n, b, n);
    blasObj.TRSM(Teuchos::LEFT_SIDE, Teuchos::UPPER_TRI, Teuchos::NO_TRANS, Teuchos::NON_UNIT_DIAG, n, 1, scal1, H, n, b, n);
    blasObj.AXPY(n, scalNeg1, x, 1, b, 1);
    *err = blasObj.NRM2(n, b, 1);
  }
  delete[] x;
  return choleskySuccessful;
}

template<typename TYPE>
void PrintArrayAsVector(TYPE* x, int n) {
  std::cout << "[";
  for(int i = 0; i < n; i++) {
    std::cout << " " << x[i];
  }
  std::cout << " ]" << std::endl;
}

template<typename TYPE>
void PrintArrayAsMatrix(TYPE* a, int m, int n) {
  std::cout << "[";
  for(int i = 0; i < m; i++) {
    if(i != 0) {
      std::cout << " ";
    }
    std::cout << "[";
    for(int j = 0; j < n; j++) {
      std::cout << " " << a[i + (j * m)];
    }
    std::cout << " ]";
    if(i != (m - 1)) {
      std::cout << std::endl;
    }
  }
  std::cout << "]" << std::endl;
}

#ifdef HAVE_TEUCHOS_ARPREC
template<>
void PrintArrayAsVector(mp_real* x, int n) {
  std::cout << "[ ";
  for(int i = 0; i < n; i++) {
    if(i != 0) {
      std::cout << "  ";
    }
    std::cout << x[i];
  }
  std::cout << "]" << std::endl;
}

template<>
void PrintArrayAsMatrix(mp_real* a, int m, int n) {
  std::cout << "[";
  for(int i = 0; i < m; i++) {
    if(i != 0) {
      std::cout << " ";
    }
    std::cout << "[";
    for(int j = 0; j < n; j++) {
      if(j != 0) {
  std::cout << "  ";
      }
      std::cout << " " << a[i + (j * m)];
    }
    std::cout << " ]";
    if(i != (m - 1)) {
      std::cout << std::endl;
    }
  }
  std::cout << "]" << std::endl;
}
#endif