Loading...
Searching...
No Matches
dia.h
Go to the documentation of this file.
# ifndef _SKIT_DIA_H
# define _SKIT_DIA_H
//
// This file is part of Rheolef.
//
// Copyright (C) 2000-2009 Pierre Saramito <Pierre.Saramito@imag.fr>
//
// Rheolef is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// Rheolef is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with Rheolef; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
//
// =========================================================================
# include "rheolef/vec.h"
# include "rheolef/csr.h"
namespace rheolef {
/*Class:dia
NAME: @code{dia} - diagonal matrix
@clindex dia
@clindex vec
@cindex diagonal matrix
DESCRIPTION:
The class implements a diagonal matrix.
A declaration without any parametrers correspond to a null size matrix:
@example
dia<Float> d;
@end example
@noindent
The constructor can be invocated with a @code{ownership} parameter
(see @ref{distributor class}):
@example
dia<Float> d(ownership);
@end example
@noindent
or an initialiser, either a vector (see @ref{vec class}):
@example
dia<Float> d(v);
@end example
@noindent
or a csr matrix (see @ref{csr class}):
@example
dia<Float> d(a);
@end example
@noindent
The conversion from @code{dia} to @code{vec} or @code{csr} is explicit.
@noindent
When a diagonal matrix is constructed from a @code{csr} matrix,
the definition of the diagonal of matrix is @emph{always} a vector of size
@var{row_ownership} which contains the elements in rows 1 to @var{nrow} of
the matrix that are contained in the diagonal.
If the diagonal element falls outside the matrix,
i.e. @var{ncol} < @var{nrow} then it is defined as a zero entry.
PRECONDITIONER INTERFACE:
The class presents a preconditioner interface,
as the @ref{solver class},
so that it can be used as preconditioner to the iterative solvers
suite (see @ref{cg algorithm}).
End:
*/
//<dia:
template<class T, class M = rheo_default_memory_model>
class dia : public vec<T,M> {
public:
// typedefs:
typedef typename vec<T,M>::size_type size_type;
typedef typename vec<T,M>::iterator iterator;
typedef typename vec<T,M>::const_iterator const_iterator;
// allocators/deallocators:
explicit dia (const distributor& ownership = distributor(),
const T& init_val = std::numeric_limits<T>::max());
explicit dia (const vec<T,M>& u);
explicit dia (const csr<T,M>& a);
dia<T,M>& operator= (const T& lambda);
// preconditioner interface: solves d*x=b
vec<T,M> solve (const vec<T,M>& b) const;
vec<T,M> trans_solve (const vec<T,M>& b) const;
};
template <class T, class M>
dia<T,M> operator/ (const T& lambda, const dia<T,M>& d);
template <class T, class M>
vec<T,M> operator* (const dia<T,M>& d, const vec<T,M>& x);
//>dia:
// =============== inline'd =====================================
template <class T, class M>
inline
dia<T,M>::dia (const distributor& ownership, const T& init_val)
: vec<T,M>(ownership, init_val)
{
}
template <class T, class M>
inline
dia<T,M>::dia (const vec<T,M>& u)
: vec<T,M>(u)
{
}
template <class T, class M>
inline
dia<T,M>::dia (const csr<T,M>& a)
: vec<T,M>(a.row_ownership())
{
size_type i = 0;
typename csr<T,M>::const_iterator dia_ia = a.begin();
for (iterator iter = vec<T,M>::begin(), last = vec<T,M>::end(); iter < last; ++iter, ++i) {
T a_ii = 0;
for (typename csr<T,M>::const_data_iterator p = dia_ia[i]; p < dia_ia[i+1]; p++) {
size_type j = (*p).first;
if (j == i) { a_ii = (*p).second; break; }
}
*iter = a_ii;
}
}
template <class T, class M>
inline
dia<T,M>&
dia<T,M>::operator= (const T& lambda)
{
std::fill (vec<T,M>::begin(), vec<T,M>::end(), lambda);
return *this;
}
template <class T, class M>
inline
dia<T,M>
operator/ (const T& lambda, const dia<T,M>& d)
{
return dia<T,M> (lambda/vec<T,M>(d));
}
template <class T, class M>
inline
vec<T,M>
operator* (const dia<T,M>& d, const vec<T,M>& x)
{
return vec<T,M>(d) * x;
}
template <class T, class M>
inline
vec<T,M>
dia<T,M>::solve (const vec<T,M>& b) const
{
return b / vec<T,M>(*this);
}
template <class T, class M>
inline
vec<T,M>
dia<T,M>::trans_solve (const vec<T,M>& b) const
{
return b / vec<T,M>(*this);
}
}// namespace rheolef
# endif /* _SKIT_DIA_H */
rheolef::csr
see the csr page for the full documentation
Definition csr.h:317
rheolef::dia
Definition dia.h:74
rheolef::dia::size_type
vec< T, M >::size_type size_type
Definition dia.h:79
rheolef::dia::solve
vec< T, M > solve(const vec< T, M > &b) const
Definition dia.h:159
rheolef::dia::trans_solve
vec< T, M > trans_solve(const vec< T, M > &b) const
Definition dia.h:166
rheolef::dia::dia
dia(const distributor &ownership=distributor(), const T &init_val=std::numeric_limits< T >::max())
Definition dia.h:108
rheolef::dia::iterator
vec< T, M >::iterator iterator
Definition dia.h:80
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::distributor
see the distributor page for the full documentation
Definition distributor.h:69
rheolef::vec
see the vec page for the full documentation
Definition vec.h:79
rheolef::vec::iterator
base::iterator iterator
Definition vec.h:91
rheolef::vec::size_type
base::size_type size_type
Definition vec.h:86
rheolef::vec::const_iterator
base::const_iterator const_iterator
Definition vec.h:92
rheolef::vec< T, rheo_default_memory_model >::vec
vec(const vec< T, rheo_default_memory_model > &)
Definition vec.h:168
T
Expr1::float_type T
Definition field_expr.h:230
rheolef
This file is part of Rheolef.
Definition compiler_eigen.h:39
rheolef::operator/
dia< T, M > operator/(const T &lambda, const dia< T, M > &d)
Definition dia.h:145
rheolef::operator*
csr< T, sequential > operator*(const T &lambda, const csr< T, sequential > &a)
Definition csr.h:437
d
size_t d
Definition neumann-laplace-lambda.cc:28
lambda
Definition yield_slip_circle.h:34
p
Definition sphere.icc:25
u
Definition leveque.h:25
M
Expr1::memory_type M
Definition vec_expr_v2.h:416