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

# ifndef _SKIT_SSK_H
# define _SKIT_SSK_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:ssk
NAME:   @code{ssk} - symmetric skyline matrix format 
@clindex ssk
@clindex csr
@clindex vec
@clindex permutation
@toindex spooles, multifrontal solver library
@toindex umfpack, multifrontal solver library
@toindex taucs, out-of-core sparse solver library
@cindex Choleski factorization
@cindex direct solver
@cindex multifrontal factorization
@cindex skyline data structure
@cindex chevron data structure
@cindex Gibbs renumbering
DESCRIPTION:       
 @noindent
 The class implements a symmetric matrix Choleski factorization.
 Let @var{a} be a square invertible matrix in 
 @code{csr} format (@pxref{csr class}).
 @example
     	csr<Float> a;
 @end example
 @noindent
 We get the factorization by:
 @example
     	ssk<Float> m = ldlt(a);
 @end example
 @noindent
 Each call to the direct solver for a*x = b writes either:
 @example
     	vec<Float> x = m.solve(b);
 @end example
 @noindent
 or
 @example
     	m.solve(b,x);
 @end example

DATA STRUCTURE:       
 @noindent
 The storage is either skyline, multi-file or
 chevron.
 This alternative depends upon the configuration
 (@pxref{Installing}).
 The chevron data structure is related to the multifrontal 
 algorithm and is implemented by using the spooles library.
 The multi-file data structure refers to the out-of-core
 algorithm and is implemented by using the taucs library.
 If such a library is not available, the @code{ssk} class
 implements a skyline storage scheme and the profil storage of the 
 @code{ssk} matrix is optimized by optimizing the renumbering,
 by using the Gibbs, Pooles and Stockmeyer algorithm available
 via the @code{permutation} class
 (@pxref{permutation class}).
 @quotation
    Algorithm 582, collected Algorithms from ACM.
    Algorithm appeared in ACM-Trans. Math. Software, vol.8, no. 2,
    Jun., 1982, p. 190.
 @end quotation
 @noindent
 When implementing the skyline data structure,
 we can go back to the @code{csr} format, for the 
 pupose of pretty-printing:
 @example
     	cout << ps << color << logscale << csr<Float>(m);
 @end example

AUTHOR: 
     Pierre Saramito
   | Pierre.Saramito@imag.fr
    LMC-IMAG, 38041 Grenoble cedex 9, France
DATE:   21 january 1997
METHODS: @ssk
End:
*/
# include "rheolef/array.h"
# include "rheolef/csr.h"
# include "rheolef/permutation.h"
namespace rheolef { 
	
#if defined(_RHEOLEF_HAVE_SPOOLES)
template<class T>
class spooles_rep : public occurence {
    	void *_p;
public:
	typedef unsigned int size_type;
	typedef T        element_type;
	spooles_rep ();
	spooles_rep (const csr<element_type>&);
	~spooles_rep();
	size_type nrow () const;
	size_type ncol () const;
	size_type nnz () const;
        void solve(const vec<element_type>& b, vec<element_type>& x) const;
	void factorize_ldlt();
};
//<ssk:
template<class T>
class ssk : smart_pointer<spooles_rep<T> > {
public:
// typedefs:

	typedef          T                          element_type;
	typedef typename spooles_rep<T>::size_type  size_type;

// allocators/deallocators:

	ssk ();
	explicit ssk<T> (const csr<T>&);

// accessors:

	size_type nrow () const;
	size_type ncol () const;
	size_type nnz  () const;
    
// factorisation and solver:

