STK++ 1.0
Public Member Functions | Protected Attributes | Private Member Functions

STK::Qr Class Reference
Algebra.

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

#include <STK_Qr.h>

List of all members.

Public Member Functions

 Qr (Matrix const &A=Matrix(), bool ref=false)
 Default constructor.
virtual ~Qr ()
 virtual destructor
Qr & operator= (const Qr &S)
 Operator = : overwrite the Qr with S.
bool isCompQ () const
 Is Q computed ?
Matrix const & Q () const
 give the matrix Q of the QR decomposition.
const MatrixUpperTriangular & R () const
 give the matrix R of the QR decomposition.
void clear ()
 clear Q_ and R_.
void run ()
 Compute the QR decomposition.
void compQ ()
 Compute Q (to use after run).
void popBackCols (Integer const &n=1)
 Delete the n last columns and update the QR decomposition.
void eraseCol (Integer const &pos)
 Delete the column pos and update the QR decomposition.
void pushBackCol (Vector const &T)
 Add a column with value T and update th QR decomposition.
void insertCol (Vector const &T, Integer const &pos)
 Add a column with value T at the position pos and update the QR decomposition.

Protected Attributes

Matrix Q_
 Q Matrix of the QR decomposition.
MatrixUpperTriangular R_
 R Matrix of th QR decomposition.
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 ()
 Compute the qr decoposition of the matrix Q_.

Detailed Description

The class Qr perform the QR decomposition of a Matrix.

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 column.

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 62 of file STK_Qr.h.

Constructor & Destructor Documentation

home iovleff Developpement workspace stkpp projects Algebra src STK_Qr cpp STK::Qr::Qr ( Matrix const &  A = Matrix(),
bool  ref = false 
)

Default constructor.

Parameters:
Athe matrix to decompose
reftrue if we overwrite A

Definition at line 50 of file STK_Qr.cpp.

      : Q_(A, ref)      // Creating Q
      , R_()            // Creating R
{ run();}
STK::Qr::~Qr ( ) [virtual]

virtual destructor

Definition at line 193 of file STK_Qr.cpp.

Member Function Documentation

Qr & STK::Qr::operator= ( const Qr &  S)

Operator = : overwrite the Qr with S.

Definition at line 204 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]

Is Q computed ?

Returns:
true if Q_ is computed, false otherwise

Definition at line 94 of file STK_Qr.h.

References compq_.

    { return compq_;}
Matrix const& STK::Qr::Q ( ) const [inline]

give the matrix Q of the QR decomposition.

Returns:
the matrix Q of the QR decomposition

Definition at line 99 of file STK_Qr.h.

References Q_.

{ return Q_;}
const MatrixUpperTriangular& STK::Qr::R ( ) const [inline]

give the matrix R of the QR decomposition.

Returns:
the matrix R of the QR decomposition

Definition at line 103 of file STK_Qr.h.

References R_.

{ return R_;}
void STK::Qr::clear ( )

clear Q_ and R_.

Definition at line 196 of file STK_Qr.cpp.

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

{
  Q_.clear();
  R_.clear();
}
void STK::Qr::run ( )

Compute the QR decomposition.

Definition at line 56 of file STK_Qr.cpp.

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

{
  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 column 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 run).

After the run process, Q_ store the householder vector in its column. Call compQ, if you want to obtain Q in its true form. Without effect if (compq_ == true)

Definition at line 132 of file STK_Qr.cpp.

References beta(), compq_, STK::dot(), STK::min(), ncolq_, nrowq_, STK::Funct::P, 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_[Range(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_, Range(iter,ncolq_), iter);
      // ref of the column iter
      Vector X(Q_, Range(iter1,nrowq_), iter);
      // Update the cols from iter+1 to ncol
      for (Integer j=iter1; j<=ncolq_; j++)
      { Real aux;
        Vector Y(Q_, Range(iter1,nrowq_), j); // ref on the column 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_(Range(iter1,nrowq_), iter) = 0.0;
    }
    // update the column iter
    Q_(Range(1,iter-1), iter) = 0.0;
  }
  // Q_ is now computed
  compq_ = true;
}
void STK::Qr::popBackCols ( Integer const &  n = 1)

Delete the n last columns and update the QR decomposition.

Parameters:
nnumber of column to delete

Definition at line 221 of file STK_Qr.cpp.

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

{
  // delete n cols
  R_.popBackCols(n);
  ncolr_ -= n;
}
void STK::Qr::eraseCol ( Integer const &  pos)

Delete the column pos and update the QR decomposition.

Parameters:
posthe position of the column to delete

Definition at line 228 of file STK_Qr.cpp.

