--- a/Eigen/src/Householder/Householder.h +++ a/Eigen/src/Householder/Householder.h @@ -87,16 +87,126 @@ void MatrixBase::makeHouseholde beta = sqrt(numext::abs2(c0) + tailSqNorm); if (numext::real(c0)>=RealScalar(0)) beta = -beta; essential = tail / (c0 - beta); tau = conj((beta - c0) / beta); } } +/** Computes the stable elementary reflector H in a stable way such that: + * \f$ H *this = [ beta 0 ... 0]^T \f$ + * where the transformation H is: + * \f$ H = I - tau v v^*\f$ + * and the vector v is: + * \f$ v^T = [1 essential^T] \f$ + * + * \note This routine provides, not only a stable elementary reflector, + * but a compatible reflector with Lapack routines. + * For instance here, beta is ensured to be greater than zero. + * + * On output: + * \param essential the essential part of the vector \c v + * \param tau the scaling factor of the Householder transformation + * \param beta ( > 0) the result of H * \c *this + * + * \sa MatrixBase::makeHouseholder(), MatrixBase::makeHouseholderInPlace(), + * MatrixBase::applyHouseholderOnTheLeft(), MatrixBase::applyHouseholderOnTheRight() + */ +template +void MatrixBase::makeStableHouseholderInPlace(Scalar& tau, RealScalar& beta) +{ + VectorBlock::ret> essentialPart(derived(), 1, size()-1); + makeStableHouseholder(essentialPart, tau, beta); +} + + +/** Computes the stable elementary reflector H in a stable way such that: + * \f$ H *this = [ beta 0 ... 0]^T \f$ + * where the transformation H is: + * \f$ H = I - tau v v^*\f$ + * and the vector v is: + * \f$ v^T = [1 essential^T] \f$ + * + * \note This routine provides, not only a stable elementary reflector, + * but a compatible reflector with Lapack routines. + * For instance here, beta is ensured to be greater than zero. + * + * On output: + * \param essential the essential part of the vector \c v + * \param tau the scaling factor of the Householder transformation + * \param beta ( > 0) the result of H * \c *this + * + * \sa MatrixBase::makeHouseholder(), MatrixBase::makeHouseholderInPlace(), + * MatrixBase::applyHouseholderOnTheLeft(), MatrixBase::applyHouseholderOnTheRight() + */ +template +template +void MatrixBase::makeStableHouseholder( + EssentialPart& essential, + Scalar& tau, + RealScalar& beta) const +{ + using std::sqrt; + using numext::conj; + + EIGEN_STATIC_ASSERT_VECTOR_ONLY(EssentialPart) + VectorBlock tail(derived(), 1, size()-1); + RealScalar tailNorm = size()==1 ? RealScalar(0) : tail.blueNorm(); + Scalar c0 = coeff(0); + if(tailNorm == RealScalar(0) && numext::imag(c0)==RealScalar(0)) + { + beta = numext::real(c0); + if(beta >= 0) + tau = RealScalar(0); + else + { + tau = 2; + essential.setZero(); + beta = -beta; + } + } + else + { + // Compute beta = sqrt(c0*c0 + tailNorm * tailNorm) in a stable way + RealScalar c0_abs = (std::abs)(c0); + RealScalar max_abs = (std::max)(c0_abs, tailNorm); + RealScalar min_abs = (std::min)(c0_abs, tailNorm); + if(min_abs == 0) beta = max_abs; + else beta = max_abs * sqrt(1. + (min_abs / max_abs)*(min_abs / max_abs)); + if (numext::real(c0) < RealScalar(0)) beta = -beta; + + if(beta < 0) + { + beta = -beta; + tau = numext::conj(beta-c0)/beta; + essential = tail / (c0-beta); + } + else + { + RealScalar gamma = numext::real(c0) + beta; + RealScalar delta = (numext::imag(c0)*(numext::imag(c0)/gamma) + tailNorm * (tailNorm/gamma)); + + if(NumTraits::IsComplex) + { +// std::complex delta_out(delta, numext::imag(c0)); + Scalar delta_out = c0; + numext::real_ref(delta_out) = delta; + tau = - delta_out/beta; + essential = tail/(delta_out); + } + else + { + tau = delta / beta; + essential = tail / (-delta); + } + } + + } +} /** Apply the elementary reflector H given by * \f$ H = I - tau v v^*\f$ * with * \f$ v^T = [1 essential^T] \f$ * from the left to a vector or matrix. * * On input: * \param essential the essential part of the vector \c v --- a/Eigen/src/QR/ColPivHouseholderQR.h +++ a/Eigen/src/QR/ColPivHouseholderQR.h @@ -407,112 +407,127 @@ typename MatrixType::RealScalar ColPivHo template typename MatrixType::RealScalar ColPivHouseholderQR::logAbsDeterminant() const { eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!"); return m_qr.diagonal().cwiseAbs().array().log().sum(); } +template +void colpivhouseholder_qr_inplace(MatrixType& qr, HCoeffsType& hCoeffs, RealRowVectorType& colSqNorms, + IntRowVectorType& colsTranspositions, PermutationType& colsPermutation, + RowVectorType& temp, typename MatrixType::RealScalar& maxpivot, + typename MatrixType::Index& nonzero_pivots, typename MatrixType::Index& det_pq) +{ + typedef typename MatrixType::Index Index; + typedef typename MatrixType::Scalar Scalar; + typedef typename MatrixType::RealScalar RealScalar; + typedef typename PermutationType::Index PermIndexType; + using std::abs; + Index rows = qr.rows(); + Index cols = qr.cols(); + Index size = qr.diagonalSize(); + + // the column permutation is stored as int indices, so just to be sure: + eigen_assert(cols<=NumTraits::highest()); + + hCoeffs.resize(size); + + temp.resize(cols); //FIXME Why is this not a local variable + + colsTranspositions.resize(qr.cols()); + Index number_of_transpositions = 0; + + colSqNorms.resize(cols); + for(Index k = 0; k < cols; ++k) + colSqNorms.coeffRef(k) = qr.col(k).squaredNorm(); + + RealScalar threshold_helper = colSqNorms.maxCoeff() * numext::abs2(NumTraits::epsilon()) / RealScalar(rows); + + nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case) + maxpivot = RealScalar(0); + + for(Index k = 0; k < size; ++k) + { + // first, we look up in our table colSqNorms which column has the biggest squared norm + Index biggest_col_index; + RealScalar biggest_col_sq_norm = colSqNorms.tail(cols-k).maxCoeff(&biggest_col_index); + biggest_col_index += k; + + // since our table colSqNorms accumulates imprecision at every step, we must now recompute + // the actual squared norm of the selected column. + // Note that not doing so does result in solve() sometimes returning inf/nan values + // when running the unit test with 1000 repetitions. + biggest_col_sq_norm = qr.col(biggest_col_index).tail(rows-k).squaredNorm(); + + // we store that back into our table: it can't hurt to correct our table. + colSqNorms.coeffRef(biggest_col_index) = biggest_col_sq_norm; + + // if the current biggest column is smaller than epsilon times the initial biggest column, + // terminate to avoid generating nan/inf values. + // Note that here, if we test instead for "biggest == 0", we get a failure every 1000 (or so) + // repetitions of the unit test, with the result of solve() filled with large values of the order + // of 1/(size*epsilon). + if(biggest_col_sq_norm < threshold_helper * RealScalar(rows-k)) + { + nonzero_pivots = k; + hCoeffs.tail(size-k).setZero(); + qr.bottomRightCorner(rows-k,cols-k) + .template triangularView() + .setZero(); + break; + } + + // apply the transposition to the columns + colsTranspositions.coeffRef(k) = biggest_col_index; + if(k != biggest_col_index) { + qr.col(k).swap(qr.col(biggest_col_index)); + std::swap(colSqNorms.coeffRef(k), colSqNorms.coeffRef(biggest_col_index)); + ++number_of_transpositions; + } + + // generate the householder vector, store it below the diagonal + RealScalar beta; + qr.col(k).tail(rows-k).makeHouseholderInPlace(hCoeffs.coeffRef(k), beta); +// qr.col(k).tail(rows-k).makeStableHouseholderInPlace(hCoeffs.coeffRef(k), beta); + + // apply the householder transformation to the diagonal coefficient + qr.coeffRef(k,k) = beta; + + // remember the maximum absolute value of diagonal coefficients + if(abs(beta) > maxpivot) maxpivot = abs(beta); + + // apply the householder transformation + qr.bottomRightCorner(rows-k, cols-k-1) + .applyHouseholderOnTheLeft(qr.col(k).tail(rows-k-1), hCoeffs.coeffRef(k), &temp.coeffRef(k+1)); + + // update our table of squared norms of the columns + colSqNorms.tail(cols-k-1) -= qr.row(k).tail(cols-k-1).cwiseAbs2(); + } + + colsPermutation.setIdentity(PermIndexType(cols)); + for(PermIndexType k = 0; k < nonzero_pivots; ++k) + colsPermutation.applyTranspositionOnTheRight(k, PermIndexType(colsTranspositions.coeff(k))); + + det_pq = (number_of_transpositions%2) ? -1 : 1; +} /** Performs the QR factorization of the given matrix \a matrix. The result of * the factorization is stored into \c *this, and a reference to \c *this * is returned. * * \sa class ColPivHouseholderQR, ColPivHouseholderQR(const MatrixType&) */ template ColPivHouseholderQR& ColPivHouseholderQR::compute(const MatrixType& matrix) { - using std::abs; - Index rows = matrix.rows(); - Index cols = matrix.cols(); - Index size = matrix.diagonalSize(); + m_qr = matrix; + colpivhouseholder_qr_inplace(m_qr, m_hCoeffs, m_colSqNorms, m_colsTranspositions, m_colsPermutation, + m_temp, m_maxpivot, m_nonzero_pivots, m_det_pq); - // the column permutation is stored as int indices, so just to be sure: - eigen_assert(cols<=NumTraits::highest()); - - m_qr = matrix; - m_hCoeffs.resize(size); - - m_temp.resize(cols); - - m_colsTranspositions.resize(matrix.cols()); - Index number_of_transpositions = 0; - - m_colSqNorms.resize(cols); - for(Index k = 0; k < cols; ++k) - m_colSqNorms.coeffRef(k) = m_qr.col(k).squaredNorm(); - - RealScalar threshold_helper = m_colSqNorms.maxCoeff() * numext::abs2(NumTraits::epsilon()) / RealScalar(rows); - - m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case) - m_maxpivot = RealScalar(0); - - for(Index k = 0; k < size; ++k) - { - // first, we look up in our table m_colSqNorms which column has the biggest squared norm - Index biggest_col_index; - RealScalar biggest_col_sq_norm = m_colSqNorms.tail(cols-k).maxCoeff(&biggest_col_index); - biggest_col_index += k; - - // since our table m_colSqNorms accumulates imprecision at every step, we must now recompute - // the actual squared norm of the selected column. - // Note that not doing so does result in solve() sometimes returning inf/nan values - // when running the unit test with 1000 repetitions. - biggest_col_sq_norm = m_qr.col(biggest_col_index).tail(rows-k).squaredNorm(); - - // we store that back into our table: it can't hurt to correct our table. - m_colSqNorms.coeffRef(biggest_col_index) = biggest_col_sq_norm; - - // if the current biggest column is smaller than epsilon times the initial biggest column, - // terminate to avoid generating nan/inf values. - // Note that here, if we test instead for "biggest == 0", we get a failure every 1000 (or so) - // repetitions of the unit test, with the result of solve() filled with large values of the order - // of 1/(size*epsilon). - if(biggest_col_sq_norm < threshold_helper * RealScalar(rows-k)) - { - m_nonzero_pivots = k; - m_hCoeffs.tail(size-k).setZero(); - m_qr.bottomRightCorner(rows-k,cols-k) - .template triangularView() - .setZero(); - break; - } - - // apply the transposition to the columns - m_colsTranspositions.coeffRef(k) = biggest_col_index; - if(k != biggest_col_index) { - m_qr.col(k).swap(m_qr.col(biggest_col_index)); - std::swap(m_colSqNorms.coeffRef(k), m_colSqNorms.coeffRef(biggest_col_index)); - ++number_of_transpositions; - } - - // generate the householder vector, store it below the diagonal - RealScalar beta; - m_qr.col(k).tail(rows-k).makeHouseholderInPlace(m_hCoeffs.coeffRef(k), beta); - - // apply the householder transformation to the diagonal coefficient - m_qr.coeffRef(k,k) = beta; - - // remember the maximum absolute value of diagonal coefficients - if(abs(beta) > m_maxpivot) m_maxpivot = abs(beta); - - // apply the householder transformation - m_qr.bottomRightCorner(rows-k, cols-k-1) - .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), m_hCoeffs.coeffRef(k), &m_temp.coeffRef(k+1)); - - // update our table of squared norms of the columns - m_colSqNorms.tail(cols-k-1) -= m_qr.row(k).tail(cols-k-1).cwiseAbs2(); - } - - m_colsPermutation.setIdentity(PermIndexType(cols)); - for(PermIndexType k = 0; k < m_nonzero_pivots; ++k) - m_colsPermutation.applyTranspositionOnTheRight(k, PermIndexType(m_colsTranspositions.coeff(k))); - - m_det_pq = (number_of_transpositions%2) ? -1 : 1; m_isInitialized = true; return *this; } namespace internal { template --- a/lapack/complex_double.cpp +++ a/lapack/complex_double.cpp @@ -10,8 +10,9 @@ #define SCALAR std::complex #define SCALAR_SUFFIX z #define SCALAR_SUFFIX_UP "Z" #define REAL_SCALAR_SUFFIX d #define ISCOMPLEX 1 #include "cholesky.cpp" #include "lu.cpp" +#include "qr.cpp" --- a/lapack/complex_single.cpp +++ a/lapack/complex_single.cpp @@ -10,8 +10,9 @@ #define SCALAR std::complex #define SCALAR_SUFFIX c #define SCALAR_SUFFIX_UP "C" #define REAL_SCALAR_SUFFIX s #define ISCOMPLEX 1 #include "cholesky.cpp" #include "lu.cpp" +#include "qr.cpp" --- a/lapack/double.cpp +++ a/lapack/double.cpp @@ -10,8 +10,9 @@ #define SCALAR double #define SCALAR_SUFFIX d #define SCALAR_SUFFIX_UP "D" #define ISCOMPLEX 0 #include "cholesky.cpp" #include "lu.cpp" #include "eigenvalues.cpp" +#include "qr.cpp" --- a/lapack/eigenvalues.cpp +++ a/lapack/eigenvalues.cpp @@ -1,21 +1,22 @@ // This file is part of Eigen, a lightweight C++ template library // for linear algebra. // // Copyright (C) 2011 Gael Guennebaud +// Copyright (C) 2013 Désiré Nuentsa-Wakam // // This Source Code Form is subject to the terms of the Mozilla // Public License v. 2.0. If a copy of the MPL was not distributed // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. #include "common.h" #include -// computes an LU factorization of a general M-by-N matrix A using partial pivoting with row interchanges +// computes all eigenvalues and, optionally, eigenvectors of a real symmetric matrix A EIGEN_LAPACK_FUNC(syev,(char *jobz, char *uplo, int* n, Scalar* a, int *lda, Scalar* w, Scalar* /*work*/, int* lwork, int *info)) { // TODO exploit the work buffer bool query_size = *lwork==-1; *info = 0; if(*jobz!='N' && *jobz!='V') *info = -1; else if(UPLO(*uplo)==INVALID) *info = -2; @@ -72,8 +73,64 @@ EIGEN_LAPACK_FUNC(syev,(char *jobz, char } vector(w,*n) = eig.eigenvalues(); if(computeVectors) matrix(a,*n,*n,*lda) = eig.eigenvectors(); return 0; } + +// computes all the eigenvalues, and optionally, the eigenvectors of a real generalized symmetric-definite eigenproblem, of the form A*x=(lambda)*B*x, A*Bx=(lambda)*x, or B*A*x=(lambda)*x + +EIGEN_LAPACK_FUNC(sygv, (int *itype, char *jobz, char *uplo, int *n, Scalar *a, int *lda, Scalar *b, int *ldb, Scalar *w, Scalar *work, int *lwork, int *info)) +{ + bool query_size = (*lwork == -1); + *info = 0; + + if (*itype < 1 || *itype > 3) *info = -1; + else if (!(*jobz=='V' || *jobz == 'N')) *info = -2; + else if (UPLO(*uplo)==INVALID) *info = -3; + else if (*n < 0) *info = -4; + else if (*lda == (std::max)(1, *n)) *info = -6; + else if (*ldb == (std::max)(1, *n)) *info = -8; + + // TODO Estimate workspace and return if asked by query_size + + if (*info != 0) + { + int e = -*info; + return xerbla_(SCALAR_SUFFIX_UP"SYGV", &e, 5); + } + + if(*n == 0) return 0; + + // Call the GeneralizedSelfAdjointEigenSolver class + PlainMatrixType matA(*n, *n), matB(*n, *n); + int option; + if (*jobz=='V' || *jobz=='v') option = ComputeEigenvectors; + if (*itype == 1) option |= Ax_lBx; + else if (*itype == 2) option |= ABx_lx; + else option |= BAx_lx; + + if(UPLO(*uplo)==UP) matA = matrix(a,*n,*n,*lda).adjoint(); + else matA = matrix(a,*n,*n,*lda); + + if(UPLO(*uplo)==UP) matB = matrix(b,*n,*n,*ldb).adjoint(); + else matB = matrix(b,*n,*n,*ldb); + + GeneralizedSelfAdjointEigenSolver es(matA, matB, option); + if(es.info() == NoConvergence) + { + vector(w,*n).setZero(); + if((option&ComputeEigenvectors) == ComputeEigenvectors) + matrix(a, *n,*n, *lda).setIdentity(); + //FIXME Set the correct return code in *info + return 0; + } + + // Get the eigenvalues + vector(w, *n) = es.eigenvalues(); + //Get the eigenvectors + matrix(a, *n, *n, *lda) = es.eigenvectors(); + + return 0; +} --- a/lapack/lu.cpp +++ a/lapack/lu.cpp @@ -32,16 +32,17 @@ EIGEN_LAPACK_FUNC(getrf,(int *m, int *n, ::blocked_lu(*m, *n, a, *lda, ipiv, nb_transpositions)); for(int i=0; i=0) *info = ret+1; + return 0; } //GETRS solves a system of linear equations // A * X = B or A' * X = B // with a general N-by-N matrix A using the LU factorization computed by GETRF EIGEN_LAPACK_FUNC(getrs,(char *trans, int *n, int *nrhs, RealScalar *pa, int *lda, int *ipiv, RealScalar *pb, int *ldb, int *info)) { --- a/lapack/qr.cpp +++ a/lapack/qr.cpp @@ -0,0 +1,227 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2013 Desire Nuentsa +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/ + +#include "lapack_common.h" +#include + +// computes a QR factorization of a general M-by-N matrix A +EIGEN_LAPACK_FUNC(geqrf, (int *m, int *n, Scalar *a, int *lda, Scalar *tau, Scalar *work, + int *lwork, int *info)) +{ + bool query_size = (*lwork == -1); + *info = 0; + + if(*m < 0) *info = -1; + else if(*n < 0) *info = -2; + else if(*lda < (std::max)(1 + , *m)) *info = -4; + else if((*lwork < (std::max)(1, *n)) && !query_size) *info = -7; + if (*info != 0) + { + int e = -*info; + return xerbla_(SCALAR_SUFFIX_UP"GEQRF", &e, 6); + } +// if (query_size) +// { +// *lwork = 0; +// return 0; +// } + if (*m == 0 || *n == 0) + { + work[0] = 1; + return 0; + } + MatrixType mat(a, *m, *n, OuterStride<>(*lda)); + CompactVectorType hcoeffs(tau, (std::min)(*m,*n)); + householder_qr_inplace_blocked(mat, hcoeffs, 32, work); + hcoeffs = hcoeffs.conjugate(); + return 0; +} + + /* Compute the QR factorization with column pivoting */ + /* FIXME xgeqpf is deprecated in LAPACK, We should define xgeqp3 instead + * However, the implementation of column pivoting QR in Eigen is more similar to xgeqpf than xgeqp3, in terms of performance (BLAS3 operations) + */ + EIGEN_LAPACK_FUNC(geqpf, (int *m, int *n, Scalar *a, int *lda, int *jpvt, Scalar *tau, + Scalar *work, int *info)) + { + *info = 0; + if(*m < 0) *info = -1; + else if(*n < 0) *info = -2; + else if(*lda < (std::max)(1, *m)) *info = -4; + + if (*info != 0) + { + int e = -*info; + return xerbla_(SCALAR_SUFFIX_UP"GEQPF", &e, 6); + } + if((std::min)(*m, *n) == 0) return 0; + + MatrixType qr(a, *m, *n, OuterStride<>(*lda)); + CompactVectorType hcoeffs(tau, (std::min)(*m, *n)); + + Matrix colSqNorms; + Matrix colsTranspositions; + CompactVectorType temp(work, *n); + RealScalar maxpivot; + typename MatrixType::Index nonzero_pivots, det_pq; + + //FIXME Find a way to deal with 'free' columns specified in jpvt. + PermutationMatrix colsPerm(*n); + + // Call ColPivHouseholder kernel + colpivhouseholder_qr_inplace(qr, hcoeffs, colSqNorms, colsTranspositions, + colsPerm, temp, maxpivot, nonzero_pivots, det_pq); + hcoeffs = hcoeffs.conjugate(); + for(int i = 0; i < *n; i++) jpvt[i] = colsPerm.indices()(i)+1; + + return 0; + } + template + int apply_reflectors(const char *method, char *side, char *trans, int *m, int *n, int *k, + T *a, int *lda, T *tau, T *c, int *ldc, T *work, int *lwork, int *info) + { + + *info = 0; + bool query_size = (*lwork==-1); + int nq; // order of Q + int nw; // minimum dimension of Work + if(*side == 'L') + { + nq = *m; + nw = *n; + } else { + nq = *n; + nw = *m; + } + if( (*side != 'L') && (*side != 'R') ) *info = -1; + else if( (*trans != 'N') && (*trans != 'T') && (*trans != 'C') ) *info = -2; + else if (*m < 0) *info = -3; + else if (*n < 0) *info = -4; + else if (*k < 0 || *k > nq) *info = -5; + else if (*lda < (std::max)(1, nq)) *info = -7; + else if (*ldc < (std::max)(1, *m)) *info = -10; + else if (*lwork < (std::max)(1, nw) && !query_size) *info = -12; + + if(*info != 0) + { + int e = -*info; + return xerbla_(method, &e, 6); + } + if(query_size) return 0; + if ( *m == 0 || *n == 0 || *k == 0) + { + work[0] = 1; + return 0; + } + + //Map input arrays to Eigen matrices + //NOTE Here, k should be the the number of nonzero pivots as returned by xgeqpf + int row = (*side == 'L') ? *m : *n; + Map, 0, OuterStride<> > qr(a, row, *k, OuterStride<>(*lda)); + Map, 0, OuterStride<> > dst(c, *m, *n, OuterStride<>(*ldc)); + Map > hcoeffs(tau, *k); + + Matrix matrixQ(*m, *n); + matrixQ.setIdentity(); + matrixQ = householderSequence(qr, hcoeffs); + + // Apply the Householder reflectors + if (*side == 'L') + { + if(*trans == 'N') dst = matrixQ * dst; + //dst.applyOnTheLeft(householderSequence(qr, hcoeffs)); + else if(*trans =='T') dst = matrixQ.transpose() * dst; + //dst.applyOnTheLeft(householderSequence(qr, hcoeffs).transpose()); + else dst = matrixQ.adjoint() * dst; + //dst.applyOnTheLeft(householderSequence(qr, hcoeffs.conjugate())/*.adjoint()*/); + } + else + { + if(*trans == 'N') dst = dst * matrixQ; + //dst.applyOnTheRight(householderSequence(qr, hcoeffs)); + else if(*trans =='T') dst = dst * matrixQ.transpose(); + //dst.applyOnTheRight(householderSequence(qr, hcoeffs).transpose()); + else dst = dst * matrixQ.adjoint(); + //dst.applyOnTheRight(householderSequence(qr, hcoeffs.conjugate())/*.adjoint()*/); + } + return 0; + } + + /* Apply the product of k elementary reflectors to a real matrix + */ + EIGEN_LAPACK_FUNC(ormqr, (char *side, char *trans, int *m, int *n, int *k, RealScalar *a, + int *lda, RealScalar *tau, RealScalar *c, int *ldc, RealScalar *work, + int *lwork, int *info)) + { + int err = apply_reflectors(SCALAR_SUFFIX_UP"ORMQR", side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info); + return err; + } + + /* Apply the product of k elementary reflectors to a complex matrix */ + EIGEN_LAPACK_FUNC(unmqr, (char *side, char *trans, int *m, int *n, int *k, Complex *a, + int *lda, Complex *tau, Complex *c, int *ldc, Complex *work, + int *lwork, int *info)) + { + int err = apply_reflectors(SCALAR_SUFFIX_UP"UNMQR", side, trans, m, n, k, a, lda, tau, c, ldc, + work, lwork, info); + return err; + } + +template +int generate_Q(const char *method, int *m, int *n, int *k, T *a, int *lda, T *tau, + T *work, int *lwork, int *info) +{ + bool query_size = (*lwork == -1); + *info = 0; + if (*m < 0) *info = -1; + else if (*n < 0 || *n > *m) *info = -2; + else if (*k < 0 || *k > *n) *info = -3; + else if (*lda < (std::max)(1, *m)) *info = -5; + else if (*lwork < (std::max)(1, *n) && !query_size) *info = -8; + if (*info != 0) + { + int e = -*info; + return xerbla_(method, &e, 6); + } + if (query_size) + { + work[0] = 0; + return 0; + } + if (*n <= 0) + { + work[0] = 1; + return 0; + } + + Map, 0, OuterStride<> > qr(a, *m, *n, OuterStride<>(*lda)); + Matrix tA(*m, *k); + tA = matrix(a, *m, *k, *lda); + Map > hcoeffs(tau, *k); + qr.setIdentity(); + qr = householderSequence(tA, hcoeffs) * qr; + + return 0; +} +/* Generates a m-by-n real matrix Q from the k elementary reflectors */ +EIGEN_LAPACK_FUNC(orgqr, (int *m, int *n, int *k, RealScalar *a, int *lda, RealScalar *tau, + RealScalar *work, int *lwork, int *info)) +{ + int err = generate_Q(SCALAR_SUFFIX_UP"ORGQR", m, n, k, a, lda, tau, work, lwork, info); + return err; +} + +/* Generates a m-by-n complex matrix Q from the k elementary reflectors */ +EIGEN_LAPACK_FUNC(ungqr, (int *m, int *n, int *k, Complex *a, int *lda, Complex *tau, + Complex *work, int *lwork, int *info)) +{ + int err = generate_Q(SCALAR_SUFFIX_UP"UNGQR", m, n, k, a, lda, tau, work, lwork, info); + return err; +} --- a/lapack/single.cpp +++ a/lapack/single.cpp @@ -10,8 +10,9 @@ #define SCALAR float #define SCALAR_SUFFIX s #define SCALAR_SUFFIX_UP "S" #define ISCOMPLEX 0 #include "cholesky.cpp" #include "lu.cpp" #include "eigenvalues.cpp" +#include "qr.cpp"