sparseqr/givensRotations.hh

Go to the documentation of this file.

// Givens rotations

#ifndef givensRotations_hh
#define givensRotations_hh

#include "sparseqr/sparseqr.hh"
#include "operator/matrix.hh"

namespace sparseqr {

  // ******************************************************* GivensRotations **

  template<typename F>
  class GivensRotations : public concepts::Operator<F> {
  public:
    template<class G>
    GivensRotations(const concepts::Space<G>& space,
                    const Qmatrix& q, bool transpose)
      : concepts::Operator<F>(space.dim(), space.dim())
      , q_(q), transpose_(transpose) {}
    GivensRotations(uint dim, const Qmatrix& q, bool transpose)
      : concepts::Operator<F>(dim, dim), q_(q), transpose_(transpose) {}
    virtual void operator()(const concepts::Function<F>& fncY,
                            concepts::Function<F>& fncX);
    void multiply(const concepts::Matrix<F>& fact,
                  concepts::Matrix<F>& dest) const;
    void multiplyFirst(const concepts::Matrix<F>& fact,
                       concepts::Matrix<F>& dest) const;
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  private:
    const Qmatrix& q_;
    const bool transpose_;
  };

  template<typename F>
  void GivensRotations<F>::operator()
    (const concepts::Function<F>& fncY, concepts::Function<F>& fncX) {
    conceptsAssert(fncY.dim() == this->dimY(), concepts::Assertion());
    conceptsAssert(fncX.dim() == this->dimX(), concepts::Assertion());
    fncX = fncY;
    const J* j = q_.qend;
    if (transpose_) j = q_.qbegin;
    while (j != 0) {
      conceptsAssert((uint)j->row1 < this->dimY(), concepts::Assertion());
      conceptsAssert((uint)j->row2 < this->dimY(), concepts::Assertion());
      F t = fncX(j->row1);
      if (transpose_) {
        fncX(j->row1) = j->c*fncX(j->row1) + j->s*fncX(j->row2);
        fncX(j->row2) = j->c*fncX(j->row2) - j->s*t;
      } else {
        fncX(j->row1) = j->c*fncX(j->row1) - j->s*fncX(j->row2);
        fncX(j->row2) = j->c*fncX(j->row2) + j->s*t;
      }
      if (transpose_) j = j->next;
      else  j = j->previous;
    }
  }

  /* The following routine does the same for a matrix as operator() above
     does for a vector.

     It has to compute \c dest = \c fact times Given's rotations whereas
     operator() computes \c fncX = Given's rotations times \c fncY.
     \c dest^T = Given's^T times \c fact^T => each column of \c fact^T is
     treated the same way as \c fncY. A column of \c fact^T is a row
     of \c fact.

     The precondition of the routine is \c fact = \c dest, therefore, we
     can work on \c dest only.
  */
  template<typename F>
  void GivensRotations<F>::multiply
  (const concepts::Matrix<F>& fact, concepts::Matrix<F>& dest) const {
    // check sizes of the matrices
    conceptsAssert(fact.dimX() == dest.dimX(), concepts::Assertion());
    conceptsAssert(fact.dimY() == dest.dimY(), concepts::Assertion());
    conceptsAssert(fact.dimY() == this->dimX(), concepts::Assertion());
    // loop over rows of dest
    for (uint i = 0; i < dest.dimX(); ++i) {
      const J* j = q_.qbegin;
      if (transpose_) j = q_.qend;
      while (j != 0) { // loop over Given's rotations
        conceptsAssert((uint)j->row1 < this->dimX(), concepts::Assertion());
        conceptsAssert((uint)j->row2 < this->dimX(), concepts::Assertion());
        F t = dest(i, j->row1);
        if (transpose_) {
          dest(i, j->row1) = j->c*dest(i, j->row1) - j->s*dest(i, j->row2);
          dest(i, j->row2) = j->c*dest(i, j->row2) + j->s*t;
        } else {
          dest(i, j->row1) = j->c*dest(i, j->row1) + j->s*dest(i, j->row2);
          dest(i, j->row2) = j->c*dest(i, j->row2) - j->s*t;
        }
        if (transpose_) j = j->previous;
        else j = j->next;
      } // while j
    } // for i
  }

  template<typename F>
  void GivensRotations<F>::multiplyFirst
  (const concepts::Matrix<F>& fact, concepts::Matrix<F>& dest) const {
    // check sizes of the matrices
    conceptsAssert(fact.dimX() == dest.dimX(), concepts::Assertion());
    conceptsAssert(fact.dimY() == dest.dimY(), concepts::Assertion());
    conceptsAssert(fact.dimX() == this->dimX(),  concepts::Assertion());
    // loop over columns of dest
    for (uint i = 0; i < dest.dimY(); ++i) {
      const J* j = q_.qend;
      if (transpose_) j = q_.qbegin;
      while (j != 0) {
        conceptsAssert((uint)j->row1 < this->dimX(), concepts::Assertion());
        conceptsAssert((uint)j->row2 < this->dimX(), concepts::Assertion());
        F t = dest(j->row1, i);
        if (transpose_) {
          dest(j->row1, i) = j->c*dest(j->row1, i) + j->s*dest(j->row2, i);
          dest(j->row2, i) = j->c*dest(j->row2, i) - j->s*t;
        } else {
          dest(j->row1, i) = j->c*dest(j->row1, i) - j->s*dest(j->row2, i);
          dest(j->row2, i) = j->c*dest(j->row2, i) + j->s*t;
        }
        if (transpose_) j = j->next;
        else j = j->previous;
      }
    }
  }

  template<typename F>
  std::ostream& GivensRotations<F>::info(std::ostream& os) const {
    return os << "GivensRotations(" << this->dimX() << ", " << q_ << ')';
  }

} // namespace sparseqr

#endif // givensRotations_hh