SVD3x3.hh

/* 
This file is part of the BIAS library (Basic ImageAlgorithmS).

Copyright (C) 2003-2009    (see file CONTACT for details)
  Multimediale Systeme der Informationsverarbeitung
  Institut fuer Informatik
  Christian-Albrechts-Universitaet Kiel


BIAS is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation; either version 2.1 of the License, or
(at your option) any later version.

BIAS is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU Lesser General Public License for more details.

You should have received a copy of the GNU Lesser General Public License
along with BIAS; if not, write to the Free Software
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*/

#ifndef __SVD3x3_hh__
#define __SVD3x3_hh__

#include "SVD.hh"

namespace BIAS {
  /** @class SVD3x3
      @ingroup g_mathalgo
      @brief singular value decomposition for 3x3 matrices

      koeser: There is an example matrix
      F[0][0] =0.00075284066373225704;    
      F[0][1] =-0.099906426438442303;      
      F[0][2] =0.16267766990357244;
      F[1][0] =0.10333163045064053;
      F[1][1] =0.0014727467665194025;
      F[1][2] =0.6803859609102173;
      F[2][0] =-0.16234020366691171;      
      F[2][1] =-0.68087420277635458;     
      F[2][2] =0.0021896762792443826;
      where the svd fails (at least) on p4 systems with gcc 3.3.1 and 3.3.3
      with -O3 (independent of march-flag), but where it does not fail with O2.
      SVD3x3 is implemented as an algorithm, which computes each singular value
      and corresponding vector by iterating until some error measure e becomes
      small compared to a reference value ref:
      if (ref+e==ref) then converged=true;
      With O3 this condition never becomes true for the second singular value,
      because the difference of lhs and rhs stays at 1e-18, while in O2 
      this converges after two iterations to zero.
      We have problems here in the order of machine precision, the question is 
      now whether this different behaviour of O2 and O3 is due to a compiler 
      error or an algorithmic weakness.
      When tracing back the two execution branches of ExampleSVD3x3 (O2 and 
      O3), which uses the above matrix one finds that in the very beginning, 
      s and a[k][i] are absolutely the same before the line:
      s += a[k][i]*a[k][i];
      and after this line s differs by the last digit. Such effects sum up and
      result in the failing of the convergence of the SVD.
      Debugging in O3 is complicated because every cout causes a side-effect 
      and changes execution order and results, so I found that even adding one
      simple cout at some line made the algorithm converge for that example.

      In my opinion this means that the svd algorihtm makes some convergence 
      assumptions about values in the order of machine precision which are not
      always fulfilled, therefore it does not converge. This is just one 
      example which happens with O3 but there may also be one for O2, which we
      have not found ? It may however also be an optimization "error" of gcc.

      @author woelk 08/2004 */
  class BIASMathAlgo_EXPORT SVD3x3 : public SVD
  {
  public:
    SVD3x3();

    /** constructor
        solve the general eigenproblem of the rectangular matrix M
        @author woelk 07/2004 
        @bug does not always converge when general svd does **/
    SVD3x3(const Matrix <double> & M, 
           double ZeroThreshold = DEFAULT_DOUBLE_ZERO_THRESHOLD );

    inline virtual ~SVD3x3() {};
  
    /** set a new matrix and compute its decomposition.
        solves the general eigenproblem of the rectangular matrix M 
        @bug does not always converge when general svd does
        @author woelk 07/2004 **/
    int Compute(const Matrix<double>& M, 
                double ZeroThreshold=DEFAULT_DOUBLE_ZERO_THRESHOLD);

  protected:
    // swaps two singular values in S_ and the according columns/rows in V_/U_
    inline void _Swap(int a, int  b)
    {
      double tmp;
      // swap singular values
      tmp=S_[a]; S_[a]=S_[b]; S_[b]=tmp;
      // swap rows in VT_
      tmp=VT_[a][0]; VT_[a][0]=VT_[b][0]; VT_[b][0]=tmp; 
      tmp=VT_[a][1]; VT_[a][1]=VT_[b][1]; VT_[b][1]=tmp; 
      tmp=VT_[a][2]; VT_[a][2]=VT_[b][2]; VT_[b][2]=tmp;
      // swap columns in U_
      tmp=U_[0][a]; U_[0][a]=U_[0][b]; U_[0][b]=tmp; 
      tmp=U_[1][a]; U_[1][a]=U_[1][b]; U_[1][b]=tmp; 
      tmp=U_[2][a]; U_[2][a]=U_[2][b]; U_[2][b]=tmp;
    }
  }; 


} // namespace


#endif // __SVD3x3_hh__