STK::Qr Class Reference
[Subproject STKAlgebra::Algebra]

The class Qr perform the QR decomposition of a Matrix. More...

#include <STK_Qr.h>

List of all members.

Public Member Functions
	Qr (const Matrix &A=Matrix(), const bool &ref=false)
	~Qr ()
void	clearQ ()
void	clear ()
void	compQr ()
void	compQ ()
void	NewQr (const Matrix &A)
Qr &	operator= (const Qr &S)
bool	isCompQ () const
	< Is Q computed ?
const Matrix &	getQ () const
const MatrixUpperTriangular &	getR () const
Qr &	popBackCols (const Integer &n=1)
Qr &	eraseCol (const Integer &pos)
Qr &	pushBackCol (const Vector &T)
Qr &	insertCol (const Vector &T, const Integer &pos)
Protected Attributes
Matrix	Q_
	Array2D Q.
MatrixUpperTriangular	R_
	Array2D R.
Integer	ncolr_
	Number of cols of R actually computed.
Integer	ncolq_
	Number of cols used by Q and Number of rows of R_.
Integer	nrowq_
	Number of rows used by Q.
bool	compq_
	is Q computed ?
Private Member Functions
void	qr ()
	Private function for computing the qr.

---

Detailed Description

The user have to take care to the behavior of the class if withq_ is false: in this case the QR decompostion cannot be updated if he want to add a col.

Also the Addcol or insertCols methods will fail if ncolmax is not large enough.

• Input: A matrix (nrow,ncolq)
• Output:
  1. Q matrix (nrow,ncol) orthonormal.
  2. R matrix square (ncol,ncolq) upper triangular
• A=QR.

Definition at line 71 of file STK_Qr.h.

---

Constructor & Destructor Documentation

STK::Qr::Qr	(	const Matrix &	A = Matrix(),
		const bool &	ref = false
	)

Default constructor.

Parameters:

	A	the matrix to decompose
	ref	true if we overwrite A

Definition at line 72 of file STK_Qr.cpp.

References compQr().

      : Q_(A, ref)      // Creating Q
      , R_()            // Creating R
{ compQr();}

STK::Qr::~Qr

(

)

Dtor

Definition at line 227 of file STK_Qr.cpp.

{ ;}

---

Member Function Documentation

void STK::Qr::qr

(

)

[private]

Definition at line 115 of file STK_Qr.cpp.

References beta(), STK::dot(), STK::ITContainer1D< TYPE, TContainer1D >::front(), STK::house(), STK::min(), ncolq_, ncolr_, nrowq_, Q_, R_, and STK::ITContainer2D< TYPE, TContainerHo, TContainerVe, TContainer2D >::resize().

Referenced by compQr().

{
  Integer niter = min(nrowq_,ncolr_);   // number of iterations
  R_.resize(nrowq_, ncolr_);
  R_ = 0.0;                   // initialize to 0.0
  
  /*------------------------------------------------------------------*/
  /* Main loop.                                                       */
  for (Integer iter=1, iter1=2; iter<=niter; iter++, iter1++)
  { 
    // A ref on the row of the matrix R_
    Point  Rrow0(R_, Inx(iter, ncolr_), iter);
    // A ref on the row of the matrix R_
    Point  Qrow0(Q_, Inx(iter, ncolq_), iter);
    // A ref on the col iter of the matrix Q_ : will contain the
    // Householder vector
    Vector u(Q_, Inx(iter, nrowq_), iter);
    // compute the Householder vector of the current col
    Rrow0[iter] = house(u);
    // get beta
    Real beta = u.front();
    if (beta)
    {
      // ref on the essential part of the Householder vector
      Vector v(u, Inx(iter1, nrowq_));
      // Apply Householder to next cols
      for (Integer j=iter1; j<=ncolr_; j++)
      {
        // Auxiliary data
        Real aux;
        // a ref on the jth column of Q_
        Vector Z(Q_, Inx(iter1, nrowq_), j);
        // save the  current row of R_
        Rrow0[j] = Qrow0[j] + (aux = ( Qrow0[j] + dot(v, Z)) * beta);
        // update the next cols of Q_ 
        Z += v * aux;
      }
    }
    else
    {
      // just copy the row iter in R_
      for (Integer j=iter1; j<=ncolr_; j++)
      {
        Rrow0[j] = Qrow0[j];
      }
    }
  }
}

