librheolef-dev 5.93-2 / usr / include / rheolef / diag.h

# ifndef _SKIT_DIAG_H
# define _SKIT_DIAG_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
/// 
/// =========================================================================
/*Class:basic_diag
NAME: @code{basic_diag} - diagonal matrix 
@clindex basic_diag
@clindex vec
@cindex diagonal matrix
DESCRIPTION:       
 The class implements a diagonal matrix.
 A declaration whithout any parametrers correspond to a null size matrix:
 @example
     	basic_diag<Float> d;
 @end example
 @noindent
 The constructor can be invocated whith a size parameter:
 @example
     	basic_diag<Float> d(n);
 @end example
 @noindent
 or an initialiser, either a vector (@pxref{vec class}):
 @example
     	basic_diag<Float> d = basic_diag(v);
 @end example
 @noindent
 or a csr matrix @pxref{csr class}:
 @example
     	basic_diag<Float> d = basic_diag(a);
 @end example
 @noindent
 The conversion from @code{basic_diag} 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{nrow}
 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. 

NOTE:
 Since the @code{basic_diag} class derives from the 
 @code{vec},  the @code{basic_diag} class
 present also a STL-like interface.
 
AUTHOR: 
    Pierre Saramito
    | Pierre.Saramito@imag.fr
    LMC-IMAG, 38041 Grenoble cedex 9, France
DATE:   28 january 1997

End:
*/
# include "rheolef/vec.h"
# include "rheolef/csr.h"
# include "rheolef/blas1-dense.h"
namespace rheolef { 
//<basic_diag:
template<class T>
class basic_diag : public vec<T> {
public:

// typedefs:

    typedef typename vec<T>::element_type element_type;
    typedef typename vec<T>::size_type    size_type;
    typedef typename vec<T>::iterator     iterator;

// allocators/deallocators:

    explicit basic_diag (size_type sz = 0);
    explicit basic_diag (const vec<T>& u);
    explicit basic_diag (const csr<T>& a);

// assignment:
  
    basic_diag<T> operator = (const T& lambda);

// accessors:
	
    size_type nrow () const { return vec<T>::size(); }
    size_type ncol () const { return vec<T>::size(); }

// basic_diag as a preconditionner: solves D.x=b

    vec<T> solve (const vec<T>& b) const;
    vec<T> trans_solve (const vec<T>& b) const;
};
template <class T> 
basic_diag<T> dcat (const basic_diag<T>& a1, const basic_diag<T>& a2);

template <class T>
basic_diag<T> operator / (const T& lambda, const basic_diag<T>& d);

template <class T>
vec<T>
operator * (const basic_diag<T>& d, const vec<T>& x);

template<class T>
vec<T> left_div (const vec<T>& x, const basic_diag<T>& d);
//>basic_diag:

// =============== inline'd =====================================

template <class T>
inline
basic_diag<T>::basic_diag (size_type sz)
  : vec<T>(sz)
{
}
template <class T>
inline
basic_diag<T>::basic_diag (const vec<T>& u)
  : vec<T>(u)  
{
}
template <class T>
inline
basic_diag<T>::basic_diag (const csr<T>& a)
  : vec<T>(a.nrow())
{
    size_type i = 0;
    iterator iter = vec<T>::begin();
    iterator last = vec<T>::end();
    while (iter < last) {
	*iter = a(i,i);
	++iter;
	++i;
    }
}
template <class T>
inline
basic_diag<T>
basic_diag<T>::operator = (const T& lambda)
{
    fill (vec<T>::begin(), vec<T>::end(), lambda);
    return *this;
}
template <class T>
inline
basic_diag<T>
dcat (const basic_diag<T>& a1, const basic_diag<T>& a2)
{
  basic_diag<T> a(a1.size()+a2.size()) ;
  Cat(a1.begin(),a1.size(),a2.begin(),a2.size(),a.begin()) ; 
  return a ;
}
template <class T>
inline
basic_diag<T>
operator / (const T& lambda, const basic_diag<T>& d)
{
    return basic_diag<T> (lambda/(vec<T>(d)));
}
template <class T>
inline
vec<T>
operator *(const basic_diag<T>& d, const vec<T>& x)
{
    return (vec<T>(d)) * x;
}
template <class T>
inline
vec<T>
basic_diag<T>::solve (const vec<T>& b) const
{
    return b / (vec<T>(*this));
}
template <class T>
inline
vec<T>
basic_diag<T>::trans_solve (const vec<T>& b) const
{
    return b / (vec<T>(*this));
}
template<class T>
inline
vec<T>
left_div (const vec<T>& x, const basic_diag<T>& d)
{
  return x/(vec<T>(d)) ; 
}
}// namespace rheolef
# endif /* _SKIT_DIAG_H */