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

# ifndef _SKIT_CSR_ALGO_AMUB_H
# define _SKIT_CSR_ALGO_AMUB_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
/// 
/// =========================================================================
//
// CSR: Compressed Sparse Row format
//
// algorithm-oriented generic library
// inspired from sparskit2 fortran library
//  
// author: Pierre.Saramito@imag.fr
//  
// date: 12 november 1997
//  
//@!\vfill\listofalgorithms
/*@! 
 \vfill \pagebreak \mbox{} \vfill \begin{algorithm}[h]
  \caption{{\tt xmub\_size}: sparse size of $y=x^T*b$, where $x$ is a sparse vector.}
  \begin{algorithmic}
    \INPUT {the sparse matrix and vector patterns}
      jx(0:nnzx-1), ib(0:nrowb), jb(0:nnzb-1)
    \ENDINPUT
    \OUTPUT {number of non-null elements in $y=x^T*b$}
      nnzy
    \ENDOUTPUT
    \BEGIN 
      set {\cal y} := $\emptyset$ \\
      \FORTO {p := 0} {nnzx-1}
        i := jx(p) \\
        \FORTO {q := ib(i)} {ib(i+1)-1} 
          {\cal y} := {\cal y} $\cup$ \{ jb(q) \}
        \ENDFOR
      \ENDFOR
      nnzy := {\rm card} ({\cal y})
    \END
 \end{algorithmic} \end{algorithm}
 \vfill \pagebreak \mbox{} \vfill
*/
namespace rheolef { 

template <
    class InputIterator1,
    class InputIterator2,
    class RandomAcessIterator,
    class Size>
inline
Size
xmub_size (
    InputIterator1       jx,
    InputIterator1       last_jx,
    RandomAcessIterator  ib,
    InputIterator2       jb,
    const Size&)
{
    std::set<Size,std::less<Size> > y;
    while (jx != last_jx) {
        Size i = *jx++;
	InputIterator2 jb1 = jb + ib[i];
	InputIterator2 jb2 = jb + ib[i+1];
	typename std::set<Size,std::less<Size> >::iterator iter_y = y.begin();
	while (jb1 != jb2)
	    iter_y = y.insert(iter_y, *jb1++);
    }
    return y.size();
}
/*@! 
 \vfill \pagebreak \mbox{} \vfill \begin{algorithm}[h]
  \caption{{\tt xmub}: compute $y=x^T*b$ where $x$ is a sparse vector.}
  \begin{algorithmic}
    \INPUT {the sparse vector and matrix}
      jx(0:nnzx-1), x(0:nnzx-1), ib(0:nrowb), jb(0:nnzb-1), b(0:nnzb-1)
    \ENDINPUT
    \OUTPUT {the sparse vector $y=x^T*b$}
      jy(0:nnzy-1), y(0:nnzy-1)
    \ENDOUTPUT
    \NOTE {complexity in ${\cal O}(nnzy \log (nnzy))$}
       Could be decreased in ${\cal O}(nnzy)$, with an additional
       working array of size max(nrowa,nrowb).
       Usefull during the sparse product $a*b$.
       See sparskit2 implementation.
    \ENDNOTE
    \BEGIN
      map y := $\emptyset$ \\
      \FORTO {p := 0} {nnzx-1}
	i := jx(p) \\
        \FORTO {q := ib(i)} {ib(i+1)-1}
            y := y + \{ (jb(q), x(p)*b(q)) \}
        \ENDFOR
      \ENDFOR
      \FORTO {r := 0} {card(y)-1}
	  jy (r) := y(r).first \\
	  y (r) := y(r).second
      \ENDFOR
    \END
 \end{algorithmic} \end{algorithm}
 \vfill \pagebreak \mbox{} \vfill
*/
template <
    class InputIterator1,
    class InputIterator2,
    class InputIterator3,
    class InputIterator4,
    class RandomAccessIterator,
    class OutputIterator1,
    class OutputIterator2,
    class Size,
    class T>
inline
Size
xmub (
    InputIterator1       jx,
    InputIterator1       last_jx,
    InputIterator2       x,
    RandomAccessIterator ib,
    InputIterator3       jb,
    InputIterator4       b,
    OutputIterator1      jy,
    OutputIterator2      y,
    const Size&,
    const T&)
{
    typename std::map<Size,T,std::less<Size> > yy;
    typename std::map<Size,T,std::less<Size> >::iterator pred_y, iter_y;
    while (jx != last_jx) {
	Size i = *jx++;
	T lambda = *x++;
	iter_y = yy.begin();
	InputIterator3 jb1 = jb + ib[i];
	InputIterator3 jb2 = jb + ib[i+1];
	InputIterator4  b1 =  b + ib[i];
	while (jb1 != jb2) {
	    Size j = *jb1++;
	    T    yval = lambda*(*b1++);
            pred_y = yy.find(j);
	    if (pred_y == yy.end()) {
		yy.insert(iter_y, std::pair<const Size, T>(j, yval));
	    } else {
		iter_y = pred_y;
		(*iter_y).second += yval;
	    }
        }
    }
    for (iter_y = yy.begin(); iter_y != yy.end(); iter_y++) {
	*jy++ = (*iter_y).first;
	*y++ = (*iter_y).second;
    }
    return yy.size();
}
/*@! 
 \vfill \pagebreak \mbox{} \vfill \begin{algorithm}[h]
  \caption{{\tt amub\_size}: sparse size of $c=a*b$.}
  \begin{algorithmic}
    \INPUT {the sparse matrix patterns}
      ia(0:nrowa), ja(0:nnza-1), ib(0:nrowb), jb(0:nnzb-1)
    \ENDINPUT
    \OUTPUT {sparse size of $c=a*b$}
      nnzc
    \ENDOUTPUT
    \BEGIN 
      nnzc := 0 \\
      \FORTO {i := 0} {nrowa-1}
        nnzc += xmua\_size(ja(ia(i):ia(i+1)-1), ib, jb)
      \ENDFOR
    \END
 \end{algorithmic} \end{algorithm}
 \vfill \pagebreak \mbox{} \vfill
*/
template <
    class InputIterator1,
    class InputIterator2, 
    class InputIterator3,
    class RandomAcessIterator,
    class Size>
Size
amub_size (
    InputIterator1      ia,
    InputIterator1      last_ia,
    InputIterator2      ja,
    RandomAcessIterator ib,
    InputIterator3      jb,
    const Size&)
{
    Size nnzc = 0;
    InputIterator2 first_ja = ja + *ia++;
    while (ia != last_ia) {
        InputIterator2 last_ja = ja + *ia++;
	nnzc += xmub_size (first_ja, last_ja, ib, jb, nnzc);
	first_ja = last_ja;
    }
    return nnzc;
}
/*@! 
 \vfill \pagebreak \mbox{} \vfill \begin{algorithm}[h]
  \caption{{\tt amub}: compute $c=a*b$.}
  \begin{algorithmic}
    \INPUT {the sparse matrix $a,b$ and $ic$ pointers}
      ia(0:nrowa), ja(0:nnza-1), a(0:nnza-1), 
      ib(0:nrowb), jb(0:nnzb-1), b(0:nnzb-1)
      ic(0:nrowa)
    \ENDINPUT
    \OUTPUT {the sparse matrix $c=a*b$}
      jc(0:nnzc-1), c(0:nnzc-1)
    \ENDOUTPUT
    \BEGIN
      ic (0) := 0;
      \FORTO {i := 0} {nrowa-1}
	 jx := ja(ia(i):ia(i+1)-1) \\
	 x := a(ia(i):ia(i+1)-1) \\
	 jz := jc(ic(i):ic(i+1)-1) \\
	 z := c(ic(i):ic(i+1)-1) \\
         ic(i+1) := ic(i) + xmub (jx, x, ib, ib, b, jz, z)
      \ENDFOR
    \END
 \end{algorithmic} \end{algorithm}
 \vfill \pagebreak \mbox{} \vfill
*/
template <
    class InputIterator1,
    class InputIterator2,
    class InputIterator3,
    class InputIterator4,
    class InputIterator5,
    class RandomAcessIterator,
    class OutputIterator1,
    class OutputIterator2,
    class OutputIterator3,
    class Size,
    class T>
Size
amub (
    InputIterator1       ia,
    InputIterator1       last_ia,
    InputIterator2       ja,
    InputIterator3       a,
    RandomAcessIterator  ib,
    InputIterator4       jb,
    InputIterator5       b,
    OutputIterator1      ic,
    OutputIterator2      jc,
    OutputIterator3      c,
    const Size&,
    const T&)
{
    Size a_offset = *ia++;
    InputIterator2 first_ja = ja + a_offset;
    InputIterator3 first_a  = a + a_offset;
    Size nnzc = *ic++ = 0;
    OutputIterator2 first_jc = jc + nnzc;
    OutputIterator3 first_c  = c + nnzc;
    while (ia != last_ia) {
        a_offset = *ia++;
        InputIterator2 last_ja = ja + a_offset;
	nnzc += xmub (first_ja, last_ja, first_a, ib, jb, b, 
	      first_jc, first_c, Size(), T());
        *ic++ = nnzc;
	first_ja = last_ja;
	first_a  = a + a_offset;
	first_jc = jc + nnzc;
	first_c  =  c + nnzc;
    }
    return nnzc;
}
//@!\vfill
}// namespace rheolef
#endif // _SKIT_CSR_ALGO_AMUB_H