void STK::Qr::clearQ

(

)

clear Q_.

Definition at line 231 of file STK_Qr.cpp.

References STK::IArray2D< TYPE, TArray2D >::clear(), compq_, nrowq_, and Q_.

{ Q_.clear(); nrowq_ = 0; compq_ = false;}

void STK::Qr::clear

(

)

clear Q_ and R_.

Definition at line 236 of file STK_Qr.cpp.

References STK::IArray2D< TYPE, TArray2D >::clear(), Q_, and R_.

Referenced by NewQr().

{
  Q_.clear();
  R_.clear();
}

void STK::Qr::compQr

(

)

Compute the QR decomposition.

Definition at line 90 of file STK_Qr.cpp.

References compq_, STK::ITContainer2D< TYPE, TContainerHo, TContainerVe, TContainer2D >::empty(), ncolq_, ncolr_, nrowq_, Q_, qr(), STK::IArray2D< TYPE, TArray2D >::shift(), STK::ITContainer2D< TYPE, TContainerHo, TContainerVe, TContainer2D >::sizeHo(), and STK::ITContainer2D< TYPE, TContainerHo, TContainerVe, TContainer2D >::sizeVe().

Referenced by NewQr(), and Qr().

{
  if (Q_.empty())     // if the container is empty
  {
    ncolr_ =0;
    ncolq_ =0;
    nrowq_ =0;
    compq_  = true;  // Q_ is computed

    return;
  }

  Q_.shift(1,1);      // translate the beg to 1
  ncolr_ = Q_.sizeHo(); // Number of cols of R
  ncolq_ = Q_.sizeHo(); // Number of col of Q
  nrowq_ = Q_.sizeVe(); // Number of rows of Q

  compq_  = false;    // Q_ is not computed

  // compute QR decomposition
  qr();
}

void STK::Qr::compQ

(

)

Compute Q (to use after compQr). Without effect compq_ == true

Definition at line 166 of file STK_Qr.cpp.

References beta(), compq_, STK::dot(), STK::min(), ncolq_, nrowq_, STK::IArray2D< TYPE, TArray2D >::popBackCols(), STK::IArray2D< TYPE, TArray2D >::pushBackCols(), and Q_.

Referenced by eraseCol(), insertCol(), and pushBackCol().

{
  // if Q_ is computed yet
  if (compq_) return;
  // number of non zero cols of Q_  
  Integer ncol  = min(nrowq_, ncolq_);
  // Q_ is square
  if (nrowq_ < ncolq_)
  {
     Q_.popBackCols(ncolq_-nrowq_);
  }
  else
  { 
    Q_.pushBackCols(nrowq_-ncolq_);
    // initialization of the remaining cols of Q_ to 0.0
    Q_[Inx(ncol+1,nrowq_)] = 0.0;
  }
  // the number of col_ is equal to the number of row
  ncolq_ = nrowq_;
  // Computation of Q_.
  // compute added col
  for (Integer iter=ncolq_; iter> ncol; --iter)
  { Q_(iter, iter) = 1.0;}
  // compute other cols
  for (Integer iter=ncol, iter1=ncol+1; iter>=1; iter--, iter1--)
  {
    // Get beta and test
    Real beta = Q_(iter,iter);
    if (beta)
    {
      // ref of the row iter
      Point P(Q_, Inx(iter,ncolq_), iter);
      // ref of the col iter
      Vector X(Q_, Inx(iter1,nrowq_), iter);
      // Update the cols from iter+1 to ncol
      for (Integer j=iter1; j<=ncolq_; j++)
      { Real aux;
        Vector Y(Q_, Inx(iter1,nrowq_), j); // ref on the col j
        // Q_(iter, j) = beta * X'Y
        P[j] = (aux = dot( X, Y) * beta);
        // Q^j += aux * Q^iter
        Y += X * aux;
      }
      P[iter] = 1.0 + beta;
      // Q^iter *= beta
      X *= beta;
    }
    else // Q^iter = identity
    {
      Q_(iter, iter) = 1.0;
      Q_(Inx(iter1,nrowq_), iter) = 0.0;
    }
    // update the col iter
    Q_(Inx(1,iter-1), iter) = 0.0;
  }
  // Q_ is now computed
  compq_ = true;
}

