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

# ifndef _SKIT_CSR_ALGORITHM_H
# define _SKIT_CSR_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
/// 
/// =========================================================================
//
// CSR: Compressed Sparse Row format
//
// algorithm-oriented generic library
// inspired from sparskit2 fortran library
//
// triangular part extraction
// and sort by increasing column order
//
// author: Pierre.Saramito@imag.fr
//
// date: 4 september 1997
//
//@!\vfill\listofalgorithms

/*@! 
 \vfill \pagebreak \mbox{} \vfill \begin{algorithm}[h]
  \caption{{\tt tril, triu}: get lower, upper triangular part}
  \begin{algorithmic}
    \NOTE {this algorithm is in-place: {\tt tril} (ia, ja, a, ia, ja, a).}
      Current implementation uses use a predicate operator as parameter.
    \ENDNOTE
    \INPUT {the matrix in CSR format}
      ia(0:nrow), ja(0:nnz-1), a(0:nnz-1)
    \ENDINPUT
    \OUTPUT {non-null elements in selected part}
      il(0:nrow), jl(0:nnzl-1), l(0:nnzl-1)
    \ENDOUTPUT
    \BEGIN 
      q := 0 \\
      \FORTO {i := 0} {nrow-1}
	qold := q \\
        \FORTO {p := ia(i)} {ia(i+1)-1} 
	  \IF {ja(p) $\leq$ i}
	    jl (q) := ja(p) \\
            l  (q) := a(p)  \\
	    q++
          \ENDIF
        \ENDFOR
	il (i) := qold;
      \ENDFOR
      il (n) := q
    \END
 \end{algorithmic} \end{algorithm}
 \vfill
 {\bf Implementation Note} \\
  A {\tt BinaryPredicate} argument is required for the selected part of the matrix.
  Uses {\tt less\_equal<size_t>()} for lower triangular part {\tt tril}, 
  and {\tt greater\_equal<size_t>()} for upper part {\tt triu}.
  See also {\tt <functional.h>} for others usefull predicates.

 {\bf Example} \\
\begin{verbatim}
  template <class Difference>
  class less_equal_k : binary_function<Difference, Difference, bool> {
       Difference k;
     public:
       less_equal_k (Difference k1 = 0) : k(k1) {}
       bool operator()(const Difference& i, const Difference& j) const 
                         { return i-j <= k; }
  };
  trig (ia, ia+nrow+1, ja, a, ia, ja, a, less_equal_k (k))
\end{verbatim}

  The argument k specifies how many diagonals above the main diagonal should 
  also be set to zero. If the value of k is negative, additional elements above
  are also selected.
 \vfill \pagebreak \mbox{} \vfill
*/
namespace rheolef { 
template <
    class InputIterator1,
    class InputIterator2,
    class InputIterator3,
    class OutputIterator1,
    class OutputIterator2,
    class OutputIterator3,
    class BinaryPredicate>
inline
void
trig (
    InputIterator1  iter_ia,
    InputIterator1  last_ia,
    InputIterator2  iter_ja,
    InputIterator3  iter_a,
    OutputIterator1 iter_il,
    OutputIterator2 iter_jl,
    OutputIterator3 iter_l,
    BinaryPredicate select)
{
    size_t nnzl = 0;
    size_t i = 0; 
    size_t p = *iter_ia++;
    while (iter_ia < last_ia) {
        size_t nnzlold = nnzl;
	size_t pf = *iter_ia++;
        while (p < pf) {
            size_t j = *iter_ja++;
            if (select(j, i)) {
	        *iter_jl++ = j;
	        *iter_l++  = *iter_a++;
	        nnzl++;
            }
	    p++;
        }
        *iter_il++ = nnzlold;
        i++;
    }
    *iter_il++ = nnzl;
}
/*@! 
 \vfill \pagebreak \mbox{} \vfill \begin{algorithm}[h]
  \caption{{\tt nnz\_triu, nnz\_tril}: 
	   count non-null in lower, upper triangular part}
  \begin{algorithmic}
    \NOTE {see {\tt tril} and {\tt triu}}
    \ENDNOTE
    \INPUT {the matrix in CSR format (without values)}
      ia(0:nrow), ja(0:nnz-1)
    \ENDINPUT
    \OUTPUT {number of non-null elements in selected part}
      nnzl
    \ENDOUTPUT
    \BEGIN 
      nnzl := 0 \\
      \FORTO{i := 0}{nrow-1}
        \FORTO{p := ia(i)} {ia(i+1)-1}
	  {\em here is lower; for upper, see the generic implementation}
	  \IF {ja(p) $\leq$ i}
	    nnzl++;
          \ENDIF
        \ENDFOR
      \ENDFOR
    \END
 \end{algorithmic} \end{algorithm}
 \vfill
*/
template <
    class InputIterator1,
    class InputIterator2,
    class BinaryPredicate>
inline
size_t
nnz_trig (
    InputIterator1  iter_ia,
    InputIterator1  last_ia,
    InputIterator2  iter_ja,
    BinaryPredicate select)
{
    size_t nnzg = 0;
    size_t p = *iter_ia++;
    size_t i = 0;
    while (iter_ia < last_ia) {
        size_t pf = *iter_ia++;
        while (p < pf) {
	    if (select(*iter_ja++, i))
		nnzg++;
	    ++p;
        }
        ++i;
    }
    return nnzg;
}
/*@! 
 \vfill \pagebreak \mbox{} \vfill \begin{algorithm}[h]
  \caption{{\tt perm}: permutation: $B$ := $P*A*Q^t$, i.e a(p(i), q(j)) :== a(i,j)}
  \begin{algorithmic}
    \NOTE {}
       {\sc warning: } if the original matrix is sorted by increasing column order,
	then the permuted one may not be on return. see also {\tt csr\_sort} \\ \mbox{}
    \ENDNOTE
    \INPUT {the permutations and the CSR matrix}
      p(0:nrow-1), q(0:ncol-1), ia(0:nrow), ja(0:nnz-1), a(0:nnz-1)
    \ENDINPUT
    \OUTPUT {the permuted matrix}
      iao(0:nrow), jao(0:nnz-1), ao(0:nnz-1)
    \ENDOUTPUT
    \BEGIN 
      {\em first pass on ia and iao: get row length from pointers} \\
      iao(p(i)+1) := ia(i+1) - ia(i), \ i := 0\ldots nrow-1 \\ \mbox{} \\

      {\em second pass on iao: get new pointers from lengths} \\
      iao(0) := 0 \\
      iao(i+1) += iao(i), \ i := 0\ldots nrow-1 \\ \mbox{} \\

      {\em third pass on iao, second on ia: copy indexes and values} \\
      \FORTO {i := 0} {nrow-1}
	ro := iao (p(i)) \\
        \FORTO {r := ia(i)} {ia(i+1)-1}
	  jao(ro) := q(ja(r)) \\
	  ao(ro) := a(r)  \\
	  ro++
        \ENDFOR           
      \ENDFOR           
    \END
 \end{algorithmic} \end{algorithm} \vfill
 \vfill \pagebreak \mbox{} \vfill
*/
template <
    class RandomAccessIterator,
    class InputIterator,
    class MutableRandomAccessIterator>
size_t
_length_from_pointers (
    RandomAccessIterator        p,
    InputIterator               iter_ia,
    InputIterator               last_ia,
    MutableRandomAccessIterator rand_iao)
{
    size_t nrow    = 0;
    size_t prev_ia = *iter_ia++;
    while (iter_ia != last_ia) {
	size_t curr_ia = *iter_ia++;
	rand_iao [p[nrow]+1] = curr_ia - prev_ia;
	prev_ia = curr_ia;
	nrow++;
    }
    return nrow;
}
//@!\vfill
template <class MutableForwardIterator>
void
_pointers_from_length (
    MutableForwardIterator        first_iao,
    MutableForwardIterator        last_iao)
{
    size_t prev_iao = (*first_iao++) = 0;
    while (first_iao != last_iao)
	prev_iao = ( (*first_iao++) += prev_iao);
} 
//@! \vfill \pagebreak \mbox{} \vfill
template <
    class RandomAccessIterator1,
    class RandomAccessIterator2,
    class InputIterator1,
    class InputIterator2,
    class InputIterator3,
    class RandomAccessIterator3,
    class MutableRandomAccessIterator1,
    class MutableRandomAccessIterator2>
void
_perm_copy (
    RandomAccessIterator1              p,
    RandomAccessIterator2             q,
    InputIterator1                    iter_ia,
    InputIterator1                    last_ia,
    InputIterator2                    iter_ja,
    InputIterator3                    iter_a,
    RandomAccessIterator3             rand_iao,
    MutableRandomAccessIterator1      rand_jao,
    MutableRandomAccessIterator2      rand_ao)
{
    size_t r = *iter_ia++;
    size_t i = 0;
    while (iter_ia != last_ia) {
	size_t ko = rand_iao [p[i]];
	MutableRandomAccessIterator1 iter_jao = rand_jao + ko;
	MutableRandomAccessIterator2 iter_ao  = rand_ao  + ko;
	size_t last_r = *iter_ia++;
        while (r < last_r) {
	    *iter_jao++ = q[*iter_ja++];
	    *iter_ao++  = *iter_a++;
	    r++;
        }
	i++;
    }
}
//@! \vfill \pagebreak \mbox{} \vfill
template <
    class RandomAccessIterator,
    class ForwardIterator,
    class InputIterator1,
    class InputIterator2,
    class RandomAccessIterator1,
    class RandomAccessIterator2,
    class RandomAccessIterator3>
void
perm (
    RandomAccessIterator  p,
    RandomAccessIterator  q,
    ForwardIterator       first_ia,
    ForwardIterator       last_ia,
    InputIterator1        iter_ja,
    InputIterator2        iter_a,
    RandomAccessIterator1 rand_iao,
    RandomAccessIterator2 rand_jao,
    RandomAccessIterator3 rand_ao)
{
    size_t nrow;
    nrow = _length_from_pointers (p, first_ia, last_ia, rand_iao);
    _pointers_from_length (rand_iao, rand_iao + nrow + 1);
    _perm_copy (p, q, first_ia, last_ia, iter_ja, iter_a, rand_iao, rand_jao, rand_ao);
}
/*@! 
 \vfill \pagebreak \mbox{} \vfill \begin{algorithm}[h]
  \caption{{\tt csr\_is\_sorted}: rows are sorted by increasing column order}
  \begin{algorithmic}
    \INPUT {the CSR matrix without values}
      ia(0:nrow), ja(0:nnz-1)
    \ENDINPUT
    \OUTPUT {boolean }
      is\_sorted
    \ENDOUTPUT
    \BEGIN 
      \FORTO {i := 0} {nrow-1}
        \FORTO {p := ia(i)+1} {ia(i+1)-1}
	  \IF {ja(p) <= ja(p-1)} 
	    is\_sorted := false
          \ENDIF
        \ENFOR
      \ENDFOR           
    \END
 \end{algorithmic} \end{algorithm} \vfill
 \vfill \pagebreak \mbox{} \vfill
*/
template <class InputIterator1, class InputIterator2>
bool
csr_is_sorted (
	InputIterator1 iter_ia,
	InputIterator1 last_ia,
	InputIterator2 iter_ja)
{
    size_t p = *iter_ia++;
    while (iter_ia != last_ia) {
	size_t last_p = *iter_ia++;
	if (p >= last_p) continue;
	p++;
	size_t old_j = *iter_ja++;
	size_t j;
	while (p < last_p) {
	    j = *iter_ja++;
	    p++;
	    if (old_j >= j) {
	      return false;
            }
	}
	old_j = j;
    }
    return true;
}
//@!\vfill
}// namespace rheolef
# endif // _SKIT_CSR_ALGORITHM_H