Re: [eigen] Fast QR for 2x2, 3x3

[ Thread Index | Date Index | More lists.tuxfamily.org/eigen Archives ]


Assuming m_q and m_r are 2x2 matrices to store the results, matrix is the 
input 2x2 matrix which is to be factorized, this code gives the QR 
factorisation:

  if (matrix(1,0) == 0.0) // nothing to do
  {
    m_q = MatrixType::Identity();
    m_r = matrix;
  }
  else if (matrix(0,0) == 0.0) // transposition of rows
  {
    m_q(0,0) = 0.0;
    m_q(0,1) = 1.0;
    m_q(1,0) = 1.0;
    m_q(1,1) = 0.0;
    m_r(0,0) = matrix(1,0);
    m_r(0,1) = matrix(1,1);
    m_r(1,0) = 0.0;
    m_r(1,1) = matrix(0,1);
  }
  else
  {
    // avoid overflow. is this neccesary?
    // if not, just use Scalar er = 1.0/matrix.block<2,1>(0,0).norm();
    Scalar er;
    if (matrix(0,0) >= matrix(1,0))
    {
      Scalar ba = matrix(1,0)/matrix(0,0);
      er = 1.0/matrix(0,0) * sqrt(1+ba*ba);
    }
    else
    {
      Scalar ab = matrix(0,0)/matrix(1,0);
      er = 1.0/matrix(1,0) * sqrt(1+ab*ab);
    }
    m_q(0,0) = er * matrix(0,0);
    m_q(1,1) = er * matrix(0,0);
    m_q(0,1) = -er * matrix(1,0);
    m_q(1,0) = er * matrix(1,0);
    m_r(0,0) = 1.0/er;
    m_r(1,0) = 0.0;
    m_r(0,1) = er * (matrix(0,0)*matrix(0,1) + matrix(1,0)*matrix(1,1));
    m_r(1,1) = er * (matrix(0,0)*matrix(1,1) - matrix(1,0)*matrix(0,1));
  }

A specialisation for 3x3 will follow in a few days.
-- Markus



Mail converted by MHonArc 2.6.19+ http://listengine.tuxfamily.org/