Loading...
Searching...
No Matches
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
rheolef::tiny_matrix
Definition tiny_matvec.h:62
rheolef::tiny_matrix::resize
void resize(size_type nr, size_type nc)
Definition tiny_matvec.h:134
rheolef::tiny_vector
Definition tiny_matvec.h:39
error_macro
#define error_macro(message)
Definition dis_macros.h:49
T
Expr1::float_type T
Definition field_expr.h:230
amax
#define amax(a, b)
Definition geo_partition_scotch.h:124
rheolef
This file is part of Rheolef.
Definition compiler_eigen.h:39
rheolef::solve
void solve(tiny_matrix< T > &a, tiny_vector< size_t > &piv, const tiny_vector< T > &b, tiny_vector< T > &x)
Definition tiny_lu.h:92
rheolef::put
void put(std::ostream &out, std::string name, const tiny_matrix< T > &a)
Definition tiny_lu.h:155
rheolef::lu
void lu(tiny_matrix< T > &a, tiny_vector< size_t > &piv)
Definition tiny_lu.h:41
rheolef::invert
void invert(tiny_matrix< T > &a, tiny_matrix< T > &inv_a)
Definition tiny_lu.h:127