References STK::compGivens(), compQ(), compq_, STK::IArray2D< TYPE, TArray2D >::eraseCols(), STK::IContainer2D::firstCol(), STK::IContainer2D::lastCol(), STK::IContainer2D::lastRow(), STK::leftGivens(), STK::min(), ncolr_, Q_, R_, STK::rightGivens(), STK::IContainer2D::sizeHo(), STK::IContainer2D::sizeVe(), 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_.sizeVe(), R_.sizeHo());
  // 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_[Range(iter+1, R_.lastCol())]
                , iter-1
                , iter
                , cosinus
                , sinus
                );
      // Update the cols of Q_
      rightGivens(Q_, iter-1, iter, cosinus, sinus);
    }
  }
  // erase the column pos
  R_.eraseCols(pos);
  ncolr_--;
  
  // update the range of the remaining cols of the container
  R_.update(Range(pos, min(R_.lastRow(), R_.lastCol())));
}
void STK::Qr::pushBackCol ( Vector const &  T)

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

Parameters:
Tthe column to add

Definition at line 279 of file STK_Qr.cpp.

References STK::compGivens(), compQ(), compq_, STK::MatrixUpperTriangular::compRangeVe(), STK::dot(), STK::IContainer2D::firstCol(), STK::IContainer2D::lastCol(), ncolq_, ncolr_, STK::IArray2D< TYPE, TArray2D >::pushBackCols(), Q_, R_, STK::ITContainer1D< TYPE, TContainer1D >::range(), STK::IContainer2D::rangeVe(), and STK::rightGivens().

{
  // check conditions
  if (T.range() != Q_.rangeVe())
  { throw std::runtime_error("Qr::pushBackCol(T) "
                        "T.range() != Q_.rangeVe()");
  }
  // if Q_ is not computed yet
  if (!compq_) compQ();
  // Adding a column to R
  R_.pushBackCols();
  ncolr_++;
  // Create an auxiliary container
  Vector Rncolr(Q_.rangeVe());
  // 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_)];
}
void STK::Qr::insertCol ( Vector const &  T,
Integer const &  pos 
)

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

Parameters:
Tthe column to insert
posthe position of the column to insert

Definition at line 315 of file STK_Qr.cpp.

References STK::compGivens(), compQ(), compq_, STK::MatrixUpperTriangular::compRangeVe(), STK::dot(), STK::IContainer2D::firstCol(), STK::IArray2D< TYPE, TArray2D >::insertCols(), STK::IContainer2D::lastCol(), STK::IContainer2D::lastRow(), STK::leftGivens(), STK::min(), ncolq_, ncolr_, Q_, R_, STK::ITContainer1D< TYPE, TContainer1D >::range(), STK::IContainer2D::rangeVe(), 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_ < pos");
  }
#endif
#ifdef STK_DEBUG
  if (T.range() != Q_.rangeVe())
  { throw std::runtime_error("Qr::insertCol(T, pos) "
                             "T.range() != Q_.rangeVe()");
  }
#endif
  // if Q_ is not computed yet
  if (!compq_) compQ();
  // Adding a column to R
  R_.insertCols(pos);
  ncolr_++;
  // update the range of the remaining cols of R_
  R_.update(Range(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 column of R_
  Vector Rpos(Q_.rangeVe());
  // Multipliate T by Q'
  // we cannot use mult as we are using ncolq_ columns.
  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_[Range(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);
}
void STK::Qr::qr ( ) [private]

Compute the qr decoposition of the matrix Q_.

Definition at line 81 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::IContainer2D::resize().

Referenced by run().

{
  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_, Range(iter, ncolr_), iter);
    // A ref on the row of the matrix R_
    Point  Qrow0(Q_, Range(iter, ncolq_), iter);
    // A ref on the column iter of the matrix Q_ : will contain the
    // Householder vector
    Vector u(Q_, Range(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, Range(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_, Range(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];
      }
    }
  }
}

Member Data Documentation

Matrix STK::Qr::Q_ [protected]

Q Matrix of the QR decomposition.

Definition at line 66 of file STK_Qr.h.

Referenced by clear(), compQ(), eraseCol(), insertCol(), operator=(), pushBackCol(), Q(), qr(), and run().

MatrixUpperTriangular STK::Qr::R_ [protected]

R Matrix of th QR decomposition.

Definition at line 68 of file STK_Qr.h.

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

Integer STK::Qr::ncolr_ [protected]

Number of cols of R actually computed.

Definition at line 70 of file STK_Qr.h.

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

Integer STK::Qr::ncolq_ [protected]

Number of cols used by Q and Number of rows of R_.

Definition at line 72 of file STK_Qr.h.

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

Integer STK::Qr::nrowq_ [protected]

Number of rows used by Q.

Definition at line 74 of file STK_Qr.h.

Referenced by compQ(), operator=(), qr(), and run().

bool STK::Qr::compq_ [protected]

is Q computed ?

Definition at line 76 of file STK_Qr.h.

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

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