        void solve(const vec<T>& b, vec<T>& x) const;
	vec<T> solve(const vec<T>& b) const;
	void factorize_ldlt();
protected:
	const spooles_rep<T>& data() const {
		return smart_pointer<spooles_rep<T> >::data();
	}
	spooles_rep<T>& data() {
		return smart_pointer<spooles_rep<T> >::data();
	}
};
template <class T>
ssk<T> ldlt (const csr<T>& m);
//>ssk:
// ===============================[ INLINED ]====================================
template <class T>
inline
ssk<T>::ssk ()
 : smart_pointer<spooles_rep<T> > (new_macro(spooles_rep<T>))
{
}
template <class T>
inline
ssk<T>::ssk (const csr<T>& a)
 : smart_pointer<spooles_rep<T> > (new_macro(spooles_rep<T> (a)))
{
}
template <class T>
inline
typename ssk<T>::size_type
ssk<T>::nrow () const
{
    return data().nrow();
}
template <class T>
inline
typename ssk<T>::size_type
ssk<T>::ncol () const
{
    return data().ncol();
}
template <class T>
inline
typename ssk<T>::size_type
ssk<T>::nnz () const
{
    return data().nnz();
}
template <class T>
inline
void
ssk<T>::solve(const vec<T>& b, vec<T>& x) const
{
    data().solve(b,x);
}
template <class T>
inline
void
ssk<T>::factorize_ldlt ()
{
    data().factorize_ldlt();
}

//  ! _RHEOLEF_HAVE_SPOOLES
#elif defined(_RHEOLEF_HAVE_TAUCS)
template<class T>
class taucs_rep : public occurence {
public:
        typedef unsigned int size_type;
        typedef T            element_type;
        taucs_rep ();
        taucs_rep (const csr<element_type>&);
        taucs_rep (const taucs_rep<T>&);
        ~taucs_rep();
        size_type nrow () const;
        size_type ncol () const;
        size_type nnz () const;
        void solve(const vec<element_type>& b, vec<element_type>& x) const;
        void factorize_ldlt();
// data:
protected:
        void *_p;
private:
};
template<class T>
class ssk : smart_pointer<taucs_rep<T> > {
public:
// typedefs:

	typedef taucs_rep<T>::size_type    size_type;
	typedef T                          element_type;

// allocators/deallocators:

	ssk ();
	explicit ssk<T> (const csr<T>&);

// accessors:

	size_type nrow () const;
	size_type ncol () const;
	size_type nnz  () const;
    
// factorisation and solver:

        void solve(const vec<T>& b, vec<T>& x) const;
	vec<T> solve(const vec<T>& b) const;
	void factorize_ldlt();
protected:
	const taucs_rep<T>& data() const {
		return smart_pointer<taucs_rep<T> >::data();
	}
};
template <class T>
ssk<T> ldlt (const csr<T>& m);
// ===============================[ INLINED ]====================================
template <class T>
inline
ssk<T>::ssk ()
 : smart_pointer<taucs_rep<T> > (new_macro(taucs_rep<T>))
{
}
template <class T>
inline
ssk<T>::ssk (const csr<T>& a)
 : smart_pointer<taucs_rep<T> > (new_macro(taucs_rep<T> (a)))
{
}
template <class T>
inline
typename ssk<T>::size_type
ssk<T>::nrow () const
{
    return data().nrow();
}
template <class T>
inline
typename ssk<T>::size_type
ssk<T>::ncol () const
{
    return data().ncol();
}
template <class T>
inline
typename ssk<T>::size_type
ssk<T>::nnz () const
{
    return data().nnz();
}
template <class T>
inline
void
ssk<T>::solve(const vec<T>& b, vec<T>& x) const
{
    data().solve(b,x);
}
template <class T>
inline
void
ssk<T>::factorize_ldlt ()
{
    data().factorize_ldlt();
}
//  ! _RHEOLEF_HAVE_TAUCS
#elif defined(_RHEOLEF_HAVE_UMFPACK)

template<class T>
class umfpack_rep {
    	void *_p;
public:
	typedef unsigned int size_type;
	typedef T        element_type;
	umfpack_rep ();
	umfpack_rep (const csr<element_type>&);
	~umfpack_rep();
	size_type nrow () const;
	size_type ncol () const;
	size_type nnz () const;
        void solve(const vec<element_type>& b, vec<element_type>& x) const;
	void factorize();
};
//<ssk:
template<class T>
class ssk : smart_pointer<umfpack_rep<T> > {
public:
// typedefs:

	typedef          T                          element_type;
	typedef typename umfpack_rep<T>::size_type  size_type;

// allocators/deallocators:

	ssk ();
	explicit ssk<T> (const csr<T>&);

// accessors:

	size_type nrow () const;
	size_type ncol () const;
	size_type nnz  () const;
    
// factorisation and solver:

