librheolef-dev 6.5-1+b1 / usr / include / rheolef / dia.h

# 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 whithout any parametrers correspond to a null size matrix:
  @example
     	dia<Float> d;
  @end example
  @noindent
  The constructor can be invocated whith 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{pcg 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);

// preconditionner 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 */