/[escript]/branches/trilinos_from_5897/finley/src/Assemble_addToSystemMatrix.cpp
ViewVC logotype

Contents of /branches/trilinos_from_5897/finley/src/Assemble_addToSystemMatrix.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6058 - (show annotations)
Thu Mar 10 06:51:55 2016 UTC (3 years, 1 month ago) by caltinay
File size: 15968 byte(s)
added getPtr() to AbstractSystemMatrix so we can now use shared systemmatrix
pointers rather than circumventing them.

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 Open Software License version 3.0
9 * http://www.opensource.org/licenses/osl-3.0.php
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 USE_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 void Assemble_addToSystemMatrix_Trilinos(TrilinosMatrixAdapter* in,
50 int NN_Equa, const index_t* Nodes_Equa, int num_Equa,
51 int NN_Sol, const index_t* Nodes_Sol, int num_Sol,
52 const double* array);
53
54 void Assemble_addToSystemMatrix_CSR(paso::SystemMatrix* in, int NN_Equa,
55 const index_t* Nodes_Equa, int num_Equa,
56 int NN_Sol, const index_t* Nodes_Sol,
57 int num_Sol, const double* array);
58
59 void Assemble_addToSystemMatrix(escript::ASM_ptr in, int NN_Equa,
60 const index_t* Nodes_Equa, int num_Equa,
61 int NN_Sol, const index_t* Nodes_Sol,
62 int num_Sol, const double* array)
63 {
64 using paso::SystemMatrix;
65 SystemMatrix* pmat(dynamic_cast<SystemMatrix*>(in.get()));
66
67 if (pmat) {
68 // call the right function depending on storage type
69 if (pmat->type & MATRIX_FORMAT_CSC) {
70 Assemble_addToSystemMatrix_CSC(pmat, NN_Equa, Nodes_Equa,
71 num_Equa, NN_Sol, Nodes_Sol,
72 num_Sol, array);
73 } else { // type == CSR
74 Assemble_addToSystemMatrix_CSR(pmat, NN_Equa, Nodes_Equa,
75 num_Equa, NN_Sol, Nodes_Sol,
76 num_Sol, array);
77 }
78 return;
79 }
80 #ifdef USE_TRILINOS
81 TrilinosMatrixAdapter* tmat(dynamic_cast<TrilinosMatrixAdapter*>(in.get()));
82 if (tmat) {
83 Assemble_addToSystemMatrix_Trilinos(tmat, NN_Equa, Nodes_Equa,
84 num_Equa, NN_Sol, Nodes_Sol,
85 num_Sol, array);
86 return;
87 }
88 #endif
89 throw FinleyException("Assemble_addToSystemMatrix: unknown system matrix type.");
90 }
91
92 void Assemble_addToSystemMatrix_CSC(paso::SystemMatrix* in, int NN_Equa,
93 const index_t* Nodes_Equa, int num_Equa,
94 int NN_Sol, const index_t* Nodes_Sol,
95 int num_Sol, const double* array)
96 {
97 const int index_offset=(in->type & MATRIX_FORMAT_OFFSET1 ? 1:0);
98 const int row_block_size=in->row_block_size;
99 const int col_block_size=in->col_block_size;
100 const int block_size=in->block_size;
101 const int num_subblocks_Equa=num_Equa/row_block_size;
102 const int num_subblocks_Sol=num_Sol/col_block_size;
103 const dim_t numMyCols=in->pattern->mainPattern->numInput;
104 const dim_t numMyRows=in->pattern->mainPattern->numOutput;
105 const index_t *mainBlock_ptr=in->mainBlock->pattern->ptr;
106 const index_t *mainBlock_index=in->mainBlock->pattern->index;
107 double *mainBlock_val=in->mainBlock->val;
108 const index_t *col_coupleBlock_ptr=in->col_coupleBlock->pattern->ptr;
109 const index_t *col_coupleBlock_index=in->col_coupleBlock->pattern->index;
110 double *col_coupleBlock_val=in->col_coupleBlock->val;
111 //const index_t *row_coupleBlock_ptr=in->row_coupleBlock->pattern->ptr;
112 const index_t *row_coupleBlock_index=in->row_coupleBlock->pattern->index;
113 double *row_coupleBlock_val=in->row_coupleBlock->val;
114
115 for (int k_Sol=0; k_Sol<NN_Sol; ++k_Sol) {
116 // Down columns of array
117 const index_t j_Sol=Nodes_Sol[k_Sol];
118 for (int l_col=0; l_col<num_subblocks_Sol; ++l_col) {
119 const index_t i_col=j_Sol*num_subblocks_Sol+l_col;
120 if (i_col < numMyCols) {
121 for (int k_Equa=0;k_Equa<NN_Equa;++k_Equa) {
122 // Across cols of array
123 const index_t j_Equa=Nodes_Equa[k_Equa];
124 for (int l_row=0; l_row<num_subblocks_Equa; ++l_row) {
125 const index_t i_row=j_Equa*num_subblocks_Equa+index_offset+l_row;
126 if (i_row < numMyRows + index_offset ) {
127 for (index_t k=mainBlock_ptr[i_col]-index_offset; k<mainBlock_ptr[i_col+1]-index_offset; ++k) {
128 if (mainBlock_index[k]==i_row) {
129 // Entry array(k_Equa, j_Sol) is a block (col_block_size x col_block_size)
130 for (int ic=0; ic<col_block_size; ++ic) {
131 const int i_Sol=ic+col_block_size*l_col;
132 for (int ir=0; ir<row_block_size; ++ir) {
133 const index_t i_Equa=ir+row_block_size*l_row;
134 mainBlock_val[k*block_size+ir+row_block_size*ic]+=
135 array[INDEX4(i_Equa,i_Sol,k_Equa,k_Sol,num_Equa,num_Sol,NN_Equa)];
136 }
137 }
138 break;
139 }
140 }
141 } else {
142 for (index_t k=col_coupleBlock_ptr[i_col]-index_offset; k<col_coupleBlock_ptr[i_col+1]-index_offset; ++k) {
143 if (row_coupleBlock_index[k] == i_row-numMyRows) {
144 for (int ic=0; ic<col_block_size; ++ic) {
145 const int i_Sol=ic+col_block_size*l_col;
146 for (int ir=0; ir<row_block_size; ++ir) {
147 const index_t i_Equa=ir+row_block_size*l_row;
148 row_coupleBlock_val[k*block_size+ir+row_block_size*ic]+=
149 array[INDEX4(i_Equa,i_Sol,k_Equa,k_Sol,num_Equa,num_Sol,NN_Equa)];
150 }
151 }
152 break;
153 }
154 }
155 }
156 }
157 }
158 } else { // i_col >= numMyCols
159 for (int k_Equa=0;k_Equa<NN_Equa;++k_Equa) {
160 // Across rows of array
161 const index_t j_Equa=Nodes_Equa[k_Equa];
162 for (int l_row=0; l_row<num_subblocks_Equa; ++l_row) {
163 const index_t i_row=j_Equa*num_subblocks_Equa+index_offset+l_row;
164 if (i_row < numMyRows + index_offset ) {
165 for (index_t k=col_coupleBlock_ptr[i_col-numMyCols]-index_offset; k<col_coupleBlock_ptr[i_col-numMyCols+1]-index_offset; ++k) {
166 if (col_coupleBlock_index[k] == i_row) {
167 for (int ic=0; ic<col_block_size; ++ic) {
168 const int i_Sol=ic+col_block_size*l_col;
169 for (int ir=0; ir<row_block_size; ++ir) {
170 const int i_Equa=ir+row_block_size*l_row;
171 col_coupleBlock_val[k*block_size+ir+row_block_size*ic]+=
172 array[INDEX4(i_Equa,i_Sol,k_Equa,k_Sol,num_Equa,num_Sol,NN_Equa)];
173 }
174 }
175 break;
176 }
177 }
178 }
179 }
180 }
181 }
182 }
183 }
184 }
185
186 void Assemble_addToSystemMatrix_CSR(paso::SystemMatrix* in, int NN_Equa,
187 const index_t* Nodes_Equa, int num_Equa,
188 int NN_Sol, const index_t* Nodes_Sol,
189 int num_Sol, const double* array)
190 {
191 const int index_offset=(in->type & MATRIX_FORMAT_OFFSET1 ? 1:0);
192 const int row_block_size=in->row_block_size;
193 const int col_block_size=in->col_block_size;
194 const int block_size=in->block_size;
195 const int num_subblocks_Equa=num_Equa/row_block_size;
196 const int num_subblocks_Sol=num_Sol/col_block_size;
197 const dim_t numMyCols=in->pattern->mainPattern->numInput;
198 const dim_t numMyRows=in->pattern->mainPattern->numOutput;
199 const index_t *mainBlock_ptr=in->mainBlock->pattern->ptr;
200 const index_t *mainBlock_index=in->mainBlock->pattern->index;
201 double *mainBlock_val=in->mainBlock->val;
202 const index_t *col_coupleBlock_ptr=in->col_coupleBlock->pattern->ptr;
203 const index_t *col_coupleBlock_index=in->col_coupleBlock->pattern->index;
204 double *col_coupleBlock_val=in->col_coupleBlock->val;
205 const index_t *row_coupleBlock_ptr=in->row_coupleBlock->pattern->ptr;
206 const index_t *row_coupleBlock_index=in->row_coupleBlock->pattern->index;
207 double *row_coupleBlock_val=in->row_coupleBlock->val;
208
209 for (int k_Equa=0; k_Equa<NN_Equa; ++k_Equa) {
210 // Down columns of array
211 const index_t j_Equa=Nodes_Equa[k_Equa];
212 for (int l_row=0; l_row<num_subblocks_Equa; ++l_row) {
213 const index_t i_row=j_Equa*num_subblocks_Equa+l_row;
214 // only look at the matrix rows stored on this processor
215 if (i_row < numMyRows) {
216 for (int k_Sol=0; k_Sol<NN_Sol; ++k_Sol) {
217 // Across rows of array
218 const index_t j_Sol=Nodes_Sol[k_Sol];
219 for (int l_col=0; l_col<num_subblocks_Sol; ++l_col) {
220 // only look at the matrix rows stored on this processor
221 const index_t i_col=j_Sol*num_subblocks_Sol+index_offset+l_col;
222 if (i_col < numMyCols + index_offset ) {
223 for (index_t k=mainBlock_ptr[i_row]-index_offset; k<mainBlock_ptr[i_row+1]-index_offset; ++k) {
224 if (mainBlock_index[k]==i_col) {
225 // Entry array(k_Sol, j_Equa) is a block
226 // (row_block_size x col_block_size)
227 for (int ic=0; ic<col_block_size; ++ic) {
228 const int i_Sol=ic+col_block_size*l_col;
229 for (int ir=0; ir<row_block_size; ++ir) {
230 const index_t i_Equa=ir+row_block_size*l_row;
231 mainBlock_val[k*block_size+ir+row_block_size*ic]+=
232 array[INDEX4(i_Equa,i_Sol,k_Equa,k_Sol,num_Equa,num_Sol,NN_Equa)];
233 }
234 }
235 break;
236 }
237 }
238 } else {
239 for (index_t k=col_coupleBlock_ptr[i_row]-index_offset; k<col_coupleBlock_ptr[i_row+1]-index_offset; ++k) {
240 if (col_coupleBlock_index[k] == i_col-numMyCols) {
241 // Entry array(k_Sol, j_Equa) is a block
242 // (row_block_size x col_block_size)
243 for (int ic=0; ic<col_block_size; ++ic) {
244 const index_t i_Sol=ic+col_block_size*l_col;
245 for (int ir=0; ir<row_block_size; ++ir) {
246 const index_t i_Equa=ir+row_block_size*l_row;
247 col_coupleBlock_val[k*block_size+ir+row_block_size*ic]+=
248 array[INDEX4(i_Equa,i_Sol,k_Equa,k_Sol,num_Equa,num_Sol,NN_Equa)];
249 }
250 }
251 break;
252 }
253 }
254 }
255 }
256 }
257 } else { // i_row >= numMyRows
258 for (int k_Sol=0; k_Sol<NN_Sol; ++k_Sol) {
259 // Across rows of array
260 const index_t j_Sol=Nodes_Sol[k_Sol];
261 for (int l_col=0; l_col<num_subblocks_Sol; ++l_col) {
262 const index_t i_col=j_Sol*num_subblocks_Sol+index_offset+l_col;
263 if (i_col < numMyCols + index_offset ) {
264 for (index_t k=row_coupleBlock_ptr[i_row-numMyRows]-index_offset; k<row_coupleBlock_ptr[i_row-numMyRows+1]-index_offset; ++k) {
265 if (row_coupleBlock_index[k] == i_col) {
266 // Entry array(k_Sol, j_Equa) is a block
267 // (row_block_size x col_block_size)
268 for (int ic=0; ic<col_block_size; ++ic) {
269 const int i_Sol=ic+col_block_size*l_col;
270 for (int ir=0; ir<row_block_size; ++ir) {
271 const index_t i_Equa=ir+row_block_size*l_row;
272 row_coupleBlock_val[k*block_size+ir+row_block_size*ic]+=
273 array[INDEX4(i_Equa,i_Sol,k_Equa,k_Sol,num_Equa,num_Sol,NN_Equa)];
274 }
275 }
276 break;
277 }
278 }
279 }
280 }
281 }
282 }
283 }
284 }
285 }
286
287 #ifdef USE_TRILINOS
288 void Assemble_addToSystemMatrix_Trilinos(TrilinosMatrixAdapter* in,
289 int NN_Equa, const index_t* Nodes_Equa,
290 int num_Equa, int NN_Sol,
291 const index_t* Nodes_Sol, int num_Sol,
292 const double* array)
293 {
294 std::vector<index_t> rowIdx(Nodes_Equa, Nodes_Equa+NN_Equa);
295 //std::vector<index_t> colIdx(Nodes_Sol, Nodes_Sol+NN_Sol);
296 std::vector<double> arr(array, array+(NN_Equa*NN_Sol*num_Sol*num_Equa));
297 in->add(rowIdx, arr);
298 }
299 #endif
300
301 } // namespace finley
302

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 /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,5898-6007

  ViewVC Help
Powered by ViewVC 1.1.26