From 91b3b3aaab19fc11db18d95a28c1a0be9ae9d9cd Mon Sep 17 00:00:00 2001 From: Desire NUENTSA Date: Fri, 11 Jan 2013 17:16:14 +0100 Subject: [PATCH] Add a sparse QR factorization and update the elimination tree in SparseLU --- Eigen/SparseQR | 28 + Eigen/src/MetisSupport/MetisSupport.h | 3 +- Eigen/src/SparseLU/SparseLU.h | 14 +- Eigen/src/SparseLU/SparseLUBase.h | 4 - Eigen/src/SparseLU/SparseLU_Coletree.h | 53 +- .../src/SparseLU/SparseLU_heap_relax_snode.h | 2 +- Eigen/src/SparseQR/CMakeLists.txt | 6 + Eigen/src/SparseQR/SparseQR.h | 481 ++++++++++++++++++ test/CMakeLists.txt | 1 + test/sparseqr.cpp | 62 +++ 10 files changed, 618 insertions(+), 36 deletions(-) create mode 100644 Eigen/SparseQR create mode 100644 Eigen/src/SparseQR/CMakeLists.txt create mode 100644 Eigen/src/SparseQR/SparseQR.h create mode 100644 test/sparseqr.cpp diff --git a/Eigen/SparseQR b/Eigen/SparseQR new file mode 100644 index 000000000..d2382666a --- /dev/null +++ b/Eigen/SparseQR @@ -0,0 +1,28 @@ +#ifndef EIGEN_SPARSEQR_MODULE_H +#define EIGEN_SPARSEQR_MODULE_H + +#include "SparseCore" +#include "OrderingMethods" +#include "src/Core/util/DisableStupidWarnings.h" + +/** \defgroup SparseQR_Module SparseQR module + * + * + * This module provides a simplicial version of the left-looking Sparse QR decomposition. + * The columns of the input matrix should be reordered to limit the fill-in during the + * decomposition. Built-in methods (COLAMD, AMD) or external methods (METIS) can be used to this end. + * See \link OrderingMethods_Module OrderingMethods_Module \endlink for the list + * of built-in and external ordering methods. + * + * \code + * #include + * \endcode + * + * + */ + +#include "src/SparseQR/SparseQR.h" + +#include "src/Core/util/ReenableStupidWarnings.h" + +#endif \ No newline at end of file diff --git a/Eigen/src/MetisSupport/MetisSupport.h b/Eigen/src/MetisSupport/MetisSupport.h index a762d96f6..3a723b384 100644 --- a/Eigen/src/MetisSupport/MetisSupport.h +++ b/Eigen/src/MetisSupport/MetisSupport.h @@ -122,10 +122,9 @@ public: //NOTE: If Ap is the permuted matrix then perm and iperm vectors are defined as follows // Row (column) i of Ap is the perm(i) row(column) of A, and row (column) i of A is the iperm(i) row(column) of Ap - // To be consistent with the use of the permutation in SparseLU module, we thus keep the iperm vector matperm.resize(m); for (int j = 0; j < m; j++) - matperm.indices()(j) = iperm(j); + matperm.indices()(iperm(j)) = j; } diff --git a/Eigen/src/SparseLU/SparseLU.h b/Eigen/src/SparseLU/SparseLU.h index 2f515eb41..6003237fe 100644 --- a/Eigen/src/SparseLU/SparseLU.h +++ b/Eigen/src/SparseLU/SparseLU.h @@ -170,7 +170,7 @@ class SparseLU /** \brief Reports whether previous computation was successful. * * \returns \c Success if computation was succesful, - * \c NumericalIssue if the PaStiX reports a problem + * \c NumericalIssue if the LU factorization reports a problem, zero diagonal for instance * \c InvalidInput if the input matrix is invalid * * \sa iparm() @@ -320,15 +320,14 @@ void SparseLU::analyzePattern(const MatrixType& mat) } // Compute the column elimination tree of the permuted matrix - /*if (m_etree.size() == 0) */m_etree.resize(m_mat.cols()); - - SparseLUBase::LU_sp_coletree(m_mat, m_etree); + IndexVector firstRowElt; + internal::coletree(m_mat, m_etree,firstRowElt); // In symmetric mode, do not do postorder here if (!m_symmetricmode) { IndexVector post, iwork; // Post order etree - SparseLUBase::LU_TreePostorder(m_mat.cols(), m_etree, post); + internal::treePostorder(m_mat.cols(), m_etree, post); // Renumber etree in postorder @@ -445,13 +444,12 @@ void SparseLU::factorize(const MatrixType& matrix) // Work on one 'panel' at a time. A panel is one of the following : // (a) a relaxed supernode at the bottom of the etree, or // (b) panel_size contiguous columns, defined by the user - int jcol,kcol; + int jcol; IndexVector panel_histo(n); - Index nextu, nextlu, jsupno, fsupc, new_next; Index pivrow; // Pivotal row number in the original row matrix int nseg1; // Number of segments in U-column above panel row jcol int nseg; // Number of segments in each U-column - int irep, icol; + int irep; int i, k, jj; for (jcol = 0; jcol < n; ) { diff --git a/Eigen/src/SparseLU/SparseLUBase.h b/Eigen/src/SparseLU/SparseLUBase.h index 20338a646..893728b5e 100644 --- a/Eigen/src/SparseLU/SparseLUBase.h +++ b/Eigen/src/SparseLU/SparseLUBase.h @@ -22,10 +22,6 @@ struct SparseLUBase typedef LU_GlobalLU_t GlobalLU_t; typedef SparseMatrix MatrixType; - static int etree_find (int i, IndexVector& pp); - static int LU_sp_coletree(const MatrixType& mat, IndexVector& parent); - static void LU_nr_etdfs (int n, IndexVector& parent, IndexVector& first_kid, IndexVector& next_kid, IndexVector& post, int postnum); - static void LU_TreePostorder(int n, IndexVector& parent, IndexVector& post); template static int expand(VectorType& vec, int& length, int nbElts, int keep_prev, int& num_expansions); static int LUMemInit(int m, int n, int annz, int lwork, int fillratio, int panel_size, GlobalLU_t& glu); diff --git a/Eigen/src/SparseLU/SparseLU_Coletree.h b/Eigen/src/SparseLU/SparseLU_Coletree.h index d3bc36ea4..e373525ad 100644 --- a/Eigen/src/SparseLU/SparseLU_Coletree.h +++ b/Eigen/src/SparseLU/SparseLU_Coletree.h @@ -30,9 +30,11 @@ */ #ifndef SPARSELU_COLETREE_H #define SPARSELU_COLETREE_H + +namespace internal { /** Find the root of the tree/set containing the vertex i : Use Path halving */ -template< typename Scalar,typename Index> -int SparseLUBase::etree_find (int i, IndexVector& pp) +template +int etree_find (int i, IndexVector& pp) { int p = pp(i); // Parent int gp = pp(p); // Grand parent @@ -50,45 +52,53 @@ int SparseLUBase::etree_find (int i, IndexVector& pp) * NOTE : The matrix is supposed to be in column-major format. * */ -template -int SparseLUBase::LU_sp_coletree(const MatrixType& mat, IndexVector& parent) +template +int coletree(const MatrixType& mat, IndexVector& parent, IndexVector& firstRowElt) { - int nc = mat.cols(); // Number of columns - int nr = mat.rows(); // Number of rows - + typedef typename MatrixType::Index Index; + Index nc = mat.cols(); // Number of columns + Index m = mat.rows(); IndexVector root(nc); // root of subtree of etree root.setZero(); IndexVector pp(nc); // disjoint sets pp.setZero(); // Initialize disjoint sets - IndexVector firstcol(nr); // First nonzero column in each row - + parent.resize(mat.cols()); //Compute first nonzero column in each row int row,col; - firstcol.setConstant(nc); //for (row = 0; row < nr; firstcol(row++) = nc); + firstRowElt.resize(m); + firstRowElt.setConstant(nc); + firstRowElt.segment(0, nc).setLinSpaced(nc, 0, nc-1); + bool found_diag; for (col = 0; col < nc; col++) { for (typename MatrixType::InnerIterator it(mat, col); it; ++it) - { // Is it necessary to browse the whole matrix, the lower part should do the job ?? + { row = it.row(); - firstcol(row) = (std::min)(firstcol(row), col); + firstRowElt(row) = (std::min)(firstRowElt(row), col); } } /* Compute etree by Liu's algorithm for symmetric matrices, - except use (firstcol[r],c) in place of an edge (r,c) of A. + except use (firstRowElt[r],c) in place of an edge (r,c) of A. Thus each row clique in A'*A is replaced by a star centered at its first vertex, which has the same fill. */ int rset, cset, rroot; for (col = 0; col < nc; col++) { + found_diag = false; pp(col) = col; cset = col; root(cset) = col; parent(col) = nc; - for (typename MatrixType::InnerIterator it(mat, col); it; ++it) + /* The diagonal element is treated here even if it does not exist in the matrix + * hence the loop is executed once more */ + for (typename MatrixType::InnerIterator it(mat, col); it||!found_diag; ++it) { // A sequence of interleaved find and union is performed - row = firstcol(it.row()); + int i = col; + if(it) i = it.index(); + if (i == col) found_diag = true; + row = firstRowElt(i); if (row >= col) continue; - rset = etree_find(row, pp); // Find the name of the set containing row + rset = internal::etree_find(row, pp); // Find the name of the set containing row rroot = root(rset); if (rroot != col) { @@ -106,8 +116,8 @@ int SparseLUBase::LU_sp_coletree(const MatrixType& mat, IndexVecto * Depth-first search from vertex n. No recursion. * This routine was contributed by Cédric Doucet, CEDRAT Group, Meylan, France. */ -template -void SparseLUBase::LU_nr_etdfs (int n, IndexVector& parent, IndexVector& first_kid, IndexVector& next_kid, IndexVector& post, int postnum) +template +void nr_etdfs (int n, IndexVector& parent, IndexVector& first_kid, IndexVector& next_kid, IndexVector& post, int postnum) { int current = n, first, next; while (postnum != n) @@ -152,8 +162,8 @@ void SparseLUBase::LU_nr_etdfs (int n, IndexVector& parent, IndexV * \param parent Input tree * \param post postordered tree */ -template -void SparseLUBase::LU_TreePostorder(int n, IndexVector& parent, IndexVector& post) +template +void treePostorder(int n, IndexVector& parent, IndexVector& post) { IndexVector first_kid, next_kid; // Linked list of children int postnum; @@ -174,7 +184,8 @@ void SparseLUBase::LU_TreePostorder(int n, IndexVector& parent, In // Depth-first search from dummy root vertex #n postnum = 0; - LU_nr_etdfs(n, parent, first_kid, next_kid, post, postnum); + internal::nr_etdfs(n, parent, first_kid, next_kid, post, postnum); } +} // internal #endif \ No newline at end of file diff --git a/Eigen/src/SparseLU/SparseLU_heap_relax_snode.h b/Eigen/src/SparseLU/SparseLU_heap_relax_snode.h index 69e1d4da9..7f6b3bab2 100644 --- a/Eigen/src/SparseLU/SparseLU_heap_relax_snode.h +++ b/Eigen/src/SparseLU/SparseLU_heap_relax_snode.h @@ -44,7 +44,7 @@ void SparseLUBase::LU_heap_relax_snode (const int n, IndexVector& // The etree may not be postordered, but its heap ordered IndexVector post; - LU_TreePostorder(n, et, post); // Post order etree + internal::treePostorder(n, et, post); // Post order etree IndexVector inv_post(n+1); int i; for (i = 0; i < n+1; ++i) inv_post(post(i)) = i; // inv_post = post.inverse()??? diff --git a/Eigen/src/SparseQR/CMakeLists.txt b/Eigen/src/SparseQR/CMakeLists.txt new file mode 100644 index 000000000..15009e6b5 --- /dev/null +++ b/Eigen/src/SparseQR/CMakeLists.txt @@ -0,0 +1,6 @@ +FILE(GLOB Eigen_SparseQRSupport_SRCS "*.h") + +INSTALL(FILES + ${Eigen_SparseQRSupport_SRCS} + DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/SparseQRSupport/ COMPONENT Devel + ) diff --git a/Eigen/src/SparseQR/SparseQR.h b/Eigen/src/SparseQR/SparseQR.h new file mode 100644 index 000000000..f1e5509dd --- /dev/null +++ b/Eigen/src/SparseQR/SparseQR.h @@ -0,0 +1,481 @@ +#ifndef EIGEN_SPARSE_QR_H +#define EIGEN_SPARSE_QR_H +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2012 Desire Nuentsa +// Copyright (C) 2012 Gael Guennebaud +// +// 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/. + + +namespace Eigen { +#include "../SparseLU/SparseLU_Coletree.h" + +template class SparseQR; +template struct SparseQRMatrixQReturnType; +template struct SparseQRMatrixQTransposeReturnType; +template struct SparseQR_QProduct; +namespace internal { + template struct traits > + { + typedef typename SparseQRType::MatrixType ReturnType; + }; + template struct traits > + { + typedef typename SparseQRType::MatrixType ReturnType; + }; + template struct traits > + { + typedef typename Derived::PlainObject ReturnType; + }; +} // End namespace internal + +/** + * \ingroup SparseQRSupport_Module + * \class SparseQR + * \brief Sparse left-looking QR factorization + * + * This class is used to perform a left-looking QR decomposition + * of sparse matrices. The result is then used to solve linear leasts_square systems. + * Clearly, a QR factorization is returned such that A*P = Q*R where : + * + * P is the column permutation. Use colsPermutation() to get it. + * + * Q is the orthogonal matrix represented as Householder reflectors. + * Use matrixQ() to get an expression and matrixQ().transpose() to get the transpose. + * You can then apply it to a vector. + * + * R is the sparse triangular factor. Use matrixQR() to get it as SparseMatrix. + * + * NOTE This is not a rank-revealing QR decomposition. + * + * \tparam _MatrixType The type of the sparse matrix A, must be a column-major SparseMatrix<> + * \tparam _OrderingType The fill-reducing ordering method. See \link OrderingMethods_Module + * OrderingMethods_Module \endlink for the list of built-in and external ordering methods. + * + * + */ +template +class SparseQR +{ + public: + typedef _MatrixType MatrixType; + typedef _OrderingType OrderingType; + typedef typename MatrixType::Scalar Scalar; + typedef typename MatrixType::RealScalar RealScalar; + typedef typename MatrixType::Index Index; + typedef SparseMatrix QRMatrixType; + typedef Matrix IndexVector; + typedef Matrix ScalarVector; + typedef PermutationMatrix PermutationType; + public: + SparseQR () : m_isInitialized(false),m_analysisIsok(false) + { } + + SparseQR(const MatrixType& mat) : m_isInitialized(false),m_analysisIsok(false) + { + compute(mat); + } + void compute(/*const*/ MatrixType& mat) + { + analyzePattern(mat); + factorize(mat); + } + void analyzePattern(const MatrixType& mat); + void factorize(/*const*/ MatrixType& mat); + + /** + * Get the number of rows of the triangular matrix. + */ + inline Index rows() const { return m_R.rows(); } + + /** + * Get the number of columns of the triangular matrix. + */ + inline Index cols() const { return m_R.cols();} + + /** + * This is the number of rows in the input matrix and the Q matrix + */ + inline Index rowsQ() const {return m_pmat.rows(); } + + /// Get the sparse triangular matrix R. It is a sparse matrix + MatrixType matrixQR() const + { + return m_R; + } + /// Get an expression of the matrix Q + SparseQRMatrixQReturnType matrixQ() const + { + return SparseQRMatrixQReturnType(*this); + } + + /// Get the permutation that was applied to columns of A + PermutationType colsPermutation() const + { + eigen_assert(m_isInitialized && "Decomposition is not initialized."); + return m_perm_c; + } + template + bool _solve(const MatrixBase &B, MatrixBase &dest) const + { + eigen_assert(m_isInitialized && "The factorization should be called first, use compute()"); + eigen_assert(this->rowsQ() == B.rows() && "SparseQR::solve() : invalid number of rows in the right hand side matrix"); + Index rank = this->matrixQR().cols(); + // Compute Q^T * b; + dest = this->matrixQ().transpose() * B; + // Solve with the triangular matrix R + Dest y; + y = this->matrixQR().template triangularView().solve(dest.derived().topRows(rank)); + // Apply the column permutation + if (m_perm_c.size()) + dest.topRows(rank) = colsPermutation().inverse() * y; + else + dest = y; + m_info = Success; + return true; + } + /** \returns the solution X of \f$ A X = B \f$ using the current decomposition of A. + * + * \sa compute() + */ + template + inline const internal::solve_retval solve(const MatrixBase& B) const + { + eigen_assert(m_isInitialized && "The factorization should be called first, use compute()"); + eigen_assert(this->rowsQ() == B.rows() && "SparseQR::solve() : invalid number of rows in the right hand side matrix"); + return internal::solve_retval(*this, B.derived()); + } + /** \brief Reports whether previous computation was successful. + * + * \returns \c Success if computation was succesful, + * \c NumericalIssue if the QR factorization reports a problem + * \c InvalidInput if the input matrix is invalid + * + * \sa iparm() + */ + ComputationInfo info() const + { + eigen_assert(m_isInitialized && "Decomposition is not initialized."); + return m_info; + } + protected: + bool m_isInitialized; + bool m_analysisIsok; + bool m_factorizationIsok; + mutable ComputationInfo m_info; + QRMatrixType m_pmat; // Temporary matrix + QRMatrixType m_R; // The triangular factor matrix + QRMatrixType m_Q; // The orthogonal reflectors + ScalarVector m_hcoeffs; // The Householder coefficients + PermutationType m_perm_c; // Column permutation + PermutationType m_perm_r; // Column permutation + IndexVector m_etree; // Column elimination tree + IndexVector m_firstRowElt; // First element in each row + IndexVector m_found_diag_elem; // Existence of diagonal elements + template friend struct SparseQR_QProduct; + +}; +/** + * \brief Preprocessing step of a QR factorization + * + * In this step, the fill-reducing permutation is computed and applied to the columns of A + * and the column elimination tree is computed as well. + * \note In this step it is assumed that there is no empty row in the matrix \p mat + */ +template +void SparseQR::analyzePattern(const MatrixType& mat) +{ + // Compute the column fill reducing ordering + OrderingType ord; + ord(mat, m_perm_c); + Index n = mat.cols(); + Index m = mat.rows(); + // Permute the input matrix... only the column pointers are permuted + m_pmat = mat; + m_pmat.uncompress(); + for (int i = 0; i < n; i++) + { + Index p = m_perm_c.size() ? m_perm_c.indices()(i) : i; +// m_pmat.outerIndexPtr()[i] = mat.outerIndexPtr()[p]; +// m_pmat.innerNonZeroPtr()[i] = mat.outerIndexPtr()[p+1] - mat.outerIndexPtr()[p]; + m_pmat.outerIndexPtr()[p] = mat.outerIndexPtr()[i]; + m_pmat.innerNonZeroPtr()[p] = mat.outerIndexPtr()[i+1] - mat.outerIndexPtr()[i]; + } + // Compute the column elimination tree of the permuted matrix + internal::coletree(m_pmat, m_etree, m_firstRowElt/*, m_found_diag_elem*/); + + m_R.resize(n, n); + m_Q.resize(m, m); + // Allocate space for nonzero elements : rough estimation + m_R.reserve(2*mat.nonZeros()); //FIXME Get a more accurate estimation through symbolic factorization with the etree + m_Q.reserve(2*mat.nonZeros()); + m_hcoeffs.resize(n); + m_analysisIsok = true; +} + +/** + * \brief Perform the numerical QR factorization of the input matrix + * + * \param mat The sparse column-major matrix + */ +template +void SparseQR::factorize(MatrixType& mat) +{ + eigen_assert(m_analysisIsok && "analyzePattern() should be called before this step"); + Index m = mat.rows(); + Index n = mat.cols(); + IndexVector mark(m); mark.setConstant(-1);// Record the visited nodes + IndexVector Ridx(n),Qidx(m); // Store temporarily the row indexes for the current column of R and Q + Index nzcolR, nzcolQ; // Number of nonzero for the current column of R and Q + Index pcol; + ScalarVector tval(m); tval.setZero(); // Temporary vector + IndexVector iperm(m); + bool found_diag; + if (m_perm_c.size()) + for(int i = 0; i < m; i++) iperm(m_perm_c.indices()(i)) = i; + else + iperm.setLinSpaced(m, 0, m-1); + + // Left looking QR factorization : Compute a column of R and Q at a time + for (Index col = 0; col < n; col++) + { + m_R.startVec(col); + m_Q.startVec(col); + mark(col) = col; + Qidx(0) = col; + nzcolR = 0; nzcolQ = 1; + pcol = iperm(col); + found_diag = false; + // Find the nonzero locations of the column k of R, + // i.e All the nodes (with indexes lower than k) reachable through the col etree rooted at node k + for (typename MatrixType::InnerIterator itp(mat, pcol); itp || !found_diag; ++itp) + { + Index curIdx = col; + if (itp) curIdx = itp.row(); + if(curIdx == col) found_diag = true; + // Get the nonzeros indexes of the current column of R + Index st = m_firstRowElt(curIdx); // The traversal of the etree starts here + if (st < 0 ) + { + std::cerr << " Empty row found during Numerical factorization ... Abort \n"; + m_info = NumericalIssue; + return; + } + // Traverse the etree + Index bi = nzcolR; + for (; mark(st) != col; st = m_etree(st)) + { + Ridx(nzcolR) = st; // Add this row to the list + mark(st) = col; // Mark this row as visited + nzcolR++; + } + // Reverse the list to get the topological ordering + Index nt = nzcolR-bi; + for(int i = 0; i < nt/2; i++) std::swap(Ridx(bi+i), Ridx(nzcolR-i-1)); + + // Copy the current row value of mat + if (itp) tval(curIdx) = itp.value(); + else tval(curIdx) = Scalar(0.); + + // Compute the pattern of Q(:,k) + if (curIdx > col && mark(curIdx) < col) + { + Qidx(nzcolQ) = curIdx; // Add this row to the pattern of Q + mark(curIdx) = col; // And mark it as visited + nzcolQ++; + } + } + + //Browse all the indexes of R(:,col) in reverse order + for (Index i = nzcolR-1; i >= 0; i--) + { + Index curIdx = Ridx(i); + // Apply the householder vector to tval + Scalar tdot(0.); + //First compute q'*tval + for (typename QRMatrixType::InnerIterator itq(m_Q, curIdx); itq; ++itq) + { + tdot += internal::conj(itq.value()) * tval(itq.row()); + } + tdot *= m_hcoeffs(curIdx); + // Then compute tval = tval - q*tau + for (typename QRMatrixType::InnerIterator itq(m_Q, curIdx); itq; ++itq) + { + tval(itq.row()) -= itq.value() * tdot; + } + //With the topological ordering, updates for curIdx are fully done at this point + m_R.insertBackByOuterInnerUnordered(col, curIdx) = tval(curIdx); + tval(curIdx) = Scalar(0.); + + //Detect fill-in for the current column of Q + if(m_etree(curIdx) == col) + { + for (typename QRMatrixType::InnerIterator itq(m_Q, curIdx); itq; ++itq) + { + Index iQ = itq.row(); + if (mark(iQ) < col) + { + Qidx(nzcolQ++) = iQ; // Add this row to the pattern of Q + mark(iQ) = col; //And mark it as visited + } + } + } + } // End update current column of R + + // Record the current (unscaled) column of V. + for (Index itq = 0; itq < nzcolQ; ++itq) + { + Index iQ = Qidx(itq); + m_Q.insertBackByOuterInnerUnordered(col,iQ) = tval(iQ); + tval(iQ) = Scalar(0.); + } + // Compute the new Householder reflection + RealScalar sqrNorm =0.; + Scalar tau; RealScalar beta; + typename QRMatrixType::InnerIterator itq(m_Q, col); + Scalar c0 = (itq) ? itq.value() : Scalar(0.); + //First, the squared norm of Q((col+1):m, col) + if(itq) ++itq; + for (; itq; ++itq) + { + sqrNorm += internal::abs2(itq.value()); + } + if(sqrNorm == RealScalar(0) && internal::imag(c0) == RealScalar(0)) + { + tau = RealScalar(0); + beta = internal::real(c0); + typename QRMatrixType::InnerIterator it(m_Q,col); + it.valueRef() = 1; //FIXME A row permutation should be performed at this point + } + else + { + beta = std::sqrt(internal::abs2(c0) + sqrNorm); + if(internal::real(c0) >= RealScalar(0)) + beta = -beta; + typename QRMatrixType::InnerIterator it(m_Q,col); + it.valueRef() = 1; + for (++it; it; ++it) + { + it.valueRef() /= (c0 - beta); + } + tau = internal::conj((beta-c0) / beta); + + } + m_hcoeffs(col) = tau; + m_R.insertBackByOuterInnerUnordered(col, col) = beta; + } + // Finalize the column pointers of the sparse matrices R and Q + m_R.finalize(); m_R.makeCompressed(); + m_Q.finalize(); m_Q.makeCompressed(); + m_isInitialized = true; + m_factorizationIsok = true; + m_info = Success; + +} + +namespace internal { + +template +struct solve_retval, Rhs> + : solve_retval_base, Rhs> +{ + typedef SparseQR<_MatrixType,OrderingType> Dec; + EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs) + + template void evalTo(Dest& dst) const + { + dec()._solve(rhs(),dst); + } +}; + +} // end namespace internal + +template +struct SparseQR_QProduct : ReturnByValue > +{ + typedef typename SparseQRType::QRMatrixType MatrixType; + typedef typename SparseQRType::Scalar Scalar; + typedef typename SparseQRType::Index Index; + // Get the references + SparseQR_QProduct(const SparseQRType& qr, const Derived& other, bool transpose) : + m_qr(qr),m_other(other),m_transpose(transpose) {} + inline Index rows() const { return m_transpose ? m_qr.rowsQ() : m_qr.cols(); } + inline Index cols() const { return m_other.cols(); } + + // Assign to a vector + template + void evalTo(DesType& res) const + { + Index m = m_qr.rowsQ(); + Index n = m_qr.cols(); + if (m_transpose) + { + eigen_assert(m_qr.m_Q.rows() == m_other.rows() && "Non conforming object sizes"); + // Compute res = Q' * other : + res = m_other; + for (Index k = 0; k < n; k++) + { + Scalar tau; + // Or alternatively + tau = m_qr.m_Q.col(k).tail(m-k).dot(res.tail(m-k)); + tau = tau * m_qr.m_hcoeffs(k); + res -= tau * m_qr.m_Q.col(k); + } + } + else + { + eigen_assert(m_qr.m_Q.cols() == m_other.rows() && "Non conforming object sizes"); + // Compute res = Q * other : + res = m_other; + for (Index k = n-1; k >=0; k--) + { + Scalar tau; + tau = m_qr.m_Q.col(k).tail(m-k).dot(res.tail(m-k)); + tau = tau * m_qr.m_hcoeffs(k); + res -= tau * m_qr.m_Q.col(k); + } + } + } + + const SparseQRType& m_qr; + const Derived& m_other; + bool m_transpose; +}; + +template +struct SparseQRMatrixQReturnType{ + + SparseQRMatrixQReturnType(const SparseQRType& qr) : m_qr(qr) {} + template + SparseQR_QProduct operator*(const MatrixBase& other) + { + return SparseQR_QProduct(m_qr,other.derived(),false); + } + SparseQRMatrixQTransposeReturnType adjoint() const + { + return SparseQRMatrixQTransposeReturnType(m_qr); + } + // To use for operations with the transpose of Q + SparseQRMatrixQTransposeReturnType transpose() const + { + return SparseQRMatrixQTransposeReturnType(m_qr); + } + const SparseQRType& m_qr; +}; + +template +struct SparseQRMatrixQTransposeReturnType{ + SparseQRMatrixQTransposeReturnType(const SparseQRType& qr) : m_qr(qr) {} + template + SparseQR_QProduct operator*(const MatrixBase& other) + { + return SparseQR_QProduct(m_qr,other.derived(), true); + } + const SparseQRType& m_qr; +}; +} // end namespace Eigen +#endif \ No newline at end of file diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 067f4a83a..f676ee9eb 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -219,6 +219,7 @@ ei_add_test(simplicial_cholesky) ei_add_test(conjugate_gradient) ei_add_test(bicgstab) ei_add_test(sparselu) +ei_add_test(sparseqr) # ei_add_test(denseLM) diff --git a/test/sparseqr.cpp b/test/sparseqr.cpp new file mode 100644 index 000000000..d34f7c48d --- /dev/null +++ b/test/sparseqr.cpp @@ -0,0 +1,62 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2012 Desire 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 +#include "sparse.h" +#include + + +template +int generate_sparse_rectangular_problem(MatrixType& A, DenseMat& dA, int maxRows = 300, int maxCols = 300) +{ + eigen_assert(maxRows >= maxCols); + typedef typename MatrixType::Scalar Scalar; + int rows = internal::random(1,maxRows); + int cols = internal::random(1,rows); + double density = (std::max)(8./(rows*cols), 0.01); + + A.resize(rows,rows); + dA.resize(rows,rows); + initSparse(density, dA, A,ForceNonZeroDiag); + A.makeCompressed(); + return rows; +} + +template void test_sparseqr_scalar() +{ + typedef SparseMatrix MatrixType; + MatrixType A; + Matrix dA; + typedef Matrix DenseVector; + DenseVector refX,x,b; + SparseQR > solver; + generate_sparse_rectangular_problem(A,dA); + + int n = A.cols(); + b = DenseVector::Random(n); + solver.compute(A); + if (solver.info() != Success) + { + std::cerr << "sparse QR factorization failed\n"; + exit(0); + return; + } + x = solver.solve(b); + if (solver.info() != Success) + { + std::cerr << "sparse QR factorization failed\n"; + exit(0); + return; + } + //Compare with a dense QR solver + refX = dA.colPivHouseholderQr().solve(b); + VERIFY(x.isApprox(refX,test_precision())); +} +void test_sparseqr() +{ + CALL_SUBTEST_1(test_sparseqr_scalar()); + CALL_SUBTEST_2(test_sparseqr_scalar >()); +} \ No newline at end of file