← Back to branch summary

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

« back to all changes in this revision

Viewing changes to nfem/basis/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/basis/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