mtl2lapack.h

MTL C++ Numeric Library

  {
    int ldwork = 2*m;
    double* work = new double[ldwork];
    MTL_FCALL(dgelqf)(m, n, da, lda, dtau, work, ldwork, info);
    delete [] work;
  }

  inline void gelqf(const int& m, const int& n,
		    float da[], const int& lda,
		    float dtau[], int& info)
  {
    int ldwork = 2*m;
    float* work = new float[ldwork];
    MTL_FCALL(sgelqf)(m, n, da, lda, dtau, work, ldwork, info);
    delete [] work;
  }

  inline void gelqf(const int& m, const int& n,
		    std::complex<float> da[], const int& lda,
		    std::complex<float> dtau[], int& info)
  {
    int ldwork = 2*m;
    std::complex<float>* work = new std::complex<float>[ldwork];
    MTL_FCALL(cgelqf)(m, n, da, lda, dtau, work, ldwork, info);
    delete [] work;
  }

  inline void gelqf(const int& m, const int& n,
		    std::complex<double> da[], const int& lda,
		    std::complex<double> dtau[], int& info)
  {
    int ldwork = 2*m;
    std::complex<double>* work = new std::complex<double>[ldwork];
    MTL_FCALL(zgelqf)(m, n, da, lda, work, dtau, ldwork, info);
    delete [] work;
  }


  inline void gelss(const int& m, const int& n, const int& nrhs,
		    double da[], const int& lda, double b[], const int& ldb,
		    const double& rcond, int& rank, int& info)
  {
    int ldwork = 5*max(m,n);
    double* work = new double[ldwork];
    double* s = new double[ldwork];
    MTL_FCALL(dgelss)(m, n, nrhs, da, lda, b, ldb, s, rcond, rank, work, ldwork, info);
    delete [] work;
    delete [] s;
  }

  inline void gelss(const int& m, const int& n, const int& nrhs,
		    float da[], const int& lda, float b[], const int& ldb,
		    const float& rcond, int& rank, int& info)
  {
    int ldwork = 5*max(m,n);
    float* work = new float[ldwork];
    float* s = new float[ldwork];
    MTL_FCALL(sgelss)(m, n, nrhs, da, lda, b, ldb, s, rcond, rank, work, ldwork, info);
    delete [] work;
    delete [] s;
  }


  inline void gelss(const int& m, const int& n, const int& nrhs,
		    std::complex<float> da[], const int& lda,
		    std::complex<float> b[], const int& ldb,
		    const float& rcond, int& rank, int& info)
  {
    int ldwork = 5*max(m,n);
    std::complex<float>* work = new std::complex<float>[ldwork];
    std::complex<float>* s = new std::complex<float>[ldwork];
    MTL_FCALL(cgelss)(m, n, nrhs, da, lda, b, ldb, s, rcond, rank, work, ldwork, info);
    delete [] work;
    delete [] s;
  }


  inline void gelss(const int& m, const int& n, const int& nrhs,
		    std::complex<double> da[], const int& lda, 
		    std::complex<double> b[], const int& ldb,
		    const double& rcond, int& rank, int& info)
  {
    int ldwork = 5*max(m,n);
    std::complex<double>* work = new std::complex<double>[ldwork];
    std::complex<double>* s = new std::complex<double>[ldwork];
    MTL_FCALL(zgelss)(m, n, nrhs, da, lda, b, ldb, s, rcond, rank, work, ldwork, info);
    delete [] work;
    delete [] s;
  }


  inline void orglq(const int& m, const int& n, const int& k,
		    double da[], int& lda, double dtau[], int& info)
  {
    int ldwork = 2*m;
    double* work = new double[ldwork];
    MTL_FCALL(dorglq)(m, n, k, da, lda, dtau, work, ldwork, info);
    delete [] work;
  }

  inline void orglq(const int& m, const int& n, const int& k,
		    float da[], int& lda, float dtau[], int& info)
  {
    int ldwork = 2*m;
    float* work = new float[ldwork];
    MTL_FCALL(sorglq)(m, n, k, da, lda, dtau, work, ldwork, info);
    delete [] work;
  }

  inline void orgqr(const int& m, const int& n, const int& k,
		    double da[], int& lda, double dtau[], int& info)
  {
    int ldwork = 2*m;
    double* work = new double[ldwork];
    MTL_FCALL(dorgqr)(m, n, k, da, lda, dtau, work, ldwork, info);
    delete [] work;
  }

  inline void orgqr(const int& m, const int& n, const int& k,
		    float da[], int& lda, float dtau[], int& info)
  {
    int ldwork = 2*m;
    float* work = new float[ldwork];
    MTL_FCALL(sorgqr)(m, n, k, da, lda, dtau, work, ldwork, info);
    delete [] work;
  }

  inline void gesvd(const char& jobu, const char& jobvt, const int& m,
		    const int& n, double da[], const int& lda, double ds[],
		    double du[], const int& ldu, double dvt[], 
		    const int& ldvt, int& info)
  {
    //allocate work space
    int lwork = max(3*min(m,n)+max(m,n),5*min(m,n));
    double* work = new double[lwork];

    MTL_FCALL(dgesvd)(jobu, jobvt, m, n, da, lda, ds, du, ldu, dvt, ldvt, work, lwork, info);

    delete [] work;
  }


  inline void gesvd(const char& jobu, const char& jobvt, const int& m,
		    const int& n, float da[], const int& lda, float ds[],
		    float du[], const int& ldu, float dvt[], 
		    const int& ldvt, int& info)
  {
    //allocate work space
    int lwork = max(3*min(m,n)+max(m,n),5*min(m,n));
    float* work = new float[lwork];

    MTL_FCALL(sgesvd)(jobu, jobvt, m, n, da, lda, ds, du, ldu, dvt, ldvt, work, lwork, info);

    delete [] work;
  }

  inline void gesvd(const char& jobu, const char& jobvt, const int& m,
		    const int& n, std::complex<double> da[], const int& lda,
		    double ds[],
		    std::complex<double> du[], const int& ldu, 
		    std::complex<double> dvt[], const int& ldvt, int& info)
  {
    //allocate work space
    int lwork = 2*min(m,n)+max(m,n);
    std::complex<double>* work = new std::complex<double>[lwork];
    int lrwork = 5*min(m,n);
    double* rwork = new double[lrwork];

    MTL_FCALL(zgesvd)(jobu, jobvt, m, n, da, lda, ds, du, ldu, dvt, ldvt, work, lwork, rwork, info);

    delete [] work;
    delete [] rwork;
  }

  inline void gesvd(const char& jobu, const char& jobvt, const int& m,
		    const int& n, std::complex<float> da[], const int& lda,
		    float ds[],
		    std::complex<float> du[], const int& ldu, 
		    std::complex<float> dvt[], const int& ldvt, int& info)
  {
    //allocate work space
    int lwork = 2*min(m,n)+max(m,n);
    std::complex<float>* work = new std::complex<float>[lwork];
    int lrwork = 5*min(m,n);
    float* rwork = new float[lrwork];

    MTL_FCALL(cgesvd)(jobu, jobvt, m, n, da, lda, ds, du, ldu, dvt, ldvt, work, lwork, rwork, info);

    delete [] work;
    delete [] rwork;
  }


} /* namespace mtl lapack dispatch */

