Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
tiny_lu.h
Go to the documentation of this file.
1#ifndef _RHEO_TINY_LU_H
2#define _RHEO_TINY_LU_H
23
24#include "rheolef/tiny_matvec.h"
25
26// first step: LU factorization
27// ----------------------------
28// with partial pivoting
29//
30// references :
31// P. Lascaux, R. Theodor
32// "Analyse numerique matricielle
33// appliquee a l'art de l'ingenieur",
34// page 242,
35// Masson, 1986
36//
37
38namespace rheolef {
39template <class T>
40void
42{
43 typedef size_t size_type;
44 const size_type n = a.nrow();
45 if (n == 0) return;
46
47 // initialize permutation table
48 for (size_type i = 0; i < n; i++)
49 piv [i] = i;
50
51 // factorize in 'n' steps
52 for (size_type k = 0; k < n-1; k++) {
53
54 // we search the largest element of th k-th
55 // line, that has not yet been pivot-line
56 T amax = abs(a(piv[k],k));
57 size_type jmax = k;
58
59 for (size_type i = k+1; i < n; i++) {
60 if (abs(a(piv[i],k)) > amax) {
61 amax = abs(a(piv[i],k));
62 jmax = i;
63 }
64 }
65 // largest element is in piv[jmax] line
66 // we permut indexes
67 size_type i = piv [k];
68 piv [k] = piv [jmax];
69 piv [jmax] = i;
70
71 // and invert the pivot
72 if (1 + a(piv[k],k) == 1) { // a (piv[k],k) < zero machine
73 error_macro ("lu: unisolvence failed on pivot " << k);
74 }
75 T pivinv = 1./a(piv[k],k);
76
77 // modify lines that has not yet been
78 // pivot-lines
79 for (size_type i = k+1; i < n; i++) {
80
81 T c = a(piv[i],k) * pivinv;
82 a(piv[i],k) = c;
83 for (size_type j = k+1; j < n; j++)
84 a(piv [i],j) -= c * a(piv[k],j);
85 }
86 }
87}
88// second step: one-column resolution
89// ----------------------------------
90template <class T>
91void
93 const tiny_vector<T>& b, tiny_vector<T>& x)
94{
95 typedef size_t size_type;
96 const size_type n = a.nrow();
97 if (n == 0) return;
98
99 // solve Ly = piv(b); y is stored in x
100 for (size_type i = 0; i < n; i++) {
101
102 T c = 0;
103 for (size_type j = 0; j < i; j++)
104
105 c += a(piv[i],j) * x [j];
106
107 x [i] = b [piv[i]] - c;
108 }
109 // solve Ux = y; x contains y as input and x as output
110 for (int i = n-1; i >= 0; i--) {
111
112 T c = 0;
113 for (size_type j = i+1; j < n; j++)
114
115 c += a(piv[i],j) * x [j];
116
117 x [i] = (x [i] - c) / a(piv[i],i);
118 }
119}
120// ---------------------------------
121// third step : matrix inversion
122// NOTE: the a matrix is destroyed !
123// ---------------------------------
124
125template <class T>
126void
128{
129 typedef size_t size_type;
130 const size_type n = a.nrow();
131
132 // performs the gauss factorization: M = L.U
133 tiny_vector<size_t> piv (n);
134 lu (a, piv);
135
136 // invert M in B, column by colomn
137 tiny_vector<T> column (n);
138 tiny_vector<T> x (n);
139 inv_a.resize (n,n);
140
141 for (size_type j = 0; j < n; j++) {
142
143 for (size_type i = 0; i < n; i++)
144 column [i] = 0;
145 column [j] = 1;
146
147 solve (a, piv, column, x);
148
149 for (size_type i = 0; i < n; i++)
150 inv_a (i,j) = x [i];
151 }
152}
153template <class T>
154void
155put (std::ostream& out, std::string name, const tiny_matrix<T>& a)
156{
157 typedef size_t size_type;
158 out << name << "(" << a.nrow() << "," << a.ncol() << ")" << std::endl;
159
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;
163 }
164 }
165}
166}// namespace rheolef
167#endif // _RHEO_TINY_LU_H
168
field::size_type size_type
Definition branch.cc:430
void resize(size_type nr, size_type nc)
#define error_macro(message)
Definition dis_macros.h:49
Expr1::float_type T
Definition field_expr.h:230
#define amax(a, b)
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)
Definition tiny_lu.h:92
void put(std::ostream &out, std::string name, const tiny_matrix< T > &a)
Definition tiny_lu.h:155
void lu(tiny_matrix< T > &a, tiny_vector< size_t > &piv)
Definition tiny_lu.h:41
void invert(tiny_matrix< T > &a, tiny_matrix< T > &inv_a)
Definition tiny_lu.h:127