mtl2lapack.h

MTL C++ Numeric Library

    int              _info;

    /*make_sure tau's size is greater than or equal to min(m,n)*/
    if (int(tau.size()) < MTL_MIN(_m,_n))
      return -105;

    mtl_lapack_dispatch::geqpf(_m, _n, _a, _lda, jpivot.data(), tau.data(),
                               _info);

    a.set_contiguous(_a);

    return _info;
  }

  
  //: QR Factorization of a General Matrix
  //
  //   QR Factorization of a MxN General Matrix A.
  // <UL>
  //   <LI> a       (IN/OUT - matrix(M,N)) On entry, the coefficient matrix A. On exit , the upper triangle and diagonal is the min(M,N) by N upper triangular matrix R.  The lower triangle, together with the tau vector, is the orthogonal matrix Q as a product of min(M,N) elementary reflectors.
  //
  //   <LI> tau     (OUT - vector (min(M,N))) Vector of the same numerical type as A. The scalar factors of the elementary reflectors.
  //
  //   <LI> info    (OUT - <tt>int</tt>)
  //   0   : function completed normally
  //   < 0 : The ith argument, where i = abs(return value) had an illegal value.
  // </UL>
  //!component: function
  //!category: mtl2lapack

  template <class LapackMatA,class VectorT>
  int geqrf (LapackMatA& a,
	     VectorT & tau)
  {
    int              _m = a.nrows();
    int              _n = a.ncols();
    int              _lda = a.minor() ;
    int              _info;
 
    /*make_sure tau's size is greater than or equal to min(m,n)*/
    if (int(tau.size()) < MTL_MIN(_n, _m))
      return -104;

    /*call dispatcher*/
    mtl_lapack_dispatch::geqrf(_m, _n, a.data(), _lda, tau.data(), _info);

    return _info;
  }

 
  //: Solution to a linear system in a general matrix.
  //
  //   Computes the solution to a real system of linear equations A*X=B
  //   (simple driver).  LU decomposition with partial pivoting and row
  //   interchanges is used to solve the system.
  // <UL>
  //   <LI> a       (IN/OUT - matrix(M,N)) On entry, the coefficient matrix A, and the factors L and U from the factorization A = P*L*U on exit.
  //
  //   <LI> ipivot  (OUT - vector(N)) Integer vector. The row i of A was interchanged with row IPIV(i).
  //
  //   <LI> b        (IN/OUT - matrix(ldb,NRHS)) Matrix of same numerical type as A. On entry, the NxNRHS matrix of the right hand side matrix B. On a successful exit, it is the NxNRHS solution matrix X.
  //
  //   <LI> info     (OUT - <tt>int</tt>)
  //   0     : function completed normally
  //   < 0   : The ith argument, where i = abs(return value) had an illegal value.
  //   > 0   : U(i,i), where i = return value, is exactly zero and U is therefore singular. The LU factorization has been completed, but the solution could not be computed.
  // </UL>
  //!component: function
  //!category: mtl2lapack

  template <class LapackMatA, class LapackMatB, class VectorInt>
  int gesv (LapackMatA& a,
	    VectorInt & ipivot,
	    LapackMatB& b)
  {
    int              _lda = a.minor();  
    int              _n = a.ncols();
    int              _nrhs = b.ncols();
 
    int              _ldb = b.nrows(); 
    int              _info;

    /* make sure matrices are same size */
    if (int(a.nrows()) != int(b.nrows()))
      return -100;

    /* make sure ipivot is big enough */
    if (int(ipivot.size()) < _n)
      return -104;

    /* call dispatcher */
    mtl_lapack_dispatch::gesv(_n, _nrhs, a.data(), _lda, ipivot.data(),
			      b.data(), _ldb, _info);

    return _info;
  }

  //: LU factorization of a general matrix A.  
  //    Computes an LU factorization of a general M-by-N matrix A using
  //    partial pivoting with row interchanges. Factorization has the form
  //    A = P*L*U.
  // <UL>
  //   <LI> a       (IN/OUT - matrix(M,N)) On entry, the coefficient matrix A to be factored. On exit, the factors L and U from the factorization A = P*L*U.
  //
  //   <LI> ipivot  (OUT - vector(min(M,N))) Integer vector. The row i of A was interchanged with row IPIV(i).
  //
  //   <LI> info    (OUT - <tt>int</tt>)
  //   0   :  successful exit
  //   < 0 :  If INFO = -i, then the i-th argument had an illegal value.
  //   > 0 :  If INFO = i, then U(i,i) is exactly zero. The  factorization has been completed, but the factor U is exactly singular, and division by zero will occur if it is used to solve a system of equations.
  //
  // </UL>
  //
  //!component: function
  //!category: mtl2lapack
  template <class LapackMatrix, class VectorInt>
  int getrf (LapackMatrix& a,
	     VectorInt& ipivot)
  {
    typename LapackMatrix::value_type* _a = a.data();
    int _lda = a.minor();
    int _n = a.ncols();
    int _m = a.nrows();
    int _info;
    
    mtl_lapack_dispatch::getrf (_m, _n,	_a, _lda, ipivot.data(), _info);

    return _info;
  }

  //:  Solution to a system using LU factorization 
  //   Solves a system of linear equations A*X = B with a general NxN
  //   matrix A using the LU factorization computed by GETRF.
  // <UL>
  //   <LI> transa  (IN - char)  'T' for the transpose of A, 'N' otherwise.
  //
  //   <LI> a       (IN - matrix(M,N)) The factors L and U from the factorization A = P*L*U as computed by GETRF.
  //
  //   <LI> ipivot  (IN - vector(min(M,N))) Integer vector. The pivot indices from GETRF; row i of A was interchanged with row IPIV(i).
  //
  //   <LI> b       (IN/OUT - matrix(ldb,NRHS)) Matrix of same numerical type as A. On entry, the right hand side matrix B. On exit, the solution matrix X.
  //
  //   <LI> info    (OUT - <tt>int</tt>)
  //   0   : function completed normally
  //   < 0 : The ith argument, where i = abs(return value) had an illegal value.
  //   > 0 : if INFO =  i,  U(i,i)  is  exactly  zero;  the  matrix is singular and its inverse could not be computed.
  // </UL>
  //!component: function
  //!category: mtl2lapack

  template <class LapackMatrixA, class LapackMatrixB, class VectorInt>
  int getrs (char transa, LapackMatrixA& a,
	     VectorInt& ipivot, LapackMatrixB& b)
  {
    typename LapackMatrixA::value_type* _a = a.data();
    int _lda = a.minor();
    int a_n = a.nrows();
    int b_n = b.nrows();
    int p_n = ipivot.size();

    typename LapackMatrixB::value_type* _b = b.data();
    int _ldb = b.minor();
    int _nrhs = b.ncols(); /* B's ncols is the # of vectors on rhs */

    if (a_n != b_n) /*Test to see if AX=B has correct dimensions */
      return -101;
    if (p_n < a_n)     /*Check to see if ipivot is big enough */
      return -102;

    int _info;
    mtl_lapack_dispatch::getrs (transa, a_n, _nrhs, _a,	_lda, ipivot.data(), 
				_b, _ldb, _info);
    
    return _info;
  } 


  inline
  void error_message(int ret_val) 
  {
    if (ret_val==0)
      std::cerr << "No error" << std::endl;
    else if (ret_val > 0)
      std::cerr << "Computational error." << std::endl;
    else if (ret_val > -101)
      std::cerr << "Argument #" 
		<< -ret_val 
		<< " of lapack subroutine had an illegal value."
		<< std::endl;
    else
      std::cerr << "Argument #"
		<< -(ret_val+100)
		<< "of mtl_lapack subroutine had an illegal value."
		<< std::endl;
  } 

  //: Equilibrate and reduce condition number.
  //  Compute row and column scaling inteded to equilibrate an M-by-N
  //  matrix A and reduce its condition number.
  // <UL>
  //   <LI> a       (IN - matrix(M,N)) The M-by-N matrix whose equilibration factors are to be computed.
  //
  //   <LI> r       (OUT - vector(M)) Real vector. If INFO = 0 or INFO > M, R contains  the  row  scale factors for A.
  //
  //   <LI> c       (OUT - vector(N)) Real vector. If INFO = 0,  C contains the  column  scale  factors for A.
  //
  //   <LI> row_cond (OUT - Real number) If INFO = 0 or INFO > M, ROWCND contains  the  ratio of the smallest R(i) to the largest R(i).  If ROWCND >= 0.1 and AMAX is neither too large nor too  small, it is not worth scaling by R.
  //
  //   <LI> col_cond (OUT - Real number) If INFO = 0, COLCND contains the ratio of the  smallest C(i) to the largest C(i).  If COLCND >= 0.1, it is not worth scaling by C.
  //
  //   <LI> amax    (OUT - Real number) Absolute value of largest matrix element. If  AMAX is  very  close  to overflow or very close to underflow, the matrix should be scaled.
  //
  // </UL>
  //!component: function
  //!category: mtl2lapack

  template <class LapackMatA, class VectorReal, class Real>
  int geequ(const LapackMatA& a, VectorReal& r, VectorReal& c, 
		   Real& row_cond, Real& col_cond, Real& amax)
  {
    int lda = a.minor();
    int m = a.nrows();
    int n = a.ncols();
    int info;
    mtl_lapack_dispatch::geequ(m, n, a.data(), lda,  r.data(), c.data(),
			       row_cond, col_cond, amax, info);
    return info;
  }


  //:  Compute an LQ factorization.
  //  Compute an LQ factorization of a M-by-N matrix A.
  // <UL>
  //   <LI> a     (IN/OUT - matrix(M,N)) On entry, the M-by-N matrix A.  On  exit, the  elements on and below the diagonal of the array contain the m-by-min(m,n) lower trapezoidal matrix L  (L  is lower  triangular if m <= n); the elements above the diagonal, with the array TAU, represent the  unitary matrix  Q as a product of elementary reflectors.
  //
  //   <LI> tau   (OUT - vector(min(M,N)) The scalar factors of the elementary reflectors.
  // </UL>
  //!component: function
  //!category: mtl2lapack
  template <class LapackMatA, class VectorT>
  int gelqf(LapackMatA& a, VectorT& tau)
  {
    int lda = a.minor();
    int m = a.nrows();
    int n = a.ncols();

    if (int(tau.size()) < MTL_MIN(m,n))
      return -104;

    int info;
    mtl_lapack_dispatch::gelqf(m, n, a.data(), lda, tau.data(), info);
    return info;    
  }


  //: Compute least squares solution using svd for Ax = b
  //  Compute least squares solution using svd for Ax = b
  //  A is m by n, x is length n, b is length m
  template <class LapackMatA, class VectorT>
  int gelss(LapackMatA& a,
	    VectorT& x,
	    VectorT& b,
	    double tol=1.e-6)
  {
    int lda = a.minor();
    int m   = a.nrows();
    int n   = a.ncols();
    int ldb = b.size();

    int info, rank;

    const int max_length = max(b.size(), x.size());

    typedef typename VectorT::value_type T;
    T *rhs_data = new T[max_length];
    
    for (int i = 0; i < b.size(); ++i)
      rhs_data[i] = b[i];

    mtl_lapack_dispatch::gelss(m, n, 1, a.data(), lda,
			       rhs_data, max_length, tol, rank, info);

    for (int i = 0; i < x.size(); ++i)
      x[i] = rhs_data[i];

    delete [] rhs_data;

    return info;    
  }

  /*
  template <class LapackMatA, class VectorT>
  int gelss(LapackMatA& a, VectorT& b, double tol=1.e-6)
  {
    int lda = a.minor();
    int m = a.nrows();
    int n = a.ncols();
    int ldb = b.size();

    int info, rank;

    mtl_lapack_dispatch::gelss(m, n, 1, a.data(), lda, 
			       b.data(), ldb, tol, rank, info);

    return info;    
  }
  */



  //:  Generate a matrix Q with orthonormal rows.
  //  Generate an M-by-N real matrix Q with orthonormal rows.
  // <UL>
  //   <LI> a     (IN/OUT - matrix(M,N) On entry, the i-th row must contain the vector which defines the  elementary  reflector  H(i),  for  i  = 1,2,...,k, as returned by GELQF in the first k rows of its array argument A.  On exit, the M-by-N matrix Q.
  //
  //   <LI> tau   (IN - vector(K)) tau[i] must contain the scalar factor of the elementary reflector H(i), as returned by GELQF.
  // </UL>
  //!component: function
  //!category: mtl2lapack

  template <class LapackMatA, class VectorT>
  int orglq(LapackMatA& a, const VectorT& tau)
  {
    int lda = a.minor();
    int m = a.nrows();
    int n = a.ncols();
    int info;
    mtl_lapack_dispatch::orglq(m, n, tau.size(), a.data(), lda,
			       tau.data(), info);
    return info;
  }


