rheolef  6.5
csr_cumul_trans_mult.h
Go to the documentation of this file.
1 # ifndef _SKIT_CSR_CUMUL_TRANS_MULT_H
2 # define _SKIT_CSR_CUMUL_TRANS_MULT_H
3 //@!\vfill\listofalgorithms
4 /*@!
5  \vfill \pagebreak \mbox{} \vfill \begin{algorithm}[h]
6  \caption{{\tt trans\_mult}: sparse matrix $y += a^T*x$, where $x,y$ are dense vectors.}
7  \begin{algorithmic}
8  \INPUT {sparse matrix a and dense vector x}
9  ia(0:nrowa), ja(0:nnza-1), a(0:nnza-1),
10  x(0:nrowa)
11  \ENDINPUT
12  \OUTPUT {number of non-null elements in $z=x\pm y$}
13  y(0:ncola)
14  \ENDOUTPUT
15  \NOTE {}
16  The $y$ vector may be set to zero before the call in order to
17  compute $y := a^T*x$
18  \ENDNOTE
19  \BEGIN
20  \FORTO {i := 0}{nrowa-1}
21  \FORTO {p := ia(i)}{ia(i+1)-1}
22  y(ja(p)) += a(p) * x(i)
23  \ENDFOR
24  \ENDFOR
25  \END
26  \end{algorithmic} \end{algorithm}
27  \vfill \pagebreak \mbox{} \vfill
28 */
29 namespace rheolef {
30 
31 template <
32  class InputIterator1,
33  class InputIterator3,
34  class SetOperator,
35  class RandomAccessMutableIterator>
36 void
38  InputIterator1 ia,
39  InputIterator1 last_ia,
40  InputIterator3 x,
41  SetOperator set_op, // set_op: += or -= but not = because may cumul
42  RandomAccessMutableIterator y)
43 {
44  typedef typename std::iterator_traits<InputIterator1>::value_type InputIterator2;
45  typedef typename std::iterator_traits<RandomAccessMutableIterator>::value_type T;
46  InputIterator2 a = (*ia++);
47  while (ia != last_ia) {
48  T xi = *x++;
49  InputIterator2 last_a = (*ia++);
50  while (a != last_a) {
51  set_op (y [(*a).first], (*a).second * xi);
52  ++a;
53  }
54  }
55 }
56 //@!\vfill
57 }// namespace rheolef
58 # endif // _SKIT_CSR_CUMUL_TRANS_MULT_H
59