Repository navigation
Expand file tree
/
Copy pathMatrixInitOp.hpp
More file actions
225 lines (196 loc) · 7.08 KB
/
Copy pathMatrixInitOp.hpp
File metadata and controls
225 lines (196 loc) · 7.08 KB
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
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
/*
//@HEADER
// ************************************************************************
//
// HPCCG: Simple Conjugate Gradient Benchmark Code
// Copyright (2006) Sandia Corporation
//
// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
// license for use of this work by or on behalf of the U.S. Government.
//
// This library is free software; you can redistribute it and/or modify
// it under the terms of the GNU Lesser General Public License as
// published by the Free Software Foundation; either version 2.1 of the
// License, or (at your option) any later version.
//
// This library 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
// Lesser General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
// USA
// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
//
// ************************************************************************
//@HEADER
*/
#ifndef _MatrixInitOp_hpp_
#define _MatrixInitOp_hpp_
#include <simple_mesh_description.hpp>
#include <box_utils.hpp>
#include <CSRMatrix.hpp>
#include <ELLMatrix.hpp>
#include <algorithm>
template<typename GlobalOrdinal>
void sort_if_needed(GlobalOrdinal* list,
GlobalOrdinal list_len)
{
bool need_to_sort = false;
for(GlobalOrdinal i=list_len-1; i>=1; --i) {
if (list[i] < list[i-1]) {
need_to_sort = true;
break;
}
}
if (need_to_sort) {
std::sort(list,list+list_len);
}
}
template<typename MatrixType>
struct MatrixInitOp {
};
template<>
struct MatrixInitOp<miniFE::CSRMatrix<MINIFE_SCALAR,MINIFE_LOCAL_ORDINAL,MINIFE_GLOBAL_ORDINAL> > {
MatrixInitOp(const std::vector<MINIFE_GLOBAL_ORDINAL>& rows_vec,
const std::vector<MINIFE_LOCAL_ORDINAL>& row_offsets_vec,
const std::vector<int>& row_coords_vec,
int global_nx, int global_ny, int global_nz,
MINIFE_GLOBAL_ORDINAL global_n_rows,
const miniFE::simple_mesh_description<MINIFE_GLOBAL_ORDINAL>& input_mesh,
miniFE::CSRMatrix<MINIFE_SCALAR,MINIFE_LOCAL_ORDINAL,MINIFE_GLOBAL_ORDINAL>& matrix)
: rows(&rows_vec[0]),
row_offsets(&row_offsets_vec[0]),
row_coords(&row_coords_vec[0]),
global_nodes_x(global_nx),
global_nodes_y(global_ny),
global_nodes_z(global_nz),
global_nrows(global_n_rows),
mesh(&input_mesh),
dest_rows(&matrix.rows[0]),
dest_rowoffsets(&matrix.row_offsets[0]),
dest_cols(&matrix.packed_cols[0]),
dest_coefs(&matrix.packed_coefs[0]),
n(matrix.rows.size())
{
if (matrix.packed_cols.capacity() != matrix.packed_coefs.capacity()) {
std::cout<<"Warning, packed_cols.capacity ("<<matrix.packed_cols.capacity()<<") != "
<< "packed_coefs.capacity ("<<matrix.packed_coefs.capacity()<<")"<<std::endl;
}
size_t nnz = row_offsets_vec[n];
if (matrix.packed_cols.capacity() < nnz) {
std::cout<<"Warning, packed_cols.capacity ("<<matrix.packed_cols.capacity()<<") < "
" nnz ("<<nnz<<")"<<std::endl;
}
matrix.packed_cols.resize(nnz);
matrix.packed_coefs.resize(nnz);
dest_rowoffsets[n] = nnz;
}
typedef MINIFE_GLOBAL_ORDINAL GlobalOrdinalType;
typedef MINIFE_LOCAL_ORDINAL LocalOrdinalType;
typedef MINIFE_SCALAR ScalarType;
const GlobalOrdinalType* rows;
const LocalOrdinalType* row_offsets;
const int* row_coords;
int global_nodes_x;
int global_nodes_y;
int global_nodes_z;
GlobalOrdinalType global_nrows;
GlobalOrdinalType* dest_rows;
LocalOrdinalType* dest_rowoffsets;
GlobalOrdinalType* dest_cols;
ScalarType* dest_coefs;
int n;
const miniFE::simple_mesh_description<GlobalOrdinalType>* mesh;
inline void operator()(int i)
{
dest_rows[i] = rows[i];
int offset = row_offsets[i];
dest_rowoffsets[i] = offset;
int ix = row_coords[i*3];
int iy = row_coords[i*3+1];
int iz = row_coords[i*3+2];
GlobalOrdinalType nnz = 0;
for(int sz=-1; sz<=1; ++sz) {
for(int sy=-1; sy<=1; ++sy) {
for(int sx=-1; sx<=1; ++sx) {
GlobalOrdinalType col_id =
miniFE::get_id<GlobalOrdinalType>(global_nodes_x, global_nodes_y, global_nodes_z,
ix+sx, iy+sy, iz+sz);
if (col_id >= 0 && col_id < global_nrows) {
GlobalOrdinalType col = mesh->map_id_to_row(col_id);
dest_cols[offset+nnz] = col;
dest_coefs[offset+nnz] = 0;
++nnz;
}
}
}
}
sort_if_needed(&dest_cols[offset], nnz);
}
};
template<>
struct MatrixInitOp<miniFE::ELLMatrix<MINIFE_SCALAR,MINIFE_LOCAL_ORDINAL,MINIFE_GLOBAL_ORDINAL> > {
MatrixInitOp(const std::vector<MINIFE_GLOBAL_ORDINAL>& rows_vec,
const std::vector<MINIFE_LOCAL_ORDINAL>& /*row_offsets_vec*/,
const std::vector<int>& row_coords_vec,
int global_nx, int global_ny, int global_nz,
MINIFE_GLOBAL_ORDINAL global_n_rows,
c
6D1C
onst miniFE::simple_mesh_description<MINIFE_GLOBAL_ORDINAL>& input_mesh,
miniFE::ELLMatrix<MINIFE_SCALAR,MINIFE_LOCAL_ORDINAL,MINIFE_GLOBAL_ORDINAL>& matrix)
: rows(&rows_vec[0]),
row_coords(&row_coords_vec[0]),
global_nodes_x(global_nx),
global_nodes_y(global_ny),
global_nodes_z(global_nz),
global_nrows(global_n_rows),
mesh(&input_mesh),
dest_rows(&matrix.rows[0]),
dest_cols(&matrix.cols[0]),
dest_coefs(&matrix.coefs[0]),
n(matrix.rows.size()),
ncols_per_row(matrix.num_cols_per_row)
{
}
typedef MINIFE_GLOBAL_ORDINAL GlobalOrdinalType;
typedef MINIFE_LOCAL_ORDINAL LocalOrdinalType;
typedef MINIFE_SCALAR ScalarType;
const GlobalOrdinalType* rows;
const int* row_coords;
int global_nodes_x;
int global_nodes_y;
int global_nodes_z;
GlobalOrdinalType global_nrows;
GlobalOrdinalType* dest_rows;
GlobalOrdinalType* dest_cols;
ScalarType* dest_coefs;
int n;
int ncols_per_row;
const miniFE::simple_mesh_description<GlobalOrdinalType>* mesh;
inline void operator()(int i)
{
dest_rows[i] = rows[i];
int offset = i*ncols_per_row;
int ix = row_coords[i*3];
int iy = row_coords[i*3+1];
int iz = row_coords[i*3+2];
GlobalOrdinalType nnz = 0;
for(int sz=-1; sz<=1; ++sz)
for(int sy=-1; sy<=1; ++sy)
for(int sx=-1; sx<=1; ++sx) {
GlobalOrdinalType col_id =
miniFE::get_id<GlobalOrdinalType>(global_nodes_x, global_nodes_y, global_nodes_z,
ix+sx, iy+sy, iz+sz);
if (col_id >= 0 && col_id < global_nrows) {
GlobalOrdinalType col = mesh->map_id_to_row(col_id);
dest_cols[offset+nnz] = col;
dest_coefs[offset+nnz] = 0;
++nnz;
}
}
sort_if_needed(&dest_cols[offset], nnz);
}
};
#endif