namespace mtl2lapack {

  using namespace mtl;

  /* MTL to Lapack Interface*/


  //: Lapack Matrix
  //
  // Use this matrix type constructor to create the type of matrix to
  // use in conjunction with the mtl2lapack functions.
  //
  // <p>The vector type you use with mtl2lapack functions must be
  // contiguous in memory, and have a function data() defined which
  // returns a pointer to that memory, and a function size() with
  // gives the length.
  //
  //!tparam: T - the matrix element type, either float, double, std::complex< float >, or std::complex< double >
  //!tparam: External - Memory managed by MTL? - internal
  //!category: mtl2lapack, generators
  //!component: type
  //!example: getrf.cc, geequ.cc, gecon.cc, geev.cc, 
  template <class T, int External=mtl::internal>
  struct lapack_matrix {
    typedef typename matrix<T, rectangle<>, 
                            dense<External>, column_major>::type type;
  };


  //: Estimate the reciprocal of the condition number of a general matrix.
  //
  //   Estimate the reciprocal of the condition number of a general matrix
  //   A, with either the one-norm or the infinity-norm. GECON uses the LU
  //   factorization computed by GETRF. Currently MTL only handles column
  //   oriented matrices for this function.
  //  
  // <UL>
  // <LI>norm    (IN - <tt>char</tt>) Specifies whether to do one-norm or infinity norm. One of the following is required: 
  //  
  //   '1' the 1-norm condition number 
  //   'I' the infinity-norm condition number
  //  
  //   <LI> a       (IN - matrix(M,N)) The coefficient matrix A.
  //  
  //   <LI> anorm   (IN - Real number) The one or infinity norm of matrix A, which is of the same numerical type as A.
  //
  //  
  //   <LI> rcond  (OUT - Real number) The estimated reciprocal condition of matrix A.
  //  
  //
  //   <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 Real>
  int gecon(char  _norm,
	    const LapackMatA& a,
	    const Real& _anorm, 
	    Real& rcond)
  {
    int             _lda = a.minor();
    int             _n = a.nrows();
    int             _info;
    
    mtl_lapack_dispatch::gecon(_norm, _n, a.data(), _lda,
			       _anorm, rcond, _info);

    return _info;
  }


enum GEEV_JOBV {GEEV_CALC_LEFT,	GEEV_CALC_RIGHT, GEEV_CALC_BOTH, GEEV_CALC_NONE};


//: Compute the eigenvalues.
//   Compute for an N-by-N non-symmetric matrix A, the eigenvalues, and,  optionally, the left and/or right eigenvectors (simple driver).
// <UL>
//   <LI> jobv    (IN - <tt>int</tt>) Specifies which eigenvectors to solve for (left, right, both, or neither). One of four values:
//   GEEV_CALC_LEFT - function computes left eigenvectors.
//   GEEV_CALC_RIGHT - function computes right eigenvectors.
//   GEEV_CALC_BOTH - function computes both left and right eigenvectors.
//   GEEV_CALC_NONE - function computes neither left nor right eigenvectors.
//
//   <LI> a       (IN/OUT - matrix(M,N)) The coefficient matrix A on entry, and is overwritten on exit.
//
//   <LI> w       (OUT - vector(N))  Std::Complex Vector with same base precision as A. The computed real and imaginary parts of the eigenvectors.
//
//   <LI> vl      (OUT - matrix(N,N)) Matrix of same numerical type as A. Used to store left eigenvectors if they are computed.
//
//   <LI> vr      (OUT - matrix(N,N)) Matrix of same numerical type as A. Used to store right eigenvectors if they are computed.
//
//   <LI> info    (OUT - <tt>int</tt>)
//   0   : function completed normally
//   < 0 : The ith argument, where i = abs(return value) had an illegal value.
//   > 0 : The QR algorithm failed to compute all the eigenvalues and no eigenvectors have been computed.
// </UL>
//
  //!component: function
  //!category: mtl2lapack

