librheolef-dev 5.93-2 / usr / include / rheolef / blas-algorithm.h

#ifndef _SKIT_BLAS_ALGORITHM_H
#define _SKIT_BLAS_ALGORITHM_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
/// 
/// =========================================================================
//
// Algorithms for mat & vec expressions
//
// author: Pierre.Saramito@imag.fr
//
// date: 11 march 1997
//
# include "rheolef/skitbase.h"
namespace rheolef { 

// ============================[ CHECK SIZE ]====================================
template <class Vec1, class Vec2>
inline
void
check_length (const Vec1& x, const Vec2& y)
{
    check_macro (x.n() == y.n(), "incompatible: vec("
	<< x.n() << ") vec("
	<< y.n() << ") combination");
}
template <class A, class X>
inline
void
check_amulx_length (const A& a, const X& x)
{
    check_macro (a.ncol() == x.n(), "incompatible csr("
	<< a.nrow() << "," << a.ncol() << ")*vec("
	<< x.n()                       << ")");
}
// =================[ USEFULL CLASS FUNCTIONS ]==================================
template <class T1, class T2 = T1, class T3 = skit_promote(T1,T2)>
struct add_op : std::binary_function<T1, T2, T3> {
    T3 operator() (const T1& x, const T2& y) const { return x + y; }
};
template <class T1, class T2 = T1, class T3 = skit_promote(T1,T2)>
struct sub_op : std::binary_function<T1, T2, T3> {
    T3 operator() (const T1& x, const T2& y) const { return x - y; }
};
template <class T1, class T2 = T1, class T3 = skit_promote(T1,T2)>
struct mul_op : std::binary_function<T1, T2, T3> {
    T3 operator() (const T1& x, const T2& y) const { return x * y; }
};
template <class T1, class T2 = T1, class T3 = skit_promote(T1,T2)>
struct div_op : std::binary_function<T1, T2, T3> {
    T3 operator() (const T1& x, const T2& y) const { return x / y; }
};
template <class T1, class T2 = T1>
struct add_assign : std::binary_function<T1, T2, T1> {
    T1 operator()(T1& x, const T2& y) const { return x += y; }
};
template <class T1, class T2 = T1>
struct sub_assign : std::binary_function<T1, T2, T1> {
    T1 operator()(T1& x, const T2& y) const { return x -= y; }
};
template <class T1, class T2>
struct mul_assign : std::binary_function<T1, T2, T1> {
    T1 operator()(T1& x, const T2& y) const { return x *= y; }
};
template <class T1, class T2 = T1>
struct div_assign : std::binary_function<T1, T2, T1> {
    T1 operator()(T1& x, const T2& y) const { return x /= y; }
};
// z := op(x)
template 
   <class OutputIterator,
    class Operation,
    class InputIterator>
inline
void
zassignopx (
    OutputIterator iter_z,
    OutputIterator last_z,
    Operation      op,
    InputIterator  iter_x)
{
    while (iter_z != last_z) {
       (*iter_z) = op (*iter_x);
       ++iter_x;
       ++iter_z;
    }
}
// z := op(x,y)
template 
   <class OutputIterator,
    class BinaryOperation,
    class InputIterator1,
    class InputIterator2>
inline
void
zassignxopy (
    OutputIterator  iter_z,
    OutputIterator  last_z,
    BinaryOperation binary_op,
    InputIterator1  iter_x,
    InputIterator2  iter_y)
{
    while (iter_z != last_z) {
       (*iter_z) = binary_op (*iter_x, *iter_y);
       ++iter_x;
       ++iter_y;
       ++iter_z;
    }
}
// z op= x
template
   <class OutputIterator,
    class InputIterator,
    class ComputedAssignment>
inline
void
zopassignx (
    ComputedAssignment comp_assign,
    OutputIterator     iter_z,
    OutputIterator     last_z,
    InputIterator      iter_x)
{
    while (iter_z != last_z) {
    	comp_assign (*iter_z, *iter_x);
	++iter_z;
	++iter_x;
    }
}
// z op= lambda
template
   <class OutputIterator,
    class T,
    class ComputedAssignment>
inline
void
zopassignr (
    ComputedAssignment comp_assign,
    OutputIterator     iter_z,
    OutputIterator     last_z,
    const T&           lambda)
{
    while (iter_z != last_z) {
    	comp_assign (*iter_z, lambda);
	++iter_z;
    }
}
// ==============================[ SPARSE DOT ]==================================
// dot(sx,y); used in a*x
template <
    class ValueRandomConstIterator,
    class Pair,
    class T>
inline
void
sxdoty_one_step_cumul (
    T&                       sum,
    Pair                     sx_pair,
    ValueRandomConstIterator rand_y) 
{
    sum += (sx_pair.second) * (rand_y [sx_pair.first]);
}
template <
    class ValueRandomConstIterator, 
    class PairInputIterator, 
    class T>
T
sxdoty (
    PairInputIterator        iter_sx, 
    PairInputIterator        last_sx, 
    ValueRandomConstIterator rand_y, 
    const T&)
{
    T s = 0;
    while (iter_sx != last_sx) {
	sxdoty_one_step_cumul (s, *iter_sx, rand_y);
	++iter_sx;
    }
    return s;
}
// =======================[ SPARSE ASSIGNMENT ]==================================
// sz = op(sx); op = ident, -, abs,... unary operator
template<
    class PairInput,
    class Operation,
    class PairOutput>
inline
PairOutput
pair_apply (
    PairInput   sx_pair,
    Operation   op,
    PairOutput*)
{
    typedef typename Operation::result_type T;
    return PairOutput(sx_pair.first, op(sx_pair.second));
}
// sz := op(sx)
template<
    class Svec, 
    class Operation,
    class PairInputIterator>
void 
szassignopsx (
    Svec&              sz, 
    Operation          op,
    PairInputIterator  iter_sx, 
    PairInputIterator  last_sx)
{
    sz.reset();
    typename Svec::iterator prec_sz = sz.begin();
    typedef typename Svec::value_type pair_type;
    while (iter_sx != last_sx) {
        prec_sz = sz.insert (prec_sz, 
 	    pair_apply (*iter_sx, op, (pair_type*)(0)));
        ++iter_sx;
    }
}
// ==========================[ BLAS1 SPARSE +- ]=================================
// sx+-sy; one step logical
template <
    class PairInputIterator1, 
    class PairInputIterator2>
inline
void
sxaddopsy_one_step_iter (
    PairInputIterator1& iter_sx,
    PairInputIterator1  last_sx,
    PairInputIterator2& iter_sy,
    PairInputIterator2  last_sy)
{
    Index ind_sx = iter_sx == last_sx ? max_value((Index)(0)) : (*iter_sx).first;
    Index ind_sy = iter_sy == last_sy ? max_value((Index)(0)) : (*iter_sy).first;
    assert_macro (iter_sx != last_sx || iter_sy != last_sy, "past end.");
    if (ind_sx == ind_sy) {
        ++iter_sx;
        ++iter_sy;
    } else if (ind_sx < ind_sy) {
        ++iter_sx;
    } else {
        ++iter_sy;
    }
}
// sx+-sy; logical
template <
    class PairInputIterator1, 
    class PairInputIterator2>
Index
sxaddopsy_size (
    PairInputIterator1  iter_sx, 
    PairInputIterator1  last_sx,
    PairInputIterator2  iter_sy,
    PairInputIterator2  last_sy)
{
    Index size = 0;
    while (iter_sx != last_sx || iter_sy != last_sy) {				

	sxaddopsy_one_step_iter (iter_sx, last_sx, iter_sy, last_sy);
	size++;
    }
    return size;
}
template <
    class PairInputIterator1, 
    class PairInputIterator2,
    class T1,
    class T2,
    class T3>
inline 
Index
sxopsy_size (
    const add_op<T1, T2, T3>&,
    PairInputIterator1         iter_sx, 
    PairInputIterator1         last_sx,
    PairInputIterator2         iter_sy,
    PairInputIterator2         last_sy)
{
    return sxaddopsy_size (iter_sx, last_sx, iter_sy, last_sy);
}
template <
    class PairInputIterator1, 
    class PairInputIterator2,
    class T1,
    class T2,
    class T3>
inline 
Index
sxopsy_size (
    const sub_op<T1, T2, T3>&,
    PairInputIterator1         iter_sx, 
    PairInputIterator1         last_sx,
    PairInputIterator2         iter_sy,
    PairInputIterator2         last_sy)
{
    return sxaddopsy_size (iter_sx, last_sx, iter_sy, last_sy);
}
// sx+-sy; one full step (value and iter)
template <
    class BinaryOperation,
    class PairInputIterator1,
    class PairInputIterator2,
    class PairInput1,
    class PairInput2,
    class Index,
    class T3>
inline
std::pair<Index, T3>
sxaddopsy_one_step (
    BinaryOperation     binary_op,
    PairInputIterator1& iter_sx, 
    PairInputIterator1  last_sx,
    PairInput1          pair_sx,
    PairInputIterator2& iter_sy, 
    PairInputIterator2  last_sy,
    PairInput2          pair_sy,
    std::pair<Index, T3>*)
{
    typedef typename BinaryOperation::first_argument_type  T1;
    typedef typename BinaryOperation::second_argument_type T2;
    Index ind_sx = iter_sx == last_sx ? max_value((Index)(0)) : pair_sx.first;
    Index ind_sy = iter_sy == last_sy ? max_value((Index)(0)) : pair_sy.first;
    assert_macro (iter_sx != last_sx || iter_sy != last_sy, "past end.");
    if (ind_sx == ind_sy) {
	    std::pair<Index, T3> res (ind_sx, binary_op (pair_sx.second, pair_sy.second));
        ++iter_sx;
        ++iter_sy;
        return res;
    } else if (ind_sx < ind_sy) {
	    std::pair<Index, T3> res(ind_sx, binary_op (pair_sx.second, T2()));
        ++iter_sx;
        return res;
    } else {
	std::pair<Index, T3> res(ind_sy, binary_op (T1(), pair_sy.second));
	++iter_sy;
        return res;
    }
}
template <
    class PairInputIterator1,
    class PairInputIterator2,
    class T1,
    class T2,
    class T3>
inline
std::pair<Index, T3>
sxopsy_one_step (
    const add_op<T1, T2, T3>& add,
    PairInputIterator1&       iter_sx, 
    PairInputIterator1        last_sx,
    PairInputIterator2&       iter_sy, 
    PairInputIterator2        last_sy)
{
    return sxaddopsy_one_step (add, iter_sx, last_sx, *iter_sx,
				    iter_sy, last_sy, *iter_sy, 
				    (std::pair<Index, T3>*)(0));
}
template <
    class PairInputIterator1,
    class PairInputIterator2,
    class T1,
    class T2,
    class T3>
inline
std::pair<Index, T3>
sxopsy_one_step (
    const sub_op<T1, T2, T3>& sub,
    PairInputIterator1&       iter_sx, 
    PairInputIterator1        last_sx,
    PairInputIterator2&       iter_sy, 
    PairInputIterator2        last_sy)
{
    return sxaddopsy_one_step (sub, iter_sx, last_sx, *iter_sx,
				    iter_sy, last_sy, *iter_sy,
				    (std::pair<Index, T3>*)(0));
}
// sx +- sy; complete loop
template <
    class Svec,
    class BinaryOperation,
    class PairInputIterator1,
    class PairInputIterator2>
void
sxaddopsy_ (
    Svec&              sz,
    BinaryOperation    binary_op,
    PairInputIterator1 iter_sx,
    PairInputIterator1 last_sx,
    PairInputIterator2 iter_sy,
    PairInputIterator2 last_sy)
{
    typedef typename Svec::value_type Pair3;
    typedef typename Svec::iterator   PairIterator3;
    PairIterator3 prec_sz = sz.begin();
    while (iter_sx != last_sx && iter_sy != last_sy) {
	prec_sz = sz.insert (prec_sz,
	  sxaddopsy_one_step (binary_op, 
	    iter_sx, last_sx, *iter_sx,
	    iter_sy, last_sy, *iter_sy,
	    (Pair3*)(0)));
    }
    typedef typename std::iterator_traits<PairInputIterator2>::value_type Pair2;
    typedef typename Pair2::first_type T2;
    while (iter_sx != last_sx) {
        Pair3 res ((*iter_sx).first, binary_op ((*iter_sx).second, T2()));
	prec_sz = sz.insert (prec_sz, res);
        ++iter_sx;
    }
    typedef typename std::iterator_traits<PairInputIterator1>::value_type pair1_type;
    typedef typename pair1_type::first_type T1;
    while (iter_sy != last_sy) {
        Pair3 res ((*iter_sy).first, binary_op (T1(), (*iter_sy).second));
	prec_sz = sz.insert (prec_sz, res);
        ++iter_sy;
    }
}
template <
    class Svec,
    class PairInputIterator1,
    class PairInputIterator2,
    class T1,
    class T2,
    class T3>
inline
void
sxopsy (
    Svec&                     sz,
    const add_op<T1, T2, T3>& add,
    PairInputIterator1        iter_sx,
    PairInputIterator1        last_sx,
    PairInputIterator2        iter_sy,
    PairInputIterator2        last_sy)
{
    sxaddopsy_ (sz, add, iter_sx, last_sx, iter_sy, last_sy);
}
template <
    class Svec,
    class PairInputIterator1,
    class PairInputIterator2,
    class T1,
    class T2,
    class T3>
inline
void
sxopsy (
    Svec& sz,
    const sub_op<T1, T2, T3>& sub,
    PairInputIterator1        iter_sx,
    PairInputIterator1        last_sx,
    PairInputIterator2        iter_sy,
    PairInputIterator2        last_sy)
{
    sxaddopsy_ (sz, sub, iter_sx, last_sx, iter_sy, last_sy);
}
template <
    class PairInputIterator1,
    class PairInputIterator2,
    class T1,
    class T2,
    class T3>
inline
void
sxopsy_init (
    const add_op<T1, T2, T3>&,
    PairInputIterator1        iter_sx,
    PairInputIterator1        last_sx,
    PairInputIterator2        iter_sy,
    PairInputIterator2        last_sy)
{
    /* for compatibility with mul */
}
template <
    class PairInputIterator1,
    class PairInputIterator2,
    class T1,
    class T2,
    class T3>
inline
void
sxopsy_init (
    const sub_op<T1, T2, T3>&,
    PairInputIterator1        iter_sx,
    PairInputIterator1        last_sx,
    PairInputIterator2        iter_sy,
    PairInputIterator2        last_sy)
{
    /* for compatibility with mul */
}
}// namespace rheolef
#endif // _SKIT_BLAS_ALGORITHM_H