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

# ifndef _SKIT_CSR_H
# define _SKIT_CSR_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:csr
NAME: @code{csr} - compressed sparse row matrix format
@clindex csr
@clindex vec
@clindex basic_diag
@clindex eye
@clindex ssk
@clindex asr
@clindex dns
@clindex iorheo
@cindex sparse matrix
@cindex distributed memory 
@fiindex @file{.m} matlab matrix
@fiindex @file{.ps} postscript
@fiindex @file{.hb} Harwell-Boeing matrix
@toindex message passing interface (MPI) library
DESCRIPTION:       
 @noindent
 The class implements a matrix in compressed sparse row format.
 A declaration whithout any parametrers correspond to a matrix with null size:
 @example
     	csr<double> a;
 @end example
 @noindent
 Notes that the constructor can be invocated with an initializer:
 @example
     	csr<double> a = b;
 @end example
 
 @noindent
 Input and output, as usual 
 (@pxref{iorheo class}):
 @example
     	cin  >> a;
     	cout << a;
 @end example
 
 @noindent
 Default is the Harwell-Boeing format.
 Various others formated options are available: 
 matlab and postscript plots.

 @noindent
 Affectation from a scalar
 @example
     	a = 3.14;
 @end example
 
 @noindent
 The matrix/vector multiply:
 @example
     	a*x
 @end example
 @noindent
 and the transposed matrix/ vector multiply:
 @example
     	a.trans_mult(x);
 @end example

 @noindent
 The binary operators are:
 @example
     	a*b, a+b, a-b, lambda*a, a*lambda, a/lambda
 @end example
 @noindent
 The unary operators are sign inversion and transposition:
 @example
     	-a, trans(a);
 @end example
 @noindent
 @c **TODO**:
 The combination with a diagonal matrix
 is not yet completely available. 
 The interface would be something like:
 @example
     	basic_diag<double> d;
     	a+d, d+a, a-d, d-a
     	a*d, d*a, 
     	a/d             // aij/djj
     	left_div(a,d)   // aij/dii
 @end example
 @noindent
 @c **TODO**:
 When applied to the matrix directly:
 this feature is not yet completely available. 
 The interface would be something like:
 @example
     	a += d;		// a := a+d
     	a -= d;		// a := a-d
     	a *= d;		// a := a*d
     	a.left_mult(d);	// a := d*a
     	a /= d;		// aij := aij/djj
     	a.left_div(d);  // aij := aij/dii
 @end example
 @noindent
 @c **TODO**:
 The combination with a constant-identity matrix:
 this feature is not yet completely available. 
 The interface would be something like:
 @example
     	double k;
     	a + k*EYE, k*EYE + a, a - k*EYE, k*EYE - a,
     
     	a += e;
     	a -= e;
 @end example
 @noindent
 Get the lower triangular part:
 @example
     	csr<double> l = tril(a);
 @end example
 @noindent
 Conversely, @code{triu} get the upper triangular part.

 @noindent
 For optimizing the profile storage, I could be convenient to 
 use a renumbering algorithm 
 (see also @ref{ssk class} and @ref{permutation class}).
 @example
     	permutation p = gibbs(a);
     	csr<double> b = perm(a, p, q); // B := P*A*trans(Q)
 @end example
 @noindent
 Horizontal matrix concatenation: 
@tex
   $$ 
       a
       =
       \pmatrix{
	   a_{11} & a_{12}  \cr
       }
   $$
@end tex
 @example
     	a = hcat(a11,a12);
 @end example
 @noindent
 Vertical matrix concatenation:
@tex
 $$
   a=\pmatrix{
       a_{11} \cr
       a_{21} \cr
   }
 $$
@end tex
 @example
     	a = vcat(a11,a21);
 @end example
 @noindent
 Explicit conversion from an associative @code{asr e} matrix writes:
 @example
     	a = csr<double>(e);
 @end example
 @noindent
 from a dense @code{dns m} matrix writes:
 @example
     	a = csr<double>(m);
 @end example
 @noindent
 @c NO MORE AVAILABLE with spooles:
 @c and from a symmetric skyline @code{ssk fa} factorized 
 @c matrix writes:
 @c @example
 @c     	csr<double> fa1 = csr<double>(fa);
 @c @end example

NOTE:
  The @code{csr} class is currently under reconstruction
  for the distributed memory support by using a MPI-based
  implementation.
AUTHOR: 
     Pierre Saramito
   | Pierre.Saramito@imag.fr
    LMC-IMAG, 38041 Grenoble cedex 9, France
DATE:   21 january 1997
End:
*/

# include "rheolef/csrrep.h"

