Re: [eigen] request for comments / help

[Gael Guennebaud <gael.guennebaud@xxxxxxxxx>]

Hi,

I guess the problem is on the type of V and thus CR which have to be
of dynamic size:

enum {
  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
};

typedef Matrix<Scalar,Dynamic,Dynamic,MaxColsAtCompileTime,MaxColsAtCompileTime>
MatrixVType;

MatrixVType V;
LLT<MatrixVType> CR;


Same for the matrix y...


gael.

On Mon, Jun 11, 2012 at 11:37 AM, Helmut Jarausch
<jarausch@xxxxxxxxxxxxxxxxxxx> wrote:
> #include <Eigen/QR>
>
>
> namespace Eigen {
> template<typename _MatrixType> class QR_Cholesky : public
> ColPivHouseholderQR<_MatrixType>
> {
> public:
>  typedef _MatrixType MatrixType;
>  typedef typename ColPivHouseholderQR<MatrixType>::Index Index;
>
>  QR_Cholesky() : ColPivHouseholderQR<MatrixType>(), my_isInitialized(false)
> {}
>  QR_Cholesky(const MatrixType& matrix) : ColPivHouseholderQR<MatrixType>() {
>    compute(matrix);
>  }
>  QR_Cholesky& compute(const MatrixType& matrix);
>  MatrixType solve(const MatrixType& rhs) const;
>
> protected:
>  bool my_isInitialized;
>
>  MatrixType V;
>  LLT<MatrixType> CR;
> };
>
> template<typename MatrixType>
> QR_Cholesky<MatrixType>& QR_Cholesky<MatrixType>::compute(const MatrixType&
> matrix) {
>  ColPivHouseholderQR<MatrixType>::compute(matrix);
>  Index q= this->dimensionOfKernel();
>  Index Rank= this->rank();
>  if ( q > 0 ) {
>    V = this->matrixQR().topRightCorner(Rank,q);
>    this->matrixQR().topLeftCorner(Rank,Rank).template
> triangularView<Upper>().solveInPlace(V);
>    CR.compute(V.transpose()*V + MatrixType::Identity(q,q));
>  }
>  my_isInitialized= true;
>  return *this;
> }
>
>
> template<typename MatrixType>
> MatrixType QR_Cholesky<MatrixType>::solve(const MatrixType& rhs) const {
>  eigen_assert(my_isInitialized && "QR_Cholesky is not initialized.");
>  eigen_assert(this->rows()==rhs.rows());
>  Index Rank= this->rank();
>  MatrixType Sol(this->cols(),rhs.cols());
>  MatrixType y = (this->householderQ().inverse() *
> rhs).topLeftCorner(Rank,rhs.cols());
>  this->matrixQR().topLeftCorner(Rank,Rank).template
> triangularView<Upper>().solveInPlace(y);
>  Index q = this->dimensionOfKernel();
>  if ( q > 0 ) {
>    MatrixType z = CR.solve(V.transpose()*y);
>    Sol.topLeftCorner(Rank,rhs.cols())= y-V*z;
>    Sol.bottomLeftCorner(q,rhs.cols())=z;
>  } else
>    Sol= y;
>  return this->colsPermutation()*Sol;
> }
>
> }  // End of namespace Eigen
>
>
> // T E S T
>
> #include <iostream>
> using std::cout;  using std::endl;
> using Eigen::MatrixXd;  using Eigen::VectorXd;
>
>
>
> int main() {
>  MatrixXd A(3,4);
>  VectorXd b(3);
>  A << 1,2,3,4,
>       5,6,7,8,
>       9,10,11,12;
>
>  b << 30,70,110;
>
>  double eps=1E-8;
>  Eigen::QR_Cholesky<MatrixXd> QRChol;
> //   Eigen::QR_Cholesky<MatrixXd>  QRFac(A);
>  QRChol.setThreshold(eps);
>  QRChol.compute(A);
>  VectorXd x = QRChol.solve(b);
>  cout << "LSQ Solution: \n" << x << endl;
> }