Rheolef  7.1
an efficient C++ finite element environment
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 
26 namespace rheolef {
27 
28 /*Class:dia
29 NAME: @code{dia} - diagonal matrix
30 @clindex dia
31 @clindex vec
32 @cindex diagonal matrix
33 DESCRIPTION:
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.
65 PRECONDITIONER 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}).
70 End:
71 */
72 //<dia:
73 template<class T, class M = rheo_default_memory_model>
74 class dia : public vec<T,M> {
75 public:
76 
77 // typedefs:
78 
79  typedef typename vec<T,M>::size_type size_type;
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 };
97 template <class T, class M>
98 dia<T,M> operator/ (const T& lambda, const dia<T,M>& d);
99 
100 template <class T, class M>
101 vec<T,M> operator* (const dia<T,M>& d, const vec<T,M>& x);
102 //>dia:
103 
104 // =============== inline'd =====================================
105 
106 template <class T, class M>
107 inline
108 dia<T,M>::dia (const distributor& ownership, const T& init_val)
109  : vec<T,M>(ownership, init_val)
110 {
111 }
112 template <class T, class M>
113 inline
115  : vec<T,M>(u)
116 {
117 }
118 template <class T, class M>
119 inline
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 }
134 template <class T, class M>
135 inline
136 dia<T,M>&
138 {
139  std::fill (vec<T,M>::begin(), vec<T,M>::end(), lambda);
140  return *this;
141 }
142 template <class T, class M>
143 inline
144 dia<T,M>
145 operator/ (const T& lambda, const dia<T,M>& d)
146 {
147  return dia<T,M> (lambda/vec<T,M>(d));
148 }
149 template <class T, class M>
150 inline
151 vec<T,M>
152 operator* (const dia<T,M>& d, const vec<T,M>& x)
153 {
154  return vec<T,M>(d) * x;
155 }
156 template <class T, class M>
157 inline
158 vec<T,M>
159 dia<T,M>::solve (const vec<T,M>& b) const
160 {
161  return b / vec<T,M>(*this);
162 }
163 template <class T, class M>
164 inline
165 vec<T,M>
167 {
168  return b / vec<T,M>(*this);
169 }
170 
171 }// namespace rheolef
172 # endif /* _SKIT_DIA_H */
rheolef::dia
Definition: dia.h:74
rheolef::vec::const_iterator
base::const_iterator const_iterator
Definition: vec.h:92
mkgeo_ball.b
int b
Definition: mkgeo_ball.sh:152
rheolef::operator*
csr< T, sequential > operator*(const T &lambda, const csr< T, sequential > &a)
Definition: csr.h:437
rheolef::distributor
see the distributor page for the full documentation
Definition: distributor.h:62
rheolef::dia::trans_solve
vec< T, M > trans_solve(const vec< T, M > &b) const
Definition: dia.h:166
mkgeo_ball.d
int d
Definition: mkgeo_ball.sh:154
rheolef::operator/
dia< T, M > operator/(const T &lambda, const dia< T, M > &d)
Definition: dia.h:145
rheolef::vec
see the vec page for the full documentation
Definition: vec.h:79
p
Definition: sphere.icc:25
rheolef::dia::size_type
vec< T, M >::size_type size_type
Definition: dia.h:79
rheolef::vec::iterator
base::iterator iterator
Definition: vec.h:91
rheolef::dia::dia
dia(const distributor &ownership=distributor(), const T &init_val=std::numeric_limits< T >::max())
Definition: dia.h:108
a
Definition: diffusion_isotropic.h:25
rheolef::csr
see the csr page for the full documentation
Definition: csr.h:317
rheolef::dia::iterator
vec< T, M >::iterator iterator
Definition: dia.h:80
rheolef
This file is part of Rheolef.
Definition: compiler_eigen.h:37
u
Definition: leveque.h:25
rheolef::dia::const_iterator
vec< T, M >::const_iterator const_iterator
Definition: dia.h:81
rheolef::dia::operator=
dia< T, M > & operator=(const T &lambda)
Definition: dia.h:137
rheolef::dia::solve
vec< T, M > solve(const vec< T, M > &b) const
Definition: dia.h:159
rheolef::vec::size_type
base::size_type size_type
Definition: vec.h:86
M
Expr1::memory_type M
Definition: vec_expr_v2.h:416
T
Expr1::float_type T
Definition: field_expr.h:261
lambda
Definition: yield_slip_circle.h:34