tiny_lu.h

Go to the documentation of this file.

#ifndef _RHEO_TINY_LU_H
#define _RHEO_TINY_LU_H
 
#include "rheolef/tiny_matvec.h"
 
// first step: LU factorization
// ----------------------------
// with partial pivoting
//
// references :
//  P. Lascaux, R. Theodor
//  "Analyse numerique matricielle
//   appliquee a l'art de l'ingenieur",
//  page 242,
//  Masson, 1986
//
 
namespace rheolef { 
template <class T>
void
lu (tiny_matrix<T>& a, tiny_vector<size_t>& piv)
{
        typedef size_t size_type;
        const size_type n = a.nrow();
        if (n == 0) return;
 
        // initialize permutation table
        for (size_type i = 0; i < n; i++)
            piv [i] = i;
 
        // factorize in 'n' steps
        for (size_type k = 0; k < n-1; k++) {
 
            // we search the largest element of th k-th
            // line, that has not yet been pivot-line
            T amax = abs(a(piv[k],k));
            size_type jmax = k;
 
            for (size_type i = k+1; i < n; i++) {
                if (abs(a(piv[i],k)) > amax) {
                    amax = abs(a(piv[i],k));
                    jmax = i;
                }
            }
            // largest element is in piv[jmax] line
            // we permut indexes
            size_type i = piv [k];
            piv [k] = piv [jmax];
            piv [jmax] = i;
            
            // and invert the pivot
            if (1 + a(piv[k],k) == 1) { // a (piv[k],k) < zero machine
                error_macro ("lu: unisolvence failed on pivot " << k);
            }
            T pivinv = 1./a(piv[k],k);
 
            // modify lines that has not yet been
            // pivot-lines
            for (size_type i = k+1; i < n; i++) {
 
                T c = a(piv[i],k) * pivinv;
                a(piv[i],k) = c;
                for (size_type j = k+1; j < n; j++) 
                    a(piv [i],j)  -=  c * a(piv[k],j);
            }
        }
}
// second step: one-column resolution
// ----------------------------------
template <class T>
void
solve (tiny_matrix<T>& a, tiny_vector<size_t>& piv, 
    const tiny_vector<T>& b, tiny_vector<T>& x)
{
        typedef size_t size_type;
        const size_type n = a.nrow();
        if (n == 0) return;
 
        // solve Ly = piv(b); y is stored in x
        for (size_type i = 0; i < n; i++) {
 
            T c = 0;
            for (size_type j = 0; j < i; j++)
 
                c += a(piv[i],j) * x [j];
 
            x [i] = b [piv[i]] - c;
        }
        // solve Ux = y; x contains y as input and x as output
        for (int i = n-1; i >= 0; i--) {
 
            T c = 0;
            for (size_type j = i+1; j < n; j++)
 
                c += a(piv[i],j) * x [j];
 
            x [i] = (x [i] - c) / a(piv[i],i);
        }
}
// ---------------------------------
// third step : matrix inversion
// NOTE: the a matrix is destroyed !
// ---------------------------------
 
template <class T>
void
invert (tiny_matrix<T>& a, tiny_matrix<T>& inv_a)
{
    typedef size_t size_type;
    const size_type n = a.nrow();
 
    // performs the gauss factorization:  M = L.U
    tiny_vector<size_t> piv (n);
    lu (a, piv);
    
    // invert M in B, column by colomn
    tiny_vector<T> column (n);
    tiny_vector<T> x (n);
    inv_a.resize (n,n);
 
    for (size_type j = 0; j < n; j++) {
 
        for (size_type i = 0; i < n; i++) 
            column [i] = 0;
        column [j] = 1;
 
        solve (a, piv, column, x);
 
        for (size_type i = 0; i < n; i++) 
          inv_a (i,j) = x [i];
    }
}
template <class T>
void
put (std::ostream& out, std::string name, const tiny_matrix<T>& a)
{
    typedef size_t size_type;
    out << name << "(" << a.nrow() << "," << a.ncol() << ")" << std::endl;
 
    for (size_type i = 0; i < a.nrow(); i++) {
      for (size_type j = 0; j < a.ncol(); j++) {
         out << name << "(" << i << "," << j << ") = " << a(i,j) << std::endl;
      }
    }
}
}// namespace rheolef
#endif // _RHEO_TINY_LU_H