Rheolef  7.1
an efficient C++ finite element environment
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 
38 namespace rheolef {
39 template <class T>
40 void
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 // ----------------------------------
90 template <class T>
91 void
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 
125 template <class T>
126 void
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 }
153 template <class T>
154 void
155 put (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:425
void resize(size_type nr, size_type nc)
Definition: tiny_matvec.h:134
size_t size_type
Definition: basis_get.cc:76
#define error_macro(message)
Definition: dis_macros.h:49
Expr1::float_type T
Definition: field_expr.h:261
#define amax(a, b)
This file is part of Rheolef.
void put(std::ostream &out, std::string name, const tiny_matrix< T > &a)
Definition: tiny_lu.h:155
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 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