void STK::Qr::NewQr

(

const Matrix &

A

)

Compute the QR decomposition of the Matrix A.

Parameters:

A

the matrix to decompose

Definition at line 79 of file STK_Qr.cpp.

References clear(), compQr(), and Q_.

{
  clear();  // clearing old QR
  Q_ = A;   // Copy A in Q

  // compute QR decompositionn
  compQr();
}

Qr & STK::Qr::operator=

(

const Qr &

S

)

Operator = : overwrite the Qr with S.

Definition at line 244 of file STK_Qr.cpp.

References compq_, ncolq_, ncolr_, nrowq_, Q_, and R_.

{
  ncolr_ = S.ncolr_;       //< Number of cols of R actually computed
  ncolq_ = S.ncolq_;       //< Number of cols used by Q
  nrowq_ = S.nrowq_;       //< Number of rows used by Q
    
  compq_ = S.compq_;       //< Is Q computed ?

  Q_ = S.Q_;               //< Matrix V
  R_ = S.R_;               //< Singular values

  return *this;
}

bool STK::Qr::isCompQ

(

)

const [inline]

Definition at line 126 of file STK_Qr.h.

References compq_.

    { return compq_;}

const Matrix& STK::Qr::getQ

(

)

const [inline]

give Q.

Definition at line 130 of file STK_Qr.h.

References Q_.

    { return Q_;}

const MatrixUpperTriangular& STK::Qr::getR

(

)

const [inline]

give R.

Definition at line 134 of file STK_Qr.h.

References R_.

    { return R_;}

Qr & STK::Qr::popBackCols

(

const Integer &

n = 1

)

Delete the n last cols and update the QR decomposition.

Definition at line 262 of file STK_Qr.cpp.

References ncolr_, STK::IArray2D< TYPE, TArray2D >::popBackCols(), and R_.

{
  // delete n cols
  R_.popBackCols(n);
  ncolr_ -= n;
  return *this;
}

Qr & STK::Qr::eraseCol

(

const Integer &

pos

)

Delete the col pos and update the QR decomposition.

Definition at line 270 of file STK_Qr.cpp.

References STK::compGivens(), compQ(), compq_, STK::IArray2D< TYPE, TArray2D >::eraseCols(), STK::ITContainer2D< TYPE, TContainerHo, TContainerVe, TContainer2D >::firstCol(), STK::ITContainer2D< TYPE, TContainerHo, TContainerVe, TContainer2D >::getSizeHo(), STK::ITContainer2D< TYPE, TContainerHo, TContainerVe, TContainer2D >::getSizeVe(), STK::ITContainer2D< TYPE, TContainerHo, TContainerVe, TContainer2D >::lastCol(), STK::ITContainer2D< TYPE, TContainerHo, TContainerVe, TContainer2D >::lastRow(), STK::leftGivens(), STK::min(), ncolr_, Q_, R_, STK::rightGivens(), and STK::IArray2D< TYPE, TArray2D >::update().

