mlLinearAlgebraTools.h

Go to the documentation of this file.

/*************************************************************************************
**
** Copyright 2007, MeVis Medical Solutions AG
**
** The user may use this file in accordance with the license agreement provided with
** the Software or, alternatively, in accordance with the terms contained in a
** written agreement between the user and MeVis Medical Solutions AG.
**
** For further information use the contact form at https://www.mevislab.de/contact
**
**************************************************************************************/
 
#ifndef ML_LINEAR_ALGEBRA_TOOLS_H
#define ML_LINEAR_ALGEBRA_TOOLS_H
 
 
// Include system independent file and project settings.
#include "mlLinearAlgebraSystem.h"
#include "mlLinearAlgebraDefs.h"
 
#include <mlPrintTemplateErrors.h>
#include <sstream>
 
// 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)
  {
    OBJ_TYPE buf = obj1;
    obj1 = obj2;
    obj2 = buf;
  }
 
  //-------------------------------------------------------------------
  //-------------------------------------------------------------------
  template <class BASE_TYPE>
  BASE_TYPE MLInverseMatHelper(const BASE_TYPE                            &origMat,
                               bool                                       *isInvertible,
                               const typename BASE_TYPE::ComponentType  /*ZeroEpsilon*/,
                               const char * const                          ZeroDetErrString,
                               const BASE_TYPE                            &Identity,
                               const size_t                                Dim
                              )
  {
    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 typename BASE_TYPE::ComponentType avjj = a[j][j];
      if (MLValueIs0WOM(avjj)){
        if (isInvertible == nullptr)
        {
          std::stringstream sstream;
          sstream.precision(4);
          sstream << "Returning identity, matrix was:" << std::endl;
          for (size_t y = 0; y < Dim; ++y)
          {
            for (size_t x = 0; x < Dim; ++x)
            {
              sstream << std::fixed << origMat[x][y] << "\t";
            }
            sstream << std::endl;
          }
          printTemplateError(ZeroDetErrString, ML_BAD_PARAMETER, sstream.str().c_str());
        }
        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 != nullptr){ *isInvertible = true; }
    return b;
  }
 
ML_LA_END_NAMESPACE
 
#endif // __mlLinearAlgebraTools_H