  template <class LapackMatA, class LapackMatVL, 
            class LapackMatVR, class VectorComplex>
  int geev(int jobv,
	   LapackMatA& a,        /* N x N */
	   VectorComplex& w,     /* N x N */
	   LapackMatVL& vl,  /* LDVL x N */
	   LapackMatVR& vr)  /* LDVR x N */
  {
    char _jobvr, _jobvl;

    int  _info;
 
    /*determine _jobvr and _jobvl */
    switch (jobv) {
    case 0: /* GEEV_CALC_LEFT: */
      _jobvr = 'N';
      _jobvl = 'V';
      break;
    case 1: /* GEEV_CALC_RIGHT: */
      _jobvr = 'V';
      _jobvl = 'N';
      break;
    case 2:  /* GEEV_CALC_BOTH: */
      _jobvr = 'V';
      _jobvl = 'V';
      break;
    case 3:  /* GEEV_CALC_NONE:*/
      _jobvr = 'N';
      _jobvl = 'N';
      break;
    default:
      assert(0);
    }

    int             _lda = a.minor();
    int             a_n = a.major();
    int             _ldvl = vl.minor();
    int             _ldvr = vr.minor();
  
    /* make sure matrix is square */
    if (a.nrows() != a.ncols())
      return -100;

    /* make sure result eigenvalue vector is long enough */
    if (int(w.size()) != a_n)
      return -104;

    if (_jobvl == 'V')
      if (int(vl.nrows()) < a_n || int(vl.ncols()) < a_n)
        return -105;

    if (_jobvr == 'V')
      if (int(vr.nrows()) < a_n || int(vr.ncols()) < a_n )
	return -106;

    /* call dispatcher */
    mtl_lapack_dispatch::geev(_jobvl, _jobvr, a_n, a.data(), _lda, w.data(),
                              vl.data(), _ldvl, vr.data(), _ldvr, _info);

    return _info;

  }
  

  //:  QR Factorization with Column Pivoting.
  //  
  //   QR Factorization with Column Pivoting of a MxN General Matrix A.
  // <UL>
  //   <LI> a       (IN/OUT - matrix(M,N)) On entry, the coefficient matrix A. On exit, its upper triangle 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> jpivot  (IN/OUT - vector(N)) Integer vector. On entry, if JPVT(i) != 0, the i-th column of A is permuted to the front of A*P. If JPVT(i) = 0, the i-th column of A is a free column. On exit, if jpivot(i) = k, then the i-th column of A*P was the k-th column of A.
  //  
  //   <LI> tau     (OUT - vector (min(M,N))) Vector of 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 VectorInt, class VectorT>
  int geqpf (LapackMatA & a,
	     VectorInt & jpivot,
	     VectorT & tau)
  {
    typename LapackMatA::value_type * _a = a.get_contiguous();
    int              _m = a.nrows();
    int              _n = a.ncols();
    int              _lda = a.minor();