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 = MTL_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, class LapackMatU, class LapackMatVT>  int gesvd(const char& jobu, const char& jobvt, LapackMatA& a, VectorT& s, 	    LapackMatU& u, LapackMatVT& vt)  {    int info;    const int lda = a.minor();    const int m = a.nrows();    const int n = a.ncols();    const int ldu = u.minor();    const int ldvt = vt.minor();    mtl_lapack_dispatch::gesvd(jobu, jobvt, m, n, a.data(), lda, s.data(),			       u.data(), ldu, vt.data(), ldvt, info);    return info;  }} /* namespace mtl2lapack */#endif