[Octave-bug-tracker] [bug #29768] logm fix

[Jaroslav Hajek]

Follow-up Comment #21, bug #29768 (project octave):

Hi Marco,

I found out I have already played with this a couple of months before and had
an incomplete Fortran implementation, so I completed and pushed it. This code
stores the rotations in helper vectors so that it can process both matrices
strictly column-wise to be cache-friendly.

I invite you to verify the new function. One thing I would like to note: you
used

      s = sqrt(T(m,m-1)/(T(m-1,m)-T(m,m-1)));
      c = sqrt(T(m-1,m)/(T(m-1,m)-T(m,m-1)));

I do not think, however, that T(m-1,m) > 0 is guaranteed as an output of real
schur, is it? In that case, with c complex,
      G = [c, s;-s', c];
would not be a rotation...

I used instead the formulas

      c = sqrt(abs(T(m-1,m)/(T(m-1,m)-T(m,m-1))));
      s = sqrt(abs(T(m,m-1)/(T(m-1,m)-T(m,m-1)))) * sign(T(m-1,m));

the sign of s should ensure that the conjugate with positive imaginary part
comes first.

I will let schur() and sqrtm() take advantage of the new feature.

    _______________________________________________________

Reply to this item at:

  <http://savannah.gnu.org/bugs/?29768>

_______________________________________________
  Message sent via/by Savannah
  http://savannah.gnu.org/