MeVis/Foundation/Sources/MLLinearAlgebra/mlLinearAlgebraTools.h

Go to the documentation of this file.

// **InsertLicense** code 
//=====================================================================================

//=====================================================================================

#ifndef __mlLinearAlgebraTools_H
#define __mlLinearAlgebraTools_H

// Include system independency file and project settings.
#ifndef __mlLinearAlgebraSystem_H
#include "mlLinearAlgebraSystem.h"
#endif
#ifndef __mlLinearAlgebraDefs_H
#include "mlLinearAlgebraDefs.h"
#endif


// All declarations of this header will be in the ML_LA_NAMESPACE namespace.
ML_LA_START_NAMESPACE
    
  //-------------------------------------------------------------------
  //-------------------------------------------------------------------
  template <class OBJ_TYPE>
  void MLSwap(OBJ_TYPE &obj1, OBJ_TYPE &obj2)
  {
    ML_TRACE_IN_TIME_CRITICAL("mlLinearAlgebraTools.h: MLSwap()");

    OBJ_TYPE buf = obj1;
    obj1 = obj2;
    obj2 = buf;
  }

  //-------------------------------------------------------------------
  //-------------------------------------------------------------------
  template <class BASE_TYPE>
  BASE_TYPE MLInverseMatHelper(const BASE_TYPE                            &origMat,
                               bool                                       *isInvertible, 
                               const ML_TYPENAME BASE_TYPE::ComponentType  /*ZeroEpsilon*/,
                               const char * const                          ZeroDetErrString,
                               const BASE_TYPE                            &Identity,
                               const size_t                                Dim
                              )
  {
    ML_TRACE_IN("mlLinearAlgebraTools.h: MLInverseMatHelper()");

    BASE_TYPE retVal;
    ML_TRY
    {


      BASE_TYPE a(origMat);  // As a evolves from original mat into identity
      BASE_TYPE b(Identity); // b evolves from identity into inverse(a)
      size_t i=0, j=0, i1=0; // Index counters.

      // Loop over columns of a from left to right, eliminating above and below diagonal
      for (j=0; j<Dim; j++){
        // Find largest pivot in column j among rows j..DIM-1
        i1 = j;       // Row with largest pivot candidate
        for (i=j+1; i<Dim; i++){
          if (MLAbs(a[i][j]) > MLAbs(a[i1][j])){
            i1 = i;
          }
        }
        // Swap rows i1 and j in a and b to put pivot on diagonal
        MLSwap(a[i1], a[j]);
        MLSwap(b[i1], b[j]);      

        // Scale row j to have a identity diagonal
        const ML_TYPENAME BASE_TYPE::ComponentType avjj = a[j][j];
        if (MLValueIs0WOM(avjj)){
          if (isInvertible == NULL){
            ML_PRINT_ERROR(ZeroDetErrString,
                           ML_BAD_PARAMETER,
                           "Returning identity");
          }
          else{
            *isInvertible = false;
          }
          return Identity;
        }
        b[j] /= avjj;
        a[j] /= avjj;
      
        // Eliminate off-diagonal elements in col j of a, 
        // applying identical operations to b
        for (i=0; i<Dim; i++){
          if (i!=j){
            b[i] -= a[i][j]*b[j];
            a[i] -= a[i][j]*a[j];
          }
        }
      }

      // Set state and return inverse.
      if (isInvertible != NULL){ *isInvertible = true; }
      retVal = b;
    }
    ML_CATCH_RETHROW;

    return retVal;
  }

ML_LA_END_NAMESPACE

#endif // __mlLinearAlgebraTools_H