24#include "rheolef/tiny_matvec.h"
56 T amax = abs(a(piv[k],k));
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;
160 for (
size_type i = 0; i < a.nrow(); i++) {
161 for (
size_type j = 0; j < a.ncol(); j++) {
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 solve(tiny_matrix< T > &a, tiny_vector< size_t > &piv, const tiny_vector< T > &b, tiny_vector< T > &x)
void put(std::ostream &out, std::string name, const tiny_matrix< T > &a)
void lu(tiny_matrix< T > &a, tiny_vector< size_t > &piv)
void invert(tiny_matrix< T > &a, tiny_matrix< T > &inv_a)