/[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 6144 - (show annotations)
Wed Apr 6 05:25:13 2016 UTC (2 years, 7 months ago) by caltinay
File size: 16009 byte(s)
last round of namespacing defines.

1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2016 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 by Centre for Geoscience Computing (GeoComp)
14 *
15 *****************************************************************************/
16
17
18 /****************************************************************************
19
20 Finley: SystemMatrix
21
22 Adds the matrix array[Equa,Sol,NN,NN] to the matrix in.
23 The rows/columns are given by
24 i_Equa+Equa*Nodes_Equa[Nodes[j_Equa]] (i_Equa=0:Equa; j_Equa=0:NN_Equa).
25
26 The routine has to be called from a parallel region.
27 This routine assumes that in->Equa=in->Sol=1, i.e. array is fully packed.
28 TODO: the case in->Equa!=1
29 WARNING: MATRIX_FORMAT_CSC is not supported under MPI!
30
31 *****************************************************************************/
32
33 #include "Assemble.h"
34 #include <paso/SystemMatrix.h>
35
36 #ifdef ESYS_HAVE_TRILINOS
37 #include <trilinoswrap/TrilinosMatrixAdapter.h>
38
39 using esys_trilinos::TrilinosMatrixAdapter;
40 #endif
41
42 namespace finley {
43
44 void Assemble_addToSystemMatrix_CSC(paso::SystemMatrix* in, int NN_Equa,
45 const index_t* Nodes_Equa, int num_Equa,
46 int NN_Sol, const index_t* Nodes_Sol,
47 int num_Sol, const double* array);
48
49 #ifdef ESYS_HAVE_TRILINOS
50 void Assemble_addToSystemMatrix_Trilinos(TrilinosMatrixAdapter* in,
51 int NN_Equa, const index_t* Nodes_Equa, int num_Equa,
52 int NN_Sol, const index_t* Nodes_Sol, int num_Sol,
53 const double* array);
54 #endif
55
56 void Assemble_addToSystemMatrix_CSR(paso::SystemMatrix* in, int NN_Equa,
57 const index_t* Nodes_Equa, int num_Equa,
58 int NN_Sol, const index_t* Nodes_Sol,
59 int num_Sol, const double* array);
60
61 void Assemble_addToSystemMatrix(escript::ASM_ptr in, int NN_Equa,
62 const index_t* Nodes_Equa, int num_Equa,
63 int NN_Sol, const index_t* Nodes_Sol,
64 int num_Sol, const double* array)
65 {
66 using paso::SystemMatrix;
67 SystemMatrix* pmat(dynamic_cast<SystemMatrix*>(in.get()));
68
69 if (pmat) {
70 // call the right function depending on storage type
71 if (pmat->type & MATRIX_FORMAT_CSC) {
72 Assemble_addToSystemMatrix_CSC(pmat, NN_Equa, Nodes_Equa,
73 num_Equa, NN_Sol, Nodes_Sol,
74 num_Sol, array);
75 } else { // type == CSR
76 Assemble_addToSystemMatrix_CSR(pmat, NN_Equa, Nodes_Equa,
77 num_Equa, NN_Sol, Nodes_Sol,
78 num_Sol, array);
79 }
80 return;
81 }
82 #ifdef ESYS_HAVE_TRILINOS
83 TrilinosMatrixAdapter* tmat(dynamic_cast<TrilinosMatrixAdapter*>(in.get()));
84 if (tmat) {
85 Assemble_addToSystemMatrix_Trilinos(tmat, NN_Equa, Nodes_Equa,
86 num_Equa, NN_Sol, Nodes_Sol,
87 num_Sol, array);
88 return;
89 }
90 #endif
91 throw FinleyException("Assemble_addToSystemMatrix: unknown system matrix type.");
92 }
93
94 void Assemble_addToSystemMatrix_CSC(paso::SystemMatrix* in, int NN_Equa,
95 const index_t* Nodes_Equa, int num_Equa,
96 int NN_Sol, const index_t* Nodes_Sol,
97 int num_Sol, const double* array)
98 {
99 const int index_offset=(in->type & MATRIX_FORMAT_OFFSET1 ? 1:0);
100 const int row_block_size=in->row_block_size;
101 const int col_block_size=in->col_block_size;
102 const int block_size=in->block_size;
103 const int num_subblocks_Equa=num_Equa/row_block_size;
104 const int num_subblocks_Sol=num_Sol/col_block_size;
105 const dim_t numMyCols=in->pattern->mainPattern->numInput;
106 const dim_t numMyRows=in->pattern->mainPattern->numOutput;
107 const index_t *mainBlock_ptr=in->mainBlock->pattern->ptr;
108 const index_t *mainBlock_index=in->mainBlock->pattern->index;
109 double *mainBlock_val=in->mainBlock->val;
110 const index_t *col_coupleBlock_ptr=in->col_coupleBlock->pattern->ptr;
111 const index_t *col_coupleBlock_index=in->col_coupleBlock->pattern->index;
112 double *col_coupleBlock_val=in->col_coupleBlock->val;
113 //const index_t *row_coupleBlock_ptr=in->row_coupleBlock->pattern->ptr;
114 const index_t *row_coupleBlock_index=in->row_coupleBlock->pattern->index;
115 double *row_coupleBlock_val=in->row_coupleBlock->val;
116
117 for (int k_Sol=0; k_Sol<NN_Sol; ++k_Sol) {
118 // Down columns of array
119 const index_t j_Sol=Nodes_Sol[k_Sol];
120 for (int l_col=0; l_col<num_subblocks_Sol; ++l_col) {
121 const index_t i_col=j_Sol*num_subblocks_Sol+l_col;
122 if (i_col < numMyCols) {
123 for (int k_Equa=0;k_Equa<NN_Equa;++k_Equa) {
124 // Across cols of array
125 const index_t j_Equa=Nodes_Equa[k_Equa];
126 for (int l_row=0; l_row<num_subblocks_Equa; ++l_row) {
127 const index_t i_row=j_Equa*num_subblocks_Equa+index_offset+l_row;
128 if (i_row < numMyRows + index_offset ) {
129 for (index_t k=mainBlock_ptr[i_col]-index_offset; k<mainBlock_ptr[i_col+1]-index_offset; ++k) {
130 if (mainBlock_index[k]==i_row) {
131 // Entry array(k_Equa, j_Sol) is a block (col_block_size x col_block_size)
132 for (int ic=0; ic<col_block_size; ++ic) {
133 const int i_Sol=ic+col_block_size*l_col;
134 for (int ir=0; ir<row_block_size; ++ir) {
135 const index_t i_Equa=ir+row_block_size*l_row;
136 mainBlock_val[k*block_size+ir+row_block_size*ic]+=
137 array[INDEX4(i_Equa,i_Sol,k_Equa,k_Sol,num_Equa,num_Sol,NN_Equa)];
138 }
139 }
140 break;
141 }
142 }
143 } else {
144 for (index_t k=col_coupleBlock_ptr[i_col]-index_offset; k<col_coupleBlock_ptr[i_col+1]-index_offset; ++k) {
145 if (row_coupleBlock_index[k] == i_row-numMyRows) {
146 for (int ic=0; ic<col_block_size; ++ic) {
147 const int i_Sol=ic+col_block_size*l_col;
148 for (int ir=0; ir<row_block_size; ++ir) {
149 const index_t i_Equa=ir+row_block_size*l_row;
150 row_coupleBlock_val[k*block_size+ir+row_block_size*ic]+=
151 array[INDEX4(i_Equa,i_Sol,k_Equa,k_Sol,num_Equa,num_Sol,NN_Equa)];
152 }
153 }
154 break;
155 }
156 }
157 }
158 }
159 }
160 } else { // i_col >= numMyCols
161 for (int k_Equa=0;k_Equa<NN_Equa;++k_Equa) {
162 // Across rows of array
163 const index_t j_Equa=Nodes_Equa[k_Equa];
164 for (int l_row=0; l_row<num_subblocks_Equa; ++l_row) {
165 const index_t i_row=j_Equa*num_subblocks_Equa+index_offset+l_row;
166 if (i_row < numMyRows + index_offset ) {
167 for (index_t k=col_coupleBlock_ptr[i_col-numMyCols]-index_offset; k<col_coupleBlock_ptr[i_col-numMyCols+1]-index_offset; ++k) {
168 if (col_coupleBlock_index[k] == i_row) {
169 for (int ic=0; ic<col_block_size; ++ic) {
170 const int i_Sol=ic+col_block_size*l_col;
171 for (int ir=0; ir<row_block_size; ++ir) {
172 const int i_Equa=ir+row_block_size*l_row;
173 col_coupleBlock_val[k*block_size+ir+row_block_size*ic]+=
174 array[INDEX4(i_Equa,i_Sol,k_Equa,k_Sol,num_Equa,num_Sol,NN_Equa)];
175 }
176 }
177 break;
178 }
179 }
180 }
181 }
182 }
183 }
184 }
185 }
186 }
187
188 void Assemble_addToSystemMatrix_CSR(paso::SystemMatrix* in, int NN_Equa,
189 const index_t* Nodes_Equa, int num_Equa,
190 int NN_Sol, const index_t* Nodes_Sol,
191 int num_Sol, const double* array)
192 {
193 const int index_offset=(in->type & MATRIX_FORMAT_OFFSET1 ? 1:0);
194 const int row_block_size=in->row_block_size;
195 const int col_block_size=in->col_block_size;
196 const int block_size=in->block_size;
197 const int num_subblocks_Equa=num_Equa/row_block_size;
198 const int num_subblocks_Sol=num_Sol/col_block_size;
199 const dim_t numMyCols=in->pattern->mainPattern->numInput;
200 const dim_t numMyRows=in->pattern->mainPattern->numOutput;
201 const index_t *mainBlock_ptr=in->mainBlock->pattern->ptr;
202 const index_t *mainBlock_index=in->mainBlock->pattern->index;
203 double *mainBlock_val=in->mainBlock->val;
204 const index_t *col_coupleBlock_ptr=in->col_coupleBlock->pattern->ptr;
205 const index_t *col_coupleBlock_index=in->col_coupleBlock->pattern->index;
206 double *col_coupleBlock_val=in->col_coupleBlock->val;
207 const index_t *row_coupleBlock_ptr=in->row_coupleBlock->pattern->ptr;
208 const index_t *row_coupleBlock_index=in->row_coupleBlock->pattern->index;
209 double *row_coupleBlock_val=in->row_coupleBlock->val;
210
211 for (int k_Equa=0; k_Equa<NN_Equa; ++k_Equa) {
212 // Down columns of array
213 const index_t j_Equa=Nodes_Equa[k_Equa];
214 for (int l_row=0; l_row<num_subblocks_Equa; ++l_row) {
215 const index_t i_row=j_Equa*num_subblocks_Equa+l_row;
216 // only look at the matrix rows stored on this processor
217 if (i_row < numMyRows) {
218 for (int k_Sol=0; k_Sol<NN_Sol; ++k_Sol) {
219 // Across rows of array
220 const index_t j_Sol=Nodes_Sol[k_Sol];
221 for (int l_col=0; l_col<num_subblocks_Sol; ++l_col) {
222 // only look at the matrix rows stored on this processor
223 const index_t i_col=j_Sol*num_subblocks_Sol+index_offset+l_col;
224 if (i_col < numMyCols + index_offset ) {
225 for (index_t k=mainBlock_ptr[i_row]-index_offset; k<mainBlock_ptr[i_row+1]-index_offset; ++k) {
226 if (mainBlock_index[k]==i_col) {
227 // Entry array(k_Sol, j_Equa) is a block
228 // (row_block_size x col_block_size)
229 for (int ic=0; ic<col_block_size; ++ic) {
230 const int i_Sol=ic+col_block_size*l_col;
231 for (int ir=0; ir<row_block_size; ++ir) {
232 const index_t i_Equa=ir+row_block_size*l_row;
233 mainBlock_val[k*block_size+ir+row_block_size*ic]+=
234 array[INDEX4(i_Equa,i_Sol,k_Equa,k_Sol,num_Equa,num_Sol,NN_Equa)];
235 }
236 }
237 break;
238 }
239 }
240 } else {
241 for (index_t k=col_coupleBlock_ptr[i_row]-index_offset; k<col_coupleBlock_ptr[i_row+1]-index_offset; ++k) {
242 if (col_coupleBlock_index[k] == i_col-numMyCols) {
243 // Entry array(k_Sol, j_Equa) is a block
244 // (row_block_size x col_block_size)
245 for (int ic=0; ic<col_block_size; ++ic) {
246 const index_t i_Sol=ic+col_block_size*l_col;
247 for (int ir=0; ir<row_block_size; ++ir) {
248 const index_t i_Equa=ir+row_block_size*l_row;
249 col_coupleBlock_val[k*block_size+ir+row_block_size*ic]+=
250 array[INDEX4(i_Equa,i_Sol,k_Equa,k_Sol,num_Equa,num_Sol,NN_Equa)];
251 }
252 }
253 break;
254 }
255 }
256 }
257 }
258 }
259 } else { // i_row >= numMyRows
260 for (int k_Sol=0; k_Sol<NN_Sol; ++k_Sol) {
261 // Across rows of array
262 const index_t j_Sol=Nodes_Sol[k_Sol];
263 for (int l_col=0; l_col<num_subblocks_Sol; ++l_col) {
264 const index_t i_col=j_Sol*num_subblocks_Sol+index_offset+l_col;
265 if (i_col < numMyCols + index_offset ) {
266 for (index_t k=row_coupleBlock_ptr[i_row-numMyRows]-index_offset; k<row_coupleBlock_ptr[i_row-numMyRows+1]-index_offset; ++k) {
267 if (row_coupleBlock_index[k] == i_col) {
268 // Entry array(k_Sol, j_Equa) is a block
269 // (row_block_size x col_block_size)
270 for (int ic=0; ic<col_block_size; ++ic) {
271 const int i_Sol=ic+col_block_size*l_col;
272 for (int ir=0; ir<row_block_size; ++ir) {
273 const index_t i_Equa=ir+row_block_size*l_row;
274 row_coupleBlock_val[k*block_size+ir+row_block_size*ic]+=
275 array[INDEX4(i_Equa,i_Sol,k_Equa,k_Sol,num_Equa,num_Sol,NN_Equa)];
276 }
277 }
278 break;
279 }
280 }
281 }
282 }
283 }
284 }
285 }
286 }
287 }
288
289 #ifdef ESYS_HAVE_TRILINOS
290 void Assemble_addToSystemMatrix_Trilinos(TrilinosMatrixAdapter* in,
291 int NN_Equa, const index_t* Nodes_Equa,
292 int num_Equa, int NN_Sol,
293 const index_t* Nodes_Sol, int num_Sol,
294 const double* array)
295 {
296 std::vector<index_t> rowIdx(Nodes_Equa, Nodes_Equa+NN_Equa);
297 //std::vector<index_t> colIdx(Nodes_Sol, Nodes_Sol+NN_Sol);
298 std::vector<double> arr(array, array+(NN_Equa*NN_Sol*num_Sol*num_Equa));
299 in->add(rowIdx, arr);
300 }
301 #endif
302
303 } // namespace finley
304

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