1 #ifndef _RHEO_TINY_LU_H
2 #define _RHEO_TINY_LU_H
24 #include "rheolef/tiny_matvec.h"
60 if (abs(
a(piv[i],k)) >
amax) {
61 amax = abs(
a(piv[i],k));
72 if (1 +
a(piv[k],k) == 1) {
73 error_macro (
"lu: unisolvence failed on pivot " << k);
75 T pivinv = 1./
a(piv[k],k);
81 T c =
a(piv[i],k) * pivinv;
84 a(piv [i],j) -=
c *
a(piv[k],j);
105 c +=
a(piv[i],j) * x [j];
107 x [i] =
b [piv[i]] -
c;
110 for (
int i =
n-1; i >= 0; i--) {
115 c +=
a(piv[i],j) * x [j];
117 x [i] = (x [i] -
c) /
a(piv[i],i);
147 solve (
a, piv, column, x);
158 out <<
name <<
"(" <<
a.nrow() <<
"," <<
a.ncol() <<
")" << std::endl;
162 out <<
name <<
"(" << i <<
"," << j <<
") = " <<
a(i,j) << std::endl;
field::size_type size_type
void resize(size_type nr, size_type nc)
#define error_macro(message)
This file is part of Rheolef.
void put(std::ostream &out, std::string name, const tiny_matrix< T > &a)
void solve(tiny_matrix< T > &a, tiny_vector< size_t > &piv, const tiny_vector< T > &b, tiny_vector< T > &x)
void lu(tiny_matrix< T > &a, tiny_vector< size_t > &piv)
void invert(tiny_matrix< T > &a, tiny_matrix< T > &inv_a)