//: Generate a matrix Q with orthonormal columns.
  //  Generate an M-by-N real matrix Q with orthonormal columns.
  // <UL>
  //  <LI> a     (IN/OUT - matrix(M,N) On  entry,  the  i-th column must contain the vector which defines the elementary reflector H(i), for i = 1,2,...,k,  as  returned  by GEQRF  in the first k columns of its array argument A.  On  exit, the  M-by-N matrix Q.
  //
  //  <LI> tau   (IN - vector(K)) tau[i] must contain the scalar factor of the elementary reflector H(i), as returned by GEQRF.
  // </UL>
  //
  //!component: function
  //!category: mtl2lapack

  template <class LapackMatA, class VectorT>
  int orgqr(LapackMatA& a, const VectorT& tau)
  {
    int lda = a.minor();
    int m = a.nrows();
    int n = a.ncols();
    int info;
    mtl_lapack_dispatch::orgqr(m, n, tau.size(), a.data(), lda,
			       tau.data(), info);
    return info;
  }



  template <class LapackMatA, class VectorT>
  int gesvd(char jobu, char jobvt, LapackMatA& a, VectorT& s, 
	    LapackMatU& u, LapackMatVT& vt)
  {
    int info;
    int lda = a.minor();
    int m = a.nrows();
    int n = a.ncols();
    int ldu = u.minor();
    int ldvt = vt.minor();

    mtl_lapack_dispatch::gesvd(jobu, jobv, m, n, a.data(), lda, s.data(),
			       u.data(), ldu, vt.data(), ldvt, info);
    return info;
  }

} /* namespace mtl2lapack */

#endif