/[escript]/trunk/finley/src/Assemble_addToSystemMatrix.cpp
ViewVC logotype

Contents of /trunk/finley/src/Assemble_addToSystemMatrix.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6939 - (show annotations)
Mon Jan 20 03:37:18 2020 UTC (4 months, 2 weeks ago) by uqaeller
File size: 16500 byte(s)
Updated the copyright header.


1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2020 by The University of Queensland
5 * http://www.uq.edu.au
6 *
7 * Primary Business: Queensland, Australia
8 * Licensed under the Apache License, version 2.0
9 * http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 * Development 2012-2013 by School of Earth Sciences
13 * Development from 2014-2017 by Centre for Geoscience Computing (GeoComp)
14 * Development from 2019 by School of Earth and Environmental Sciences
15 **
16 *****************************************************************************/
17
18 #include "Assemble.h"
19
20 #ifdef ESYS_HAVE_PASO
21 #include <paso/SystemMatrix.h>
22 #endif
23
24 #ifdef ESYS_HAVE_TRILINOS
25 #include <trilinoswrap/TrilinosMatrixAdapter.h>
26
27 using esys_trilinos::TrilinosMatrixAdapter;
28 #endif
29
30 namespace finley {
31
32 using escript::DataTypes::real_t;
33 using escript::DataTypes::cplx_t;
34
35 #ifdef ESYS_HAVE_PASO
36 static void addToSystemMatrixPasoCSC(paso::SystemMatrix* S, int NN_Equa,
37 const index_t* Nodes_Equa, int num_Equa,
38 int NN_Sol, const index_t* Nodes_Sol,
39 int num_Sol, const real_t* array);
40
41 static void addToSystemMatrixPasoCSR(paso::SystemMatrix* S, int NN_Equa,
42 const index_t* Nodes_Equa, int num_Equa,
43 int NN_Sol, const index_t* Nodes_Sol,
44 int num_Sol, const real_t* array);
45 #endif
46
47 template<>
48 void Assemble_addToSystemMatrix<real_t>(escript::ASM_ptr S, int NN_Equa,
49 const index_t* Nodes_Equa, int num_Equa,
50 int NN_Sol, const index_t* Nodes_Sol,
51 int num_Sol, const real_t* array)
52 {
53 #ifdef ESYS_HAVE_PASO
54 paso::SystemMatrix* pmat = dynamic_cast<paso::SystemMatrix*>(S.get());
55 if (pmat) {
56 // call the right function depending on storage type
57 if (pmat->type & MATRIX_FORMAT_CSC) {
58 addToSystemMatrixPasoCSC(pmat, NN_Equa, Nodes_Equa,
59 num_Equa, NN_Sol, Nodes_Sol,
60 num_Sol, array);
61 } else { // type == CSR
62 addToSystemMatrixPasoCSR(pmat, NN_Equa, Nodes_Equa,
63 num_Equa, NN_Sol, Nodes_Sol,
64 num_Sol, array);
65 }
66 return;
67 }
68 #endif
69 #ifdef ESYS_HAVE_TRILINOS
70 TrilinosMatrixAdapter* tmat(dynamic_cast<TrilinosMatrixAdapter*>(S.get()));
71 if (tmat) {
72 IndexVector rowIdx(Nodes_Equa, Nodes_Equa+NN_Equa);
73 //IndexVector colIdx(Nodes_Sol, Nodes_Sol+NN_Sol);
74 std::vector<real_t> arr(array, array+(NN_Equa*NN_Sol*num_Sol*num_Equa));
75 tmat->add(rowIdx, arr);
76 return;
77 }
78 #endif
79 throw FinleyException("Assemble_addToSystemMatrix: unknown system "
80 "matrix type.");
81 }
82
83 template<>
84 void Assemble_addToSystemMatrix<cplx_t>(escript::ASM_ptr S, int NN_Equa,
85 const index_t* Nodes_Equa, int num_Equa,
86 int NN_Sol, const index_t* Nodes_Sol,
87 int num_Sol, const cplx_t* array)
88 {
89 #ifdef ESYS_HAVE_TRILINOS
90 TrilinosMatrixAdapter* tmat = dynamic_cast<TrilinosMatrixAdapter*>(S.get());
91 if (tmat) {
92 IndexVector rowIdx(Nodes_Equa, Nodes_Equa+NN_Equa);
93 //IndexVector colIdx(Nodes_Sol, Nodes_Sol+NN_Sol);
94 std::vector<cplx_t> arr(array, array+(NN_Equa*NN_Sol*num_Sol*num_Equa));
95 tmat->add(rowIdx, arr);
96 return;
97 }
98 #endif
99 throw FinleyException("addToSystemMatrix: only Trilinos matrices support "
100 "complex-valued assembly!");
101 }
102
103 #ifdef ESYS_HAVE_PASO
104 void addToSystemMatrixPasoCSC(paso::SystemMatrix* in, int NN_Equa,
105 const index_t* Nodes_Equa, int num_Equa,
106 int NN_Sol, const index_t* Nodes_Sol,
107 int num_Sol, const real_t* array)
108 {
109 const int index_offset = (in->type & MATRIX_FORMAT_OFFSET1 ? 1 : 0);
110 const int row_block_size = in->row_block_size;
111 const int col_block_size = in->col_block_size;
112 const int block_size = in->block_size;
113 const int num_subblocks_Equa = num_Equa/row_block_size;
114 const int num_subblocks_Sol = num_Sol/col_block_size;
115 const dim_t numMyCols = in->pattern->mainPattern->numInput;
116 const dim_t numMyRows = in->pattern->mainPattern->numOutput;
117
118 const index_t* mainBlock_ptr = in->mainBlock->pattern->ptr;
119 const index_t* mainBlock_index = in->mainBlock->pattern->index;
120 real_t* mainBlock_val = in->mainBlock->val;
121 const index_t* col_coupleBlock_ptr = in->col_coupleBlock->pattern->ptr;
122 const index_t* col_coupleBlock_index = in->col_coupleBlock->pattern->index;
123 real_t* col_coupleBlock_val = in->col_coupleBlock->val;
124 //const index_t* row_coupleBlock_ptr = in->row_coupleBlock->pattern->ptr;
125 const index_t* row_coupleBlock_index = in->row_coupleBlock->pattern->index;
126 real_t* row_coupleBlock_val = in->row_coupleBlock->val;
127
128 for (int k_Sol = 0; k_Sol < NN_Sol; ++k_Sol) {
129 // Down columns of array
130 const index_t j_Sol = Nodes_Sol[k_Sol];
131 for (int l_col = 0; l_col < num_subblocks_Sol; ++l_col) {
132 const index_t i_col = j_Sol * num_subblocks_Sol + l_col;
133 if (i_col < numMyCols) {
134 for (int k_Equa = 0; k_Equa < NN_Equa; ++k_Equa) {
135 // Across cols of array
136 const index_t j_Equa = Nodes_Equa[k_Equa];
137 for (int l_row = 0; l_row < num_subblocks_Equa; ++l_row) {
138 const index_t i_row = j_Equa*num_subblocks_Equa+index_offset+l_row;
139 if (i_row < numMyRows + index_offset ) {
140 for (index_t k = mainBlock_ptr[i_col]-index_offset;
141 k < mainBlock_ptr[i_col + 1]-index_offset; ++k) {
142 if (mainBlock_index[k] == i_row) {
143 // Entry array(k_Equa, j_Sol) is a block
144 // (col_block_size x col_block_size)
145 for (int ic = 0; ic < col_block_size; ++ic) {
146 const int i_Sol = ic + col_block_size * l_col;
147 for (int ir = 0; ir < row_block_size; ++ir) {
148 const int i_Eq = ir + row_block_size * l_row;
149 mainBlock_val[k*block_size + ir + row_block_size*ic] +=
150 array[INDEX4
151 (i_Eq, i_Sol, k_Equa, k_Sol, num_Equa, num_Sol, NN_Equa)];
152 }
153 }
154 break;
155 }
156 }
157 } else {
158 for (index_t k = col_coupleBlock_ptr[i_col]-index_offset;
159 k < col_coupleBlock_ptr[i_col + 1]-index_offset; ++k) {
160 if (row_coupleBlock_index[k] == i_row - numMyRows) {
161 for (int ic = 0; ic < col_block_size; ++ic) {
162 const int i_Sol = ic + col_block_size * l_col;
163 for (int ir = 0; ir < row_block_size; ++ir) {
164 const int i_Eq = ir + row_block_size * l_row;
165 row_coupleBlock_val[k*block_size + ir + row_block_size*ic] +=
166 array[INDEX4
167 (i_Eq, i_Sol, k_Equa, k_Sol, num_Equa, num_Sol, NN_Equa)];
168 }
169 }
170 break;
171 }
172 }
173 }
174 }
175 }
176 } else { // i_col >= numMyCols
177 for (int k_Equa = 0; k_Equa < NN_Equa; ++k_Equa) {
178 // Across rows of array
179 const index_t j_Equa = Nodes_Equa[k_Equa];
180 for (int l_row = 0; l_row < num_subblocks_Equa; ++l_row) {
181 const index_t i_row = j_Equa * num_subblocks_Equa + index_offset + l_row;
182 if (i_row < numMyRows + index_offset) {
183 for (index_t k = col_coupleBlock_ptr[i_col-numMyCols]-index_offset;
184 k < col_coupleBlock_ptr[i_col - numMyCols + 1] - index_offset; ++k) {
185 if (col_coupleBlock_index[k] == i_row) {
186 for (int ic = 0; ic < col_block_size; ++ic) {
187 const int i_Sol = ic + col_block_size * l_col;
188 for (int ir = 0; ir < row_block_size; ++ir) {
189 const int i_Eq = ir + row_block_size * l_row;
190 col_coupleBlock_val[k*block_size + ir + row_block_size*ic] +=
191 array[INDEX4
192 (i_Eq, i_Sol, k_Equa, k_Sol, num_Equa, num_Sol, NN_Equa)];
193 }
194 }
195 break;
196 }
197 }
198 }
199 }
200 }
201 }
202 }
203 }
204 }
205
206 void addToSystemMatrixPasoCSR(paso::SystemMatrix* in, int NN_Equa,
207 const index_t* Nodes_Equa, int num_Equa,
208 int NN_Sol, const index_t* Nodes_Sol,
209 int num_Sol, const real_t* array)
210 {
211 const int index_offset = (in->type & MATRIX_FORMAT_OFFSET1 ? 1 : 0);
212 const int row_block_size = in->row_block_size;
213 const int col_block_size = in->col_block_size;
214 const int block_size = in->block_size;
215 const int num_subblocks_Equa = num_Equa / row_block_size;
216 const int num_subblocks_Sol = num_Sol / col_block_size;
217 const dim_t numMyCols = in->pattern->mainPattern->numInput;
218 const dim_t numMyRows = in->pattern->mainPattern->numOutput;
219
220 const index_t* mainBlock_ptr = in->mainBlock->pattern->ptr;
221 const index_t* mainBlock_index = in->mainBlock->pattern->index;
222 real_t* mainBlock_val = in->mainBlock->val;
223 const index_t* col_coupleBlock_ptr = in->col_coupleBlock->pattern->ptr;
224 const index_t* col_coupleBlock_index = in->col_coupleBlock->pattern->index;
225 real_t* col_coupleBlock_val = in->col_coupleBlock->val;
226 const index_t* row_coupleBlock_ptr = in->row_coupleBlock->pattern->ptr;
227 const index_t* row_coupleBlock_index = in->row_coupleBlock->pattern->index;
228 real_t* row_coupleBlock_val = in->row_coupleBlock->val;
229
230 for (int k_Equa = 0; k_Equa < NN_Equa; ++k_Equa) {
231 // Down columns of array
232 const index_t j_Equa = Nodes_Equa[k_Equa];
233 for (int l_row = 0; l_row<num_subblocks_Equa; ++l_row) {
234 const index_t i_row = j_Equa*num_subblocks_Equa+l_row;
235 // only look at the matrix rows stored on this processor
236 if (i_row < numMyRows) {
237 for (int k_Sol=0; k_Sol<NN_Sol; ++k_Sol) {
238 // Across rows of array
239 const index_t j_Sol=Nodes_Sol[k_Sol];
240 for (int l_col=0; l_col<num_subblocks_Sol; ++l_col) {
241 // only look at the matrix rows stored on this processor
242 const index_t i_col = j_Sol * num_subblocks_Sol + index_offset + l_col;
243 if (i_col < numMyCols + index_offset) {
244 for (index_t k = mainBlock_ptr[i_row] - index_offset;
245 k < mainBlock_ptr[i_row + 1] - index_offset; ++k) {
246 if (mainBlock_index[k] == i_col) {
247 // Entry array(k_Sol, j_Equa) is a block
248 // (row_block_size x col_block_size)
249 for (int ic = 0; ic < col_block_size; ++ic) {
250 const int i_Sol = ic + col_block_size * l_col;
251 for (int ir = 0; ir < row_block_size; ++ir) {
252 const int i_Eq = ir + row_block_size * l_row;
253 mainBlock_val[k*block_size + ir + row_block_size*ic]+=
254 array[INDEX4
255 (i_Eq, i_Sol, k_Equa, k_Sol, num_Equa, num_Sol, NN_Equa)];
256 }
257 }
258 break;
259 }
260 }
261 } else {
262 for (index_t k = col_coupleBlock_ptr[i_row] - index_offset;
263 k < col_coupleBlock_ptr[i_row + 1] - index_offset; ++k) {
264 if (col_coupleBlock_index[k] == i_col - numMyCols) {
265 // Entry array(k_Sol, j_Equa) is a block
266 // (row_block_size x col_block_size)
267 for (int ic = 0; ic < col_block_size; ++ic) {
268 const int i_Sol = ic + col_block_size * l_col;
269 for (int ir = 0; ir < row_block_size; ++ir) {
270 const int i_Eq = ir+row_block_size*l_row;
271 col_coupleBlock_val[k*block_size + ir + row_block_size*ic]+=
272 array[INDEX4
273 (i_Eq, i_Sol, k_Equa, k_Sol, num_Equa, num_Sol, NN_Equa)];
274 }
275 }
276 break;
277 }
278 }
279 }
280 }
281 }
282 } else { // i_row >= numMyRows
283 for (int k_Sol = 0; k_Sol < NN_Sol; ++k_Sol) {
284 // Across rows of array
285 const index_t j_Sol = Nodes_Sol[k_Sol];
286 for (int l_col = 0; l_col < num_subblocks_Sol; ++l_col) {
287 const index_t i_col = j_Sol * num_subblocks_Sol + index_offset + l_col;
288 if (i_col < numMyCols + index_offset) {
289 for (index_t k = row_coupleBlock_ptr[i_row - numMyRows] - index_offset;
290 k < row_coupleBlock_ptr[i_row - numMyRows + 1] - index_offset; ++k) {
291 if (row_coupleBlock_index[k] == i_col) {
292 // Entry array(k_Sol, j_Equa) is a block
293 // (row_block_size x col_block_size)
294 for (int ic = 0; ic < col_block_size; ++ic) {
295 const int i_Sol = ic + col_block_size * l_col;
296 for (int ir = 0; ir < row_block_size; ++ir) {
297 const int i_Eq = ir + row_block_size * l_row;
298 row_coupleBlock_val[k*block_size + ir + row_block_size*ic]+=
299 array[INDEX4
300 (i_Eq, i_Sol, k_Equa, k_Sol, num_Equa, num_Sol, NN_Equa)];
301 }
302 }
303 break;
304 }
305 }
306 }
307 }
308 }
309 }
310 }
311 }
312 }
313 #endif // ESYS_HAVE_PASO
314
315 } // namespace finley
316

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision
svn:mergeinfo /branches/4.0fordebian/finley/src/Assemble_addToSystemMatrix.cpp:5567-5588 /branches/lapack2681/finley/src/Assemble_addToSystemMatrix.cpp:2682-2741 /branches/pasowrap/finley/src/Assemble_addToSystemMatrix.cpp:3661-3674 /branches/py3_attempt2/finley/src/Assemble_addToSystemMatrix.cpp:3871-3891 /branches/restext/finley/src/Assemble_addToSystemMatrix.cpp:2610-2624 /branches/ripleygmg_from_3668/finley/src/Assemble_addToSystemMatrix.cpp:3669-3791 /branches/stage3.0/finley/src/Assemble_addToSystemMatrix.cpp:2569-2590 /branches/symbolic_from_3470/finley/src/Assemble_addToSystemMatrix.cpp:3471-3974 /branches/trilinos_from_5897/finley/src/Assemble_addToSystemMatrix.cpp:5898-6118 /release/3.0/finley/src/Assemble_addToSystemMatrix.cpp:2591-2601 /release/4.0/finley/src/Assemble_addToSystemMatrix.cpp:5380-5406 /trunk/finley/src/Assemble_addToSystemMatrix.cpp:4257-4344

  ViewVC Help
Powered by ViewVC 1.1.26