{
#ifdef STK_BOUNDS_CHECK
  if (pos < R_.firstCol())
  { throw std::out_of_range("Qr::eraseCol(pos) "
                       "pos < R_.firstCol()");
  }
  if (R_.lastCol() < pos)
  { throw std::out_of_range("Qr::eraseCol(pos) "
                       "R_.lastCol() < pos");
  }
#endif
  // if Q_ is not computed yet
  if (!compq_) compQ();
  // compute the number of iteration for updating to zeroed
  Integer niter = R_.firstCol()-1+min(R_.getSizeVe(), R_.getSizeHo());
  // Zeroed the unwanted elements (z)
  for (Integer iter = pos+1; iter<=niter; iter++)
  {
    Real sinus, cosinus;
    // compute the Givens rotation
    R_(iter-1, iter) = compGivens( R_(iter-1, iter)
                                 , R_(iter, iter)
                                 , cosinus
                                 , sinus
                                 );
    R_(iter, iter)   = 0.0;
    // if necessary update R_ and Q_
    if (sinus)
    {
      // Update the next rows (iter1:ncolr_) of R_
      leftGivens(R_[Inx(iter+1, R_.lastCol())]
                , iter-1
                , iter
                , cosinus
                , sinus
                );
      // Update the cols of Q_
      rightGivens(Q_, iter-1, iter, cosinus, sinus);
    }
  }
  // erase the col pos
  R_.eraseCols(pos);
  ncolr_--;
  
  // update the range of the remaining cols of the container
  R_.update(Inx(pos, min(R_.lastRow(), R_.lastCol())));

  return *this;
}

Qr & STK::Qr::pushBackCol

(

const Vector &

T

)

Add a col with value T and update th QR decomposition.

Definition at line 323 of file STK_Qr.cpp.

References STK::compGivens(), compQ(), compq_, STK::MatrixUpperTriangular::compRangeVe(), STK::dot(), STK::ITContainer2D< TYPE, TContainerHo, TContainerVe, TContainer2D >::firstCol(), STK::ITContainer1D< TYPE, TContainer1D >::getRange(), STK::ITContainer2D< TYPE, TContainerHo, TContainerVe, TContainer2D >::getRangeVe(), STK::ITContainer2D< TYPE, TContainerHo, TContainerVe, TContainer2D >::lastCol(), ncolq_, ncolr_, STK::IArray2D< TYPE, TArray2D >::pushBackCols(), Q_, R_, and STK::rightGivens().

{
  // check conditions
  if (T.getRange() != Q_.getRangeVe())
  { throw std::runtime_error("Qr::pushBackCol(T) "
                        "T.range() != Q_.getRangeVe()");
  }
  // if Q_ is not computed yet
  if (!compq_) compQ();
  // Adding a col to R
  R_.pushBackCols();
  ncolr_++;
  // Create an auxiliary container
  Vector Rncolr(Q_.getRangeVe());
  // Multipliate T by Q'and put the result in Rncolr
  for (Integer j=Q_.firstCol(); j<=Q_.lastCol(); j++)
    Rncolr[j] = dot(Q_[j], T);
  // update Q_
  for (Integer iter = ncolq_-1, iter1 = ncolq_; iter>=ncolr_; iter--, iter1--)
  { 
    Real sinus, cosinus, y = Rncolr[iter], z = Rncolr[iter1] ;
    // compute the Givens rotition
    Rncolr[iter]  = compGivens(y, z, cosinus, sinus);
    // apply Givens rotation if necessary
    if (sinus)
    {
      // Update the cols of Q_
      rightGivens(Q_, iter, iter1, cosinus, sinus);
    }
  }
  // update R_
  R_[ncolr_] = Rncolr[R_.compRangeVe(ncolr_)];
  //
  return *this;
}

Qr & STK::Qr::insertCol	(	const Vector &	T,
		const Integer &	pos
	)

Add a col with value T at the position pos and update the QR decomposition.

Adding the jth col and update the QR decomposition.

Definition at line 361 of file STK_Qr.cpp.

References STK::compGivens(), compQ(), compq_, STK::MatrixUpperTriangular::compRangeVe(), STK::dot(), STK::ITContainer2D< TYPE, TContainerHo, TContainerVe, TContainer2D >::firstCol(), STK::ITContainer1D< TYPE, TContainer1D >::getRange(), STK::ITContainer2D< TYPE, TContainerHo, TContainerVe, TContainer2D >::getRangeVe(), STK::IArray2D< TYPE, TArray2D >::insertCols(), STK::ITContainer2D< TYPE, TContainerHo, TContainerVe, TContainer2D >::lastCol(), STK::ITContainer2D< TYPE, TContainerHo, TContainerVe, TContainer2D >::lastRow(), STK::leftGivens(), STK::min(), ncolq_, ncolr_, Q_, R_, STK::rightGivens(), and STK::IArray2D< TYPE, TArray2D >::update().