namespace rheolef { 

#ifdef _RHEOLEF_HAVE_EXPRESSION_TEMPLATE
template <class A> class MExpr;
#endif // _RHEOLEF_HAVE_EXPRESSION_TEMPLATE

template <class T>
class csr : public smart_pointer<csrrep<T> > {
public:

// typedefs:

    typedef typename csrrep<T>::size_type            size_type;
    typedef          T                               element_type;
    typedef typename csrrep<T>::iterator_value       iterator_value;
    typedef typename csrrep<T>::iterator_index       iterator_index;
    typedef typename csrrep<T>::const_iterator_value const_iterator_value;
    typedef typename csrrep<T>::const_iterator_index const_iterator_index;

// allocators/deallocators/conversions:

	explicit csr (size_type nrow1 = 0, size_type ncol1 = 0, size_type nnz1 = 0)
	  : smart_pointer<csrrep<T> >(new_macro(csrrep<T>(nrow1,ncol1,nnz1))) {}
	explicit csr (const basic_diag<T>& a);
	explicit csr (const dns<T>& a);
	explicit csr (const asr<T>& a);
#if !defined(_RHEOLEF_HAVE_SPOOLES) && !defined(_RHEOLEF_HAVE_TAUCS)
	explicit csr (const ssk<T>& a);
#endif // _RHEOLEF_HAVE_SPOOLES || _RHEOLEF_HAVE_TAUCS
	explicit csr (const char* filename);

// set pointers to zero and resize nnz=0; conserve nrow,ncol

        void clear();
    
// affectation from a scalar

	csr operator = (const T& lambda);

#ifdef _RHEOLEF_HAVE_EXPRESSION_TEMPLATE

// logical assign, i.e. set row pointers and resize (ind,val) arrays to nnz
// it's the first pass before to store (ind,val) values.

	template <class A>
	void logassign (const A& a) {

            size_type nnzb = 0;
            iterator iter_b = begin();
            iterator last_b = end();
            typedef typename A::const_iterator a_const_iterator;
            a_const_iterator iter_a = a.begin();
            while (iter_b != last_b) {
	        // size required when b is csr
	        size_type row_size = (*iter_a).size();
	        (*iter_b).resize (row_size);  
	        nnzb += row_size;
	        ++iter_b;
	        ++iter_a;
            }  
            (*this).resize(nrow(), ncol(), nnzb);
        }
	template <class A>
	csr(const MExpr<A>& a)
    	: smart_pointer<csrrep<Float> >(new_macro(csrrep<Float>(a.nrow(),a.ncol())))
	{
    	    (*this).logassign(a);
    	    bassigna (begin(), end(), a.begin());
	}
	template <class A>
	csr<T>& operator = (const MExpr<A>& a)
	{
    	    if (nrow() == 0 && ncol() == 0) {
	        (*this).resize (a.nrow(), a.ncol(), 0);
            } else {
                check_mat_length (*this, a);
	    }
            (*this).logassign(a);
            bassigna (begin(), end(), a.begin());
            return *this;
	}
#endif // _RHEOLEF_HAVE_EXPRESSION_TEMPLATE

// blas3 member
        csr<T> left_mult (const basic_diag<T>& d);

// read/write access and info functions
	size_type nrow () const { return data().nrow(); }
        size_type ncol () const { return data().ncol(); }
        size_type nnz () const  { return data().nnz(); }
	void resize (size_type nr, size_type nc, size_type nz)
	    { data().resize(nr, nc, nz); }

        // read direct access
	const Array<T>&       a() const { return data().a(); }
	const Array<size_type>&  ja() const { return data().ja(); }
	const Array<size_type>&  ia() const { return data().ia(); }

	// write direct access
	Array<T>&       a()     { return data().a(); }
	Array<size_type>&  ja() { return data().ja(); }
	Array<size_type>&  ia() { return data().ia(); }

	// slow read access
	T operator() (size_type i, size_type j) const;

// values interface (see also: Array<T>)
	bool any_element_is_negative () const
	    { return a().any_element_is_negative (); }
	bool any_element_is_inf_or_nan () const
	    { return a().any_element_is_inf_or_nan (); }
	bool all_elements_are_int_or_inf_or_nan () const
	    { return a().all_elements_are_int_or_inf_or_nan (); }
	T max_abs () const { return a().max_abs (); }
	T min_abs () const { return a().min_abs (); }
	T max () const { return a().max (); }
	T min () const { return a().min (); }

    // computed assignments
	vec<T> trans_mult (const vec<T>& x) const;

// internal

    // most of operations works only if matrix is sorted by increasing order of column
    // e.g. a+b or a(i,j) -- because unsorted version is less efficace
    // thus, constructors (as read on file) need to know if a matrix is sorted, and sort

	bool is_sorted () const;
	csr<T> sort ();

