← Back to branch summary

~ubuntu-branches/ubuntu/trusty/rheolef/trusty

« back to all changes in this revision

Viewing changes to nfem/geo_element/tiny_lu.h

  • Committer: Package Import Robot
  • Author(s): Pierre Saramito
  • Date: 2012-04-06 09:12:21 UTC
  • mfrom: (1.1.5)
  • Revision ID: package-import@ubuntu.com-20120406091221-m58me99p1nxqui49
Tags: 6.0-1
http://bugs.debian.org/667356
http://bugs.debian.org/661689
* New upstream release 6.0 (major changes):
  - massively distributed and parallel support
  - full FEM characteristic method (Lagrange-Gakerkin method) support
  - enhanced users documentation 
  - source code supports g++-4.7 (closes: #667356)
* debian/control: dependencies for MPI distributed solvers added
* debian/rules: build commands simplified
* debian/librheolef-dev.install: man1/* to man9/* added
* debian/changelog: package description rewritted (closes: #661689)

Show diffs side-by-side

added added

removed removed

Lines of Context:
nfem/geo_element/tiny_lu.h
 
1
#ifndef _RHEO_TINY_LU_H
 
2
#define _RHEO_TINY_LU_H
 
3
///
 
4
/// This file is part of Rheolef.
 
5
///
 
6
/// Copyright (C) 2000-2009 Pierre Saramito <Pierre.Saramito@imag.fr>
 
7
///
 
8
/// Rheolef is free software; you can redistribute it and/or modify
 
9
/// it under the terms of the GNU General Public License as published by
 
10
/// the Free Software Foundation; either version 2 of the License, or
 
11
/// (at your option) any later version.
 
12
///
 
13
/// Rheolef is distributed in the hope that it will be useful,
 
14
/// but WITHOUT ANY WARRANTY; without even the implied warranty of
 
15
/// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
16
/// GNU General Public License for more details.
 
17
///
 
18
/// You should have received a copy of the GNU General Public License
 
19
/// along with Rheolef; if not, write to the Free Software
 
20
/// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 
21
/// 
 
22
/// =========================================================================
 
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
 
41
lu (tiny_matrix<T>& a, tiny_vector<size_t>& piv)
 
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
 
92
solve (tiny_matrix<T>& a, tiny_vector<size_t>& piv, 
 
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
 
127
invert (tiny_matrix<T>& a, tiny_matrix<T>& inv_a)
 
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