        void solve(const vec<T>& b, vec<T>& x) const;
	vec<T> solve(const vec<T>& b) const;
	void factorize_ldlt();
	void factorize_lu();
protected:
	const umfpack_rep<T>& data() const {
		return smart_pointer<umfpack_rep<T> >::data();
	}
	umfpack_rep<T>& data() {
		return smart_pointer<umfpack_rep<T> >::data();
	}
};
template <class T>
ssk<T> ldlt (const csr<T>& m);
//>ssk:
// ===============================[ INLINED ]====================================
template <class T>
inline
ssk<T>::ssk ()
 : smart_pointer<umfpack_rep<T> > (new_macro(umfpack_rep<T>))
{
}
template <class T>
inline
ssk<T>::ssk (const csr<T>& a)
 : smart_pointer<umfpack_rep<T> > (new_macro(umfpack_rep<T> (a)))
{
}
template <class T>
inline
typename ssk<T>::size_type
ssk<T>::nrow () const
{
    return data().nrow();
}
template <class T>
inline
typename ssk<T>::size_type
ssk<T>::ncol () const
{
    return data().ncol();
}
template <class T>
inline
typename ssk<T>::size_type
ssk<T>::nnz () const
{
    return data().nnz();
}
template <class T>
inline
void
ssk<T>::solve(const vec<T>& b, vec<T>& x) const
{
    data().solve(b,x);
}
template <class T>
inline
void
ssk<T>::factorize_ldlt ()
{
    data().factorize();
}
template <class T>
inline
void
ssk<T>::factorize_lu ()
{
    data().factorize();
}
template <class T>
inline
ssk<T>
lu (const csr<T>& m)
{
    ssk<T> a(m);
    a.factorize_lu();
    return a;
}

//       ! _RHEOLEF_HAVE_UMFPACK
#else // ! UMF TAUCS SPOOLES...
template<class T>
class ssk {
public:
// allocators/deallocators:

	ssk ();
	explicit ssk (const csr<T>&);

// factorisation LDTt and solver:

        void solve(const vec<T>& b, vec<T>& x) const;
	vec<T> solve(const vec<T>& b) const;
	void factorize_ldlt ();

// accessors:

	unsigned int nrow () const { return _ISKY.size()-1; }
	unsigned int ncol () const { return nrow(); }
	unsigned int nnz () const { return _ASKY.size(); }
	typedef typename Array<unsigned int>::const_iterator const_iterator_isky;
	typedef typename Array<T>::const_iterator            const_iterator_asky;
	const_iterator_isky begin_isky () const { return _ISKY.begin(); }
	const_iterator_isky end_isky () const   { return _ISKY.end(); }
	const_iterator_asky begin_asky () const { return _ASKY.begin(); }
	const_iterator_asky end_asky () const   { return _ASKY.end(); }
    
// output:

	std::ostream& dump(std::ostream&) const;

// internal accessors:
private:
        // index in ASKY of the A(i,j) element, colmin(i) <= j <= i
	unsigned int index_a(unsigned int i, unsigned int j) const
	  { return _ISKY (i+1) - i + j - 1; }
        
	// value of the A(i,j) element, colmin(i) <= j <= i
	const T& operator () (unsigned int i, unsigned int j) const
	  { return _ASKY (index_a(i,j)); }
	T& operator () (unsigned int i, unsigned int j)
	  { return _ASKY (index_a(i,j)); }
	
	// index in ASKY of the first element of the i-th row
	unsigned int index_colmin(unsigned int i) const
	  { return _ISKY (i); }

	// column index j of the element asky(q) in the i-th row
	unsigned int col(unsigned int i, unsigned int q) const
	  { return q - _ISKY (i+1) + i + 1; }

	// column index j of the first element of the i-th row
	unsigned int colmin(unsigned int i) const
	  { return col(i, index_colmin(i)); }

	// index in ASKY of the i-th diagonal element
	unsigned int index_diag(unsigned int i) const
	  { return _ISKY (i+1) - 1; }
// data:
private:
	Array<T>            _ASKY ;
	Array<unsigned int> _ISKY;
	permutation         _p;
};
template<class T>
ssk<T>
ldlt (const csr<T>& m);
#endif // ! _RHEOLEF_HAVE_SPOOLES && !_RHEOLEF_HAVE_TAUCS

template <class T>
inline
vec<T>
ssk<T>::solve (const vec<T>& b) const
{ 
    vec<T> x(b.n());
    solve(b,x);
    return x;
}
template <class T>
inline
ssk<T>
ldlt (const csr<T>& m)
{
    ssk<T> a(m);
    a.factorize_ldlt();
    return a;
}
}// namespace rheolef
#endif /* _SKIT_SSK_H */