    // access to rows:
	class const_row {
          public:
	    typedef T element_type;
	    typedef std::pair<size_type,T> value_type;
	    size_type size() const { return (*(POS_IA+1))-(*POS_IA); }
	    size_type n() const { return NCOL; }
	    std::pair<size_type,T> operator() (size_type i) const;
	    explicit const_row(const_iterator_value begin_a, const_iterator_index begin_ja, 
			       const_iterator_index pos_ia, size_type ncol1)
	        : BEGIN_A(begin_a), BEGIN_JA(begin_ja), POS_IA(pos_ia), NCOL(ncol1) {}
          private:
	    const_iterator_value BEGIN_A;
	    const_iterator_index BEGIN_JA;
	    const_iterator_index POS_IA;
	    size_type        NCOL;
          public:
#ifdef TO_CLEAN
	    // random_access_iterator is GNU
            class const_iterator : public std::random_access_iterator<T, ptrdiff_t> {
#else
            class const_iterator {
#endif
              public:
		      std::pair<size_type,T> operator*() { 
		    return std::pair<size_type,T>(*ITER_I, *ITER_V); }
	        void operator++() { 
		    ++ITER_I; ++ITER_V; }
	        bool operator == (const_iterator q) { return ITER_I == q.ITER_I; }
	        bool operator != (const_iterator q) { return !(ITER_I == q.ITER_I); }
                const_iterator(const_iterator_value iter_v, const_iterator_index iter_i)
	          : ITER_I(iter_i), ITER_V(iter_v) {}
              private:
	        const_iterator_index ITER_I;
	        const_iterator_value ITER_V;
            };
	    const_iterator begin() const { return const_iterator(BEGIN_A, BEGIN_JA); }
	    const_iterator end() const { return const_iterator(BEGIN_A+size(), BEGIN_JA+size()); }
        };
#ifdef TO_CLEAN
        class const_iterator : public random_access_iterator<const_row, ptrdiff_t> {
#else
        class const_iterator {
#endif
          public:
            typedef typename csr<T>::const_row row;
	    row operator*() const { return row(ITER_A, ITER_JA, ITER_IA, NCOL); }
	    void operator++() 
	    {
	        size_type nnz_row = *(ITER_IA+1) - *ITER_IA;
	        ITER_A  += nnz_row;
	        ITER_JA += nnz_row;
	        ++ITER_IA;
	    }
	    row operator [] (size_type i) const 
	    {
	      size_type row_start = ITER_IA [i];
	      return row (ITER_A+row_start, ITER_JA+row_start, ITER_IA+i, NCOL);
	    }
	    bool operator == (const_iterator q) { return ITER_IA == q.ITER_IA; }
	    bool operator != (const_iterator q) { return !(ITER_IA == q.ITER_IA); }
	    explicit const_iterator(const_iterator_value iter_a, const_iterator_index iter_ja, 
	                            const_iterator_index iter_ia, size_type ncol1)
	        : ITER_A(iter_a), ITER_JA(iter_ja), ITER_IA(iter_ia), NCOL(ncol1) {}
          private:
	    const_iterator_value ITER_A;
	    const_iterator_index ITER_JA;
	    const_iterator_index ITER_IA;
	    size_type            NCOL;
        };
	const_iterator begin() const
	  { return const_iterator(a().begin(), ja().begin(), ia().begin(), ncol()); }
	const_iterator end() const
	  { return const_iterator(const_iterator_value(0), const_iterator_index(0), ia().end()-1, 0); }
	const_row operator() (size_type i) const // slow read access
	{
          size_type row_start = ia()(i);
          size_type i_row_nnz = ia()(i+1) - row_start;
          const_iterator_value  begin_a  =  a().begin() + row_start;
          const_iterator_index  begin_ja = ja().begin() + row_start;
          const_iterator_index  begin_ia = ia().begin() + i;
          return const_row (begin_a, begin_ja, begin_ia, ncol());
	}
	class row {
          public:
	    typedef std::pair<size_type,T> value_type;
	    typedef typename csr<T>::element_type element_type;
	    void resize(size_type nnz_row) { (*(ITER_IA+1)) = (*ITER_IA) + nnz_row; }
	    size_type size() const { return (*(ITER_IA+1))-(*ITER_IA); }
	    size_type n() const { return NCOL; }
	    void reset();
	    explicit row(iterator_value begin_a, iterator_index begin_ja, 
		iterator_index iter_ia, size_type ncol1)
	      : BEGIN_A(begin_a), BEGIN_JA(begin_ja), ITER_IA(iter_ia), NCOL(ncol1) {}
        private:
	    iterator_value BEGIN_A;
	    iterator_index BEGIN_JA;
	    iterator_index ITER_IA;
	    size_type      NCOL;
        public:
#ifdef TO_CLEAN
            class iterator : public random_access_iterator<T, ptrdiff_t> {
#else
            class iterator {
#endif
              public:
		      std::pair<size_type,T> operator*() { return std::pair<size_type,T>(*ITER_I, *ITER_V); }
	        void operator++() { ++ITER_I; ++ITER_V; }
	        bool operator == (iterator q) { return ITER_I == q.ITER_I; }
	        bool operator != (iterator q) { return !(ITER_I == q.ITER_I); }
#ifndef TOCLEAN
#endif // TOCLEAN
                iterator(iterator_value iter_v, iterator_index iter_i)
		  : ITER_I(iter_i), ITER_V(iter_v) {}
                friend class csr<T>::row;
              private:
	        iterator_index ITER_I;
	        iterator_value ITER_V;
            };
	    iterator begin() { return iterator(BEGIN_A, BEGIN_JA); }
	    iterator end() { return iterator(BEGIN_A+size(), BEGIN_JA+size()); }
	    iterator insert (iterator position, std::pair<size_type,T> iv)
	    { 
    	        *(position.ITER_I) = iv.first; 
    	        *(position.ITER_V) = iv.second;
    	        ++position;
    	        return position;
	    }
        };
#ifdef TO_CLEAN
	class iterator : public random_access_iterator<row, ptrdiff_t> {
#else
	class iterator {
#endif
          public:
            typedef typename csr<T>::row row;
	    row operator*() { return row (ITER_A, ITER_JA, ITER_IA, NCOL); }
	    void operator++() 
	    {
	        size_type nnz_row = *(ITER_IA+1) - *ITER_IA;
	        ITER_A  += nnz_row;
	        ITER_JA += nnz_row;
	        ++ITER_IA;
	    }
	    bool operator == (iterator q) { return ITER_IA == q.ITER_IA; }
	    bool operator != (iterator q) { return !(ITER_IA == q.ITER_IA); }
#ifndef TOCLEAN
#endif // TOCLEAN
	    explicit iterator(iterator_value iter_a, iterator_index iter_ja, iterator_index iter_ia, size_type ncol1)
	        : ITER_A(iter_a), ITER_JA(iter_ja), ITER_IA(iter_ia), NCOL(ncol1) {}
          private:
	    iterator_value ITER_A;
	    iterator_index ITER_JA;
	    iterator_index ITER_IA;
	    size_type      NCOL;
        };
	iterator begin()
	{
          return iterator(a().begin(), ja().begin(), ia().begin(), ncol());
        }
	iterator end()
	{
    	  return iterator(iterator_value(0), iterator_index(0), ia().end()-1,0);
	}
  protected:
	csrrep<T>& data() {
		return smart_pointer<csrrep<T> >::data(); }
	const csrrep<T>& data() const {
		return smart_pointer<csrrep<T> >::data(); }
};
// y := a*x
template <class T> vec<T> operator * (const csr<T>& a, const vec<T>& x);

// transpose, lower and upper triangular parts, gibbs reordering:
template<class T> csr<T> trans (const csr<T>&);
template<class T> csr<T> tril  (const csr<T>&, int k = 0);
template<class T> csr<T> triu  (const csr<T>&, int k = 0);
template<class T> csr<T> perm (const csr<T>& a, const permutation& p, 
                                                const permutation& q);

// the skyline area of a csr matrix
template<class T> typename csr<T>::size_type nnzsky (const csr<T>& a);

// io routines; see also iorheo manipulators
template<class T> std::ostream& operator << (std::ostream&, const csr<T>&);
template<class T> std::istream& operator >> (std::istream&, csr<T>&);

// matlab-like sparse pattern visualization
template <class T> void spy(const csr<T>&);

// concatenations
template<class T> csr<T> hcat (const csr<T>& a1, const csr<T>& a2);
template<class T> csr<T> vcat (const csr<T>& a1, const csr<T>& a2);

// get elem(i): return (i,xi) or (-1, 0.0)
// NOTE: suppose row is sorted by increasing indexes
template <class T>
std::pair<typename csr<T>::size_type,T> 
csr<T>::const_row::operator() (size_type i) const
{
    typename csr<T>::const_row::const_iterator iter = csr<T>::const_row::begin();
    typename csr<T>::const_row::const_iterator last = csr<T>::const_row::end();
    while (iter != last && (*iter).first <= i) {
	if ((*iter).first == i) return *iter;
        ++iter;
    }
    return std::pair<typename csr<T>::size_type,T>((typename csr<T>::size_type)(-1),T());
}
template <class T>
void
csr<T>::row::reset ()
{
    typename csr<T>::row::iterator iter = begin();
    typename csr<T>::row::iterator last = end();
    while (iter != last) {
	insert  (iter, std::pair<typename csr<T>::size_type,T>((*iter).first, T()));
	++iter;
    }
}
}// namespace rheolef
# endif /* _SKIT_CSR_H */