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

# ifndef _SKIT_CSR_ALGO_APLB_H
# define _SKIT_CSR_ALGO_APLB_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 xply\_size}: sparse size of $z=x\pm y$, where $x,y$ are sparse vectors.}
  \begin{algorithmic}
    \INPUT {sparse vectors patterns}
      jx(0:nnzx-1), jy(0:nnzy-1)
    \ENDINPUT
    \OUTPUT {number of non-null elements in $z=x\pm y$}
      nnzz
    \ENDOUTPUT
    \NOTE {}
      The input sparse vectors may be sorted by increasing index order (see also {\tt csr\_sort}).
    \ENDNOTE
    \BEGIN 
      nnzz := 0 \\
      p := 0 \\
      q := 0 \\
      \WHILE {p $\neq$ nnzx {\bf or } q $\neq$ nnzy}
          ix := {\bf if} p = nnzx {\bf then} $+\infty$ {\bf else} jx(p) {\bf endif} \\
          iy := {\bf if} q = nnzy {\bf then} $+\infty$ {\bf else} jy(q) {\bf endif} \\
	  \IF {ix = jy}
	      p++\\
	      q++
          \ELSIF {ix $<$ iy}
	      p++
          \ELSE
	      q++
          \ENDIF
	  nnz++
      \ENDWHILE
    \END
 \end{algorithmic} \end{algorithm}
 \vfill \pagebreak \mbox{} \vfill
