1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
|
#ifndef __DIA_H__
#define __DIA_H__
#include <algorithm>
/*
* Compute Y += A*X for DIA matrix A and dense vectors X,Y
*
*
* Input Arguments:
* I n_row - number of rows in A
* I n_col - number of columns in A
* I n_diags - number of diagonals
* I L - length of each diagonal
* I offsets[n_diags] - diagonal offsets
* T diags[n_diags,L] - nonzeros
* T Xx[n_col] - input vector
*
* Output Arguments:
* T Yx[n_row] - output vector
*
* Note:
* Output array Yx must be preallocated
* Negative offsets correspond to lower diagonals
* Positive offsets correspond to upper diagonals
*
*/
template <class I, class T>
void dia_matvec(const I n_row,
const I n_col,
const I n_diags,
const I L,
const I offsets[],
const T diags[],
const T Xx[],
T Yx[])
{
for(I i = 0; i < n_diags; i++){
const I k = offsets[i]; //diagonal offset
const I i_start = std::max<I>(0,-k);
const I j_start = std::max<I>(0, k);
const I j_end = std::min<I>(std::min<I>(n_row + k, n_col),L);
const I N = j_end - j_start; //number of elements to process
const T * diag = diags + (npy_intp)i*L + j_start;
const T * x = Xx + j_start;
T * y = Yx + i_start;
for(I n = 0; n < N; n++){
y[n] += diag[n] * x[n];
}
}
}
#endif
|
