Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
dia.h
Go to the documentation of this file.
1# ifndef _SKIT_DIA_H
2# define _SKIT_DIA_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# include "rheolef/vec.h"
24# include "rheolef/csr.h"
25
26namespace rheolef {
27
28/*Class:dia
29NAME: @code{dia} - diagonal matrix
30@clindex dia
31@clindex vec
32@cindex diagonal matrix
33DESCRIPTION:
34 The class implements a diagonal matrix.
35 A declaration without any parametrers correspond to a null size matrix:
36 @example
37 dia<Float> d;
38 @end example
39 @noindent
40 The constructor can be invocated with a @code{ownership} parameter
41 (see @ref{distributor class}):
42 @example
43 dia<Float> d(ownership);
44 @end example
45 @noindent
46 or an initialiser, either a vector (see @ref{vec class}):
47 @example
48 dia<Float> d(v);
49 @end example
50 @noindent
51 or a csr matrix (see @ref{csr class}):
52 @example
53 dia<Float> d(a);
54 @end example
55 @noindent
56 The conversion from @code{dia} to @code{vec} or @code{csr} is explicit.
57
58 @noindent
59 When a diagonal matrix is constructed from a @code{csr} matrix,
60 the definition of the diagonal of matrix is @emph{always} a vector of size
61 @var{row_ownership} which contains the elements in rows 1 to @var{nrow} of
62 the matrix that are contained in the diagonal.
63 If the diagonal element falls outside the matrix,
64 i.e. @var{ncol} < @var{nrow} then it is defined as a zero entry.
65PRECONDITIONER INTERFACE:
66 The class presents a preconditioner interface,
67 as the @ref{solver class},
68 so that it can be used as preconditioner to the iterative solvers
69 suite (see @ref{cg algorithm}).
70End:
71*/
72//<dia:
73template<class T, class M = rheo_default_memory_model>
74class dia : public vec<T,M> {
75public:
76
77// typedefs:
78
80 typedef typename vec<T,M>::iterator iterator;
82
83// allocators/deallocators:
84
85 explicit dia (const distributor& ownership = distributor(),
86 const T& init_val = std::numeric_limits<T>::max());
87
88 explicit dia (const vec<T,M>& u);
89 explicit dia (const csr<T,M>& a);
90 dia<T,M>& operator= (const T& lambda);
91
92// preconditioner interface: solves d*x=b
93
94 vec<T,M> solve (const vec<T,M>& b) const;
95 vec<T,M> trans_solve (const vec<T,M>& b) const;
96};
97template <class T, class M>
98dia<T,M> operator/ (const T& lambda, const dia<T,M>& d);
99
100template <class T, class M>
101vec<T,M> operator* (const dia<T,M>& d, const vec<T,M>& x);
102//>dia:
103
104// =============== inline'd =====================================
105
106template <class T, class M>
107inline
108dia<T,M>::dia (const distributor& ownership, const T& init_val)
109 : vec<T,M>(ownership, init_val)
110{
111}
112template <class T, class M>
113inline
115 : vec<T,M>(u)
116{
117}
118template <class T, class M>
119inline
121 : vec<T,M>(a.row_ownership())
122{
123 size_type i = 0;
124 typename csr<T,M>::const_iterator dia_ia = a.begin();
125 for (iterator iter = vec<T,M>::begin(), last = vec<T,M>::end(); iter < last; ++iter, ++i) {
126 T a_ii = 0;
127 for (typename csr<T,M>::const_data_iterator p = dia_ia[i]; p < dia_ia[i+1]; p++) {
128 size_type j = (*p).first;
129 if (j == i) { a_ii = (*p).second; break; }
130 }
131 *iter = a_ii;
132 }
133}
134template <class T, class M>
135inline
138{
139 std::fill (vec<T,M>::begin(), vec<T,M>::end(), lambda);
140 return *this;
141}
142template <class T, class M>
143inline
145operator/ (const T& lambda, const dia<T,M>& d)
146{
147 return dia<T,M> (lambda/vec<T,M>(d));
148}
149template <class T, class M>
150inline
151vec<T,M>
152operator* (const dia<T,M>& d, const vec<T,M>& x)
153{
154 return vec<T,M>(d) * x;
155}
156template <class T, class M>
157inline
158vec<T,M>
159dia<T,M>::solve (const vec<T,M>& b) const
160{
161 return b / vec<T,M>(*this);
162}
163template <class T, class M>
164inline
167{
168 return b / vec<T,M>(*this);
169}
170
171}// namespace rheolef
172# endif /* _SKIT_DIA_H */
see the csr page for the full documentation
Definition csr.h:317
vec< T, M >::size_type size_type
Definition dia.h:79
vec< T, M > solve(const vec< T, M > &b) const
Definition dia.h:159
vec< T, M > trans_solve(const vec< T, M > &b) const
Definition dia.h:166
dia(const distributor &ownership=distributor(), const T &init_val=std::numeric_limits< T >::max())
Definition dia.h:108
vec< T, M >::iterator iterator
Definition dia.h:80
vec< T, M >::const_iterator const_iterator
Definition dia.h:81
dia< T, M > & operator=(const T &lambda)
Definition dia.h:137
see the distributor page for the full documentation
Definition distributor.h:69
see the vec page for the full documentation
Definition vec.h:79
base::iterator iterator
Definition vec.h:91
base::size_type size_type
Definition vec.h:86
base::const_iterator const_iterator
Definition vec.h:92
Expr1::float_type T
Definition field_expr.h:230
This file is part of Rheolef.
dia< T, M > operator/(const T &lambda, const dia< T, M > &d)
Definition dia.h:145
csr< T, sequential > operator*(const T &lambda, const csr< T, sequential > &a)
Definition csr.h:437
Definition sphere.icc:25
Definition leveque.h:25
Expr1::memory_type M