*/
namespace rheolef { 

template <
    class InputIterator1, 
    class InputIterator2,
    class Size>
inline 
Size
xply_size (
    InputIterator1  jx, 
    InputIterator1  last_jx,
    InputIterator2  jy,
    InputIterator2  last_jy,
    const Size&)
{
    Size nnzz = 0;
    const Size infty = std::numeric_limits<Size>::max();
    while (jx != last_jx || jy != last_jy) {

        Size i1 = jx == last_jx ? infty : *jx;
        Size i2 = jy == last_jy ? infty : *jy;

        if (i1 == i2) {
            ++jx;
            ++jy;
        } else if (i1 < i2) {
            ++jx;
        } else {
            ++jy;
        }
	nnzz++;
    }
    return nnzz;
}
/*@! 
 \vfill \pagebreak \mbox{} \vfill \begin{algorithm}[h]
  \caption{{\tt xply}: compute $z=x\pm y$, where $x,y$ are sparse vectors.}
  \begin{algorithmic}
    \INPUT {sparse vectors}
      jx(0:nnzx-1), x(0:nnzx-1), \\
      jy(0:nnzy-1), y(0:nnzy-1)
    \ENDINPUT
    \OUTPUT {the sparse vector $z=x\pm y$}
      jz(0:nnzz-1), z(0:nnzz-1)
    \ENDOUTPUT
    \NOTE {}
	The input sparse vectors may be sorted by increasing index order (see also {\tt csr\_sort}).
    \ENDNOTE
    \BEGIN 
      p := q := nzzz := 0 \\
      \WHILE {p $\neq\infty$ {\bf or } q $\neq\infty$}
          ix := {\bf if} p = $+\infty$ {\bf then} $+\infty$ {\bf else} jx(p) {\bf endif} \\
          iy := {\bf if} q = $+\infty$ {\bf then} $+\infty$ {\bf else} jy(q) {\bf endif} \\
	  \IF {ix = jy}
	      jz(nzzz) := ix \\
	      z(nzzz) := x(p) $\pm$ y(q) \\
	      p++ \\
	      q++
          \ELSIF {ix $<$ iy}
	      jz(nzzz) := ix \\
	      z(nzzz) := x(p) $\pm$ 0\\
	      p++
          \ELSE
	      jz(nzzz) := iy \\
	      z(nzzz) := 0 $\pm$ y(q) \\
	      q++
          \ENDIF
	  nzzz++
      \ENDWHILE
    \END
 \end{algorithmic} \end{algorithm}
 \vfill \pagebreak \mbox{} \vfill
*/
template <
    class BinaryOperation,
    class InputIterator1,
    class InputIterator2,
    class InputIterator3,
    class InputIterator4,
    class OutputIterator1,
    class OutputIterator2,
    class Size>
inline
Size
xply (
    BinaryOperation    binary_op,
    InputIterator1     jx,
    InputIterator1     last_jx,
    InputIterator2     x,
    InputIterator3     jy,
    InputIterator3     last_jy,
    InputIterator4     y,
    OutputIterator1    jz,
    OutputIterator2    z,
    const Size&)
{
    Size nnzz = 0;
    const Size infty = std::numeric_limits<Size>::max();
    while (jx != last_jx || jy != last_jy) {

        Size i1 = jx == last_jx ? infty : *jx;
        Size i2 = jy == last_jy ? infty : *jy;
        
	if (i1 == i2) {
	    *z++ = binary_op(*x++, *y++);
	    *jz++ = i1;
            ++jx;
            ++jy;
        } else if (i1 < i2) {
	    *z++ = binary_op(*x++, 0);
	    *jz++ = i1;
            ++jx;
        } else {
	    *z++ = binary_op(0, *y++);
	    *jz++ = i2;
            ++jy;
        }
	nnzz++;
    }
    return nnzz;
}
/*@! 
 \vfill \pagebreak \mbox{} \vfill \begin{algorithm}[h]
  \caption{{\tt aplb\_size}: sparse size of $c=a\pm b$.}
  \begin{algorithmic}
    \INPUT {sparse matrix patterns}
      ia(0:nrowa), ja(0:nnza-1), \\
      ib(0:nrowb), jb(0:nnzb-1)
    \ENDINPUT
    \OUTPUT {number of non-null elements in $c=a\pm b$}
      nnzc
    \ENDOUTPUT
    \NOTE {} The input sparse matrix may be sorted by increasing column order (see also {\tt csr\_sort}).
    \ENDNOTE
    \BEGIN 
      nnzz := 0 \\
      \FORTO {i := 0} {nrowa-1}
          nnzc += xply\_size (ja(ia(i):ia(i+1)-1), jb(ib(i):ib(i+1)-1))
      \ENDFOR
    \END
 \end{algorithmic} \end{algorithm}
 \vfill \pagebreak \mbox{} \vfill
*/
template <
  class InputIterator1,
  class InputIterator2,
  class InputIterator3,
  class InputIterator4,
  class Size>
Size
aplb_size (
    InputIterator1 ia,
    InputIterator1 last_ia,
    InputIterator2 ja,
    InputIterator3 ib,
    InputIterator4 jb,
    const Size&)
{
    Size nnzc = 0;
    InputIterator2 first_ja = ja + *ia++;
    InputIterator4 first_jb = jb + *ib++;
    while (ia != last_ia) {
        InputIterator2 last_ja = ja + *ia++;
        InputIterator4 last_jb = jb + *ib++;
	nnzc += xply_size (first_ja, last_ja, first_jb, last_jb, Size());
	first_ja = last_ja;
	first_jb = last_jb;
    }
    return nnzc;
}
/*@! 
 \vfill \pagebreak \mbox{} \vfill \begin{algorithm}[h]
  \caption{{\tt aplb}: compute $c=a\pm b$.}
  \begin{algorithmic}
    \INPUT {sparse matrix in CSR format}
      ia(0:nrowa), ja(0:nnza-1), a(0:nnza-1), \\
      ib(0:nrowb), jb(0:nnzb-1), b(0:nnzb-1)
    \ENDINPUT
    \OUTPUT {the result $c=a\pm b$ in CSR format}
      ic(0:nrowa), jc(0:nnzc-1), c(0:nnzc-1)
    \ENDOUTPUT
    \NOTE {}
	The input sparse matrixes may be sorted by increasing column order (see also {\tt csr\_sort}).
    \ENDNOTE
    \BEGIN 
      nnzc := 0 \\
      \FORTO {i := 0}{nrowa-1}
	  jx := ja(ia(i):ia(i+1)-1) \\
	  x := a(ia(i):ia(i+1)-1) \\
	  jb := jb(ib(i):ib(i+1)-1) \\
	  b := b(ib(i):ib(i+1)-1) \\
	  nzzz += xply(jx, x, jy, y, jc(ic(i):ic(i+1)-1), c(ic(i):ic(i+1)-1))
      \ENDFOR
    \END
 \end{algorithmic} \end{algorithm}
 \vfill \pagebreak \mbox{} \vfill
*/
template <
    class BinaryOperation,
    class InputIterator1,
    class InputIterator2,
    class InputIterator3,
    class InputIterator4,
    class InputIterator5,
    class InputIterator6,
    class OutputIterator1,
    class OutputIterator2,
    class OutputIterator3,
    class Size>
Size
aplb (
    BinaryOperation     binary_op,
    InputIterator1      ia,
    InputIterator1      last_ia,
    InputIterator2      ja,
    InputIterator3      a,
    InputIterator4      ib,
    InputIterator5      jb,
    InputIterator6      b,
    OutputIterator1     ic,
    OutputIterator2     jc,
    OutputIterator3     c,
    const Size&)
{
    Size nnzc = *ic++ = 0;
    Size a_offset = *ia++;
    Size b_offset = *ib++;
    InputIterator2  first_ja = ja + a_offset;
    InputIterator3  first_a  =  a + a_offset;
    InputIterator5  first_jb = jb + b_offset;
    InputIterator6  first_b  =  b + b_offset;
    OutputIterator2 first_jc = jc + nnzc;
    OutputIterator3 first_c  =  c + nnzc;
    while (ia != last_ia) {
        a_offset = *ia++;
        b_offset = *ib++;
        InputIterator2 last_ja = ja + a_offset;
        InputIterator5 last_jb = jb + b_offset;
	nnzc += xply (binary_op,
		    first_ja, last_ja, first_a, 
		    first_jb, last_jb, first_b, 
		    first_jc,          first_c, 
		    Size());
	first_ja = last_ja;
	first_a  = a + a_offset;
	first_jb = last_jb;
	first_b  = b + b_offset;
	*ic++ = nnzc;
	first_jc = jc + nnzc;
	first_c  =  c + nnzc;
    }
    return nnzc;
}
//@!\vfill
}// namespace rheolef
# endif // _SKIT_CSR_ALGO_APLB_H