{
#ifdef STK_BOUNDS_CHECK
  if (pos < 1)
  { throw std::out_of_range("Qr::insertCol(T, pos) "
                       "j < 1");
  }
  if (ncolr_ < pos)
  { throw std::out_of_range("Qr::insertCol(T, pos) "
                       "ncolr_ < j");
  }
#endif
  if (T.getRange() != Q_.getRangeVe())
  { throw std::runtime_error("Qr::insertCol(T) "
                        "T.range() != Q_.getRangeVe()");
  }
  // if Q_ is not computed yet
  if (!compq_) compQ();
  // Adding a col to R
  R_.insertCols(pos);
  ncolr_++;
  // update the range of the remaining cols of R_
  R_.update(Inx(pos+1, min(R_.lastRow(), R_.lastCol())));
  for (Integer i=pos+1; i<= min(R_.lastRow(), R_.lastCol()); ++i)
    R_(i,i) = 0.0;
  // A ref on the last col of R_
  Vector Rpos(Q_.getRangeVe());
  // Multipliate T by Q'
  // we cannot use mult as we are using ncolq_ cols.
  for (Integer j=Q_.firstCol(); j<=Q_.lastCol(); j++)
    Rpos[j] = dot(Q_[j],T);
  // Zeroed the unwanted elements
  for (Integer iter= ncolq_-1, iter1= ncolq_; iter>=pos; iter--, iter1--)
  { 
    Real sinus, cosinus, y = Rpos[iter], z = Rpos[iter1] ;
    // compute the Givens rotation
    Rpos[iter]  = compGivens(y, z, cosinus, sinus);
    // apply Givens rotation if necessary
    if (sinus)
    {
      // Update the next rows (iter1:ncolr_) of R_
      leftGivens( R_[Inx(iter1, R_.lastCol())]
                       , iter
                       , iter1
                       , cosinus
                       , sinus
                       );
      // Update the cols of Q_
      rightGivens(Q_, iter, iter1, cosinus, sinus);
    }
  }
  // update R_
  R_[pos] = Rpos[R_.compRangeVe(pos)];
  R_.update(pos);

  //
  return *this;
}

---

Member Data Documentation

Matrix STK::Qr::Q_ [protected]

Definition at line 77 of file STK_Qr.h.

Referenced by clear(), clearQ(), compQ(), compQr(), eraseCol(), getQ(), insertCol(), NewQr(), operator=(), pushBackCol(), and qr().

MatrixUpperTriangular STK::Qr::R_ [protected]

Definition at line 78 of file STK_Qr.h.

Referenced by clear(), eraseCol(), getR(), insertCol(), operator=(), popBackCols(), pushBackCol(), and qr().

Integer STK::Qr::ncolr_ [protected]

Definition at line 80 of file STK_Qr.h.

Referenced by compQr(), eraseCol(), insertCol(), operator=(), popBackCols(), pushBackCol(), and qr().

Integer STK::Qr::ncolq_ [protected]

Definition at line 81 of file STK_Qr.h.

Referenced by compQ(), compQr(), insertCol(), operator=(), pushBackCol(), and qr().

Integer STK::Qr::nrowq_ [protected]

Definition at line 82 of file STK_Qr.h.

Referenced by clearQ(), compQ(), compQr(), operator=(), and qr().

bool STK::Qr::compq_ [protected]

Definition at line 84 of file STK_Qr.h.

Referenced by clearQ(), compQ(), compQr(), eraseCol(), insertCol(), isCompQ(), operator=(), and pushBackCol().

---

The documentation for this class was generated from the following files:

• STK_Qr.h
• STK_Qr.cpp