dposv_if.cc

Go to the documentation of this file.

#define LINEAR_DPOSV_IF_CC_VERSION \
  "TF_DPOSV_IF_CC   V1.1"

#include <iostream>
#include <aff/fortranshape.h>
#include <linearxx/f77lapack.h>
#include <linearxx/lapackxx.h>
#include <linearxx/error.h>

namespace linear {

namespace lapack {

  void dposv(char UPLO, TDmatrix& A, TDmatrix& B, int& INFO, const bool& debug)
  {
    LINEAR_debug(debug, "dposv", "solve system of linear equations");
    // convert to Fortran interfacing array
    typedef aff::FortranArray<TDmatrix> Tfortranarray;
    Tfortranarray FA(A);
    Tfortranarray FB(B);
    LINEAR_assert((FA.last(0)==FA.last(1)),
                  "ERROR (dposv): Matrix A is not a square matrix");
    LINEAR_assert((FB.last(0)==FA.last(0)),
                  "ERROR (dposv): Dimension of right hand side "
                  "does not match");
    f77lapack::integer N=FA.last(0);
    f77lapack::integer NRHS=FB.last(1);
    f77lapack::integer LDA=FA.dimlast(0);
    f77lapack::integer LDB=FB.dimlast(0);
    f77lapack::integer f77INFO;
    f77lapack::doublereal* PA=FA.castedpointer<f77lapack::doublereal>();
    f77lapack::doublereal* PB=FB.castedpointer<f77lapack::doublereal>();
    LINEAR_debug(debug, "dposv", "N: " << N 
               << " NRHS: " << NRHS
               << " LDA: " << LDA
               << " LDB: " << LDB);
    LINEAR_debug(debug, "dposv", "first in A: " << *PA 
               << " first in B: " << *PB);
    f77lapack::dposv_(&UPLO, &N, &NRHS, PA, &LDA, PB, &LDB, &f77INFO);
    LINEAR_debug(debug, "dposv", "INFO: " << f77INFO);
    INFO=f77INFO;
    LINEAR_LAPACK_failure((INFO<0),"dposv",
                          "the -INFO-th argument had an illegal value",
                          INFO);
    LINEAR_LAPACK_failure((INFO>0),"dposv",
                          "the leading INFO-th order is not positive definite",
                          INFO);
  } // dposv

  TDmatrix dposv(const TDmatrix::Tcoc& A, 
                 const TDmatrix::Tcoc& B, 
                 char UPLO, const bool& debug)
  {
    // factorized Matrix
    TDmatrix M(A.copyout());
    // return value
    TDmatrix X=B.copyout();
    int INFO;
    dposv(UPLO, M, X, INFO, debug);
    return(X);
  } // dposv
  
} // namespace lapack

} // namespace linear

/* ----- END OF dposv_if.cc ----- */