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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6079 - (show annotations)
Mon Mar 21 12:22:38 2016 UTC (2 years, 11 months ago) by caltinay
File size: 11323 byte(s)
Big commit - making dudley much more like finley to make it more
managable. Fixed quite a few issues that had been fixed in finley.
Disposed of all ReducedNode/ReducedDOF entities that dudley never supported.
Compiles and passes tests.

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 #include "Assemble.h"
18
19 #include <paso/SystemMatrix.h>
20
21 namespace dudley {
22
23 void Assemble_addToSystemMatrix(escript::ASM_ptr S, dim_t NN_Equa,
24 const index_t* Nodes_Equa, dim_t num_Equa,
25 dim_t NN_Sol, const index_t* Nodes_Sol,
26 dim_t num_Sol, const double *array)
27 {
28 paso::SystemMatrix* in = dynamic_cast<paso::SystemMatrix*>(S.get());
29 if (!in) {
30 throw DudleyException(
31 "Assemble_addToSystemMatrix: unsupported matrix type.");
32 }
33 index_t index_offset = (in->type & MATRIX_FORMAT_OFFSET1 ? 1 : 0);
34 dim_t k_Equa, j_Equa, j_Sol, k_Sol, i_Equa, i_Sol, l_col, l_row, ic, ir, k, i_row, i_col;
35 index_t *mainBlock_ptr, *mainBlock_index, *col_coupleBlock_ptr, *col_coupleBlock_index, *row_coupleBlock_ptr,
36 *row_coupleBlock_index;
37 double *mainBlock_val, *row_coupleBlock_val, *col_coupleBlock_val;
38 dim_t row_block_size = in->row_block_size;
39 dim_t col_block_size = in->col_block_size;
40 dim_t block_size = in->block_size;
41 dim_t num_subblocks_Equa = num_Equa / row_block_size;
42 dim_t num_subblocks_Sol = num_Sol / col_block_size;
43 dim_t numMyCols = in->pattern->mainPattern->numInput;
44 dim_t numMyRows = in->pattern->mainPattern->numOutput;
45
46 if (in->type & MATRIX_FORMAT_CSC)
47 {
48 /* MATRIX_FORMAT_CSC does not support MPI !!!!! */
49 mainBlock_ptr = in->mainBlock->pattern->ptr;
50 mainBlock_index = in->mainBlock->pattern->index;
51 mainBlock_val = in->mainBlock->val;
52 col_coupleBlock_ptr = in->col_coupleBlock->pattern->ptr;
53 col_coupleBlock_index = in->col_coupleBlock->pattern->index;
54 col_coupleBlock_val = in->col_coupleBlock->val;
55 row_coupleBlock_ptr = in->row_coupleBlock->pattern->ptr;
56 row_coupleBlock_index = in->row_coupleBlock->pattern->index;
57 row_coupleBlock_val = in->row_coupleBlock->val;
58
59 for (k_Sol = 0; k_Sol < NN_Sol; ++k_Sol)
60 { /* Down columns of array */
61 j_Sol = Nodes_Sol[k_Sol];
62 for (l_col = 0; l_col < num_subblocks_Sol; ++l_col)
63 {
64 i_col = j_Sol * num_subblocks_Sol + l_col;
65 if (i_col < numMyCols)
66 {
67 for (k_Equa = 0; k_Equa < NN_Equa; ++k_Equa)
68 { /* Across cols of array */
69 j_Equa = Nodes_Equa[k_Equa];
70 for (l_row = 0; l_row < num_subblocks_Equa; ++l_row)
71 {
72 i_row = j_Equa * num_subblocks_Equa + index_offset + l_row;
73 if (i_row < numMyRows + index_offset)
74 {
75 for (k = mainBlock_ptr[i_col] - index_offset;
76 k < mainBlock_ptr[i_col + 1] - index_offset; ++k)
77 {
78 if (mainBlock_index[k] == i_row)
79 {
80 /* Entry array(k_Equa, j_Sol) is a block (col_block_size x col_block_size) */
81 for (ic = 0; ic < col_block_size; ++ic)
82 {
83 i_Sol = ic + col_block_size * l_col;;
84 for (ir = 0; ir < row_block_size; ++ir)
85 {
86 i_Equa = ir + row_block_size * l_row;
87 mainBlock_val[k * block_size + ir + row_block_size * ic] +=
88 array[INDEX4
89 (i_Equa, i_Sol, k_Equa, k_Sol, num_Equa, num_Sol, NN_Equa)];
90 }
91 }
92 break;
93 }
94 }
95 }
96 else
97 {
98 for (k = col_coupleBlock_ptr[i_col] - index_offset;
99 k < col_coupleBlock_ptr[i_col + 1] - index_offset; ++k)
100 {
101 if (row_coupleBlock_index[k] == i_row - numMyRows)
102 {
103 for (ic = 0; ic < col_block_size; ++ic)
104 {
105 i_Sol = ic + col_block_size * l_col;
106 for (ir = 0; ir < row_block_size; ++ir)
107 {
108 i_Equa = ir + row_block_size * l_row;
109 row_coupleBlock_val[k * block_size + ir + row_block_size * ic] +=
110 array[INDEX4
111 (i_Equa, i_Sol, k_Equa, k_Sol, num_Equa, num_Sol, NN_Equa)];;
112 }
113 }
114 break;
115 }
116 }
117 }
118 }
119 }
120 }
121 else
122 {
123 for (k_Equa = 0; k_Equa < NN_Equa; ++k_Equa)
124 { /* Across rows of array */
125 j_Equa = Nodes_Equa[k_Equa];
126 for (l_row = 0; l_row < num_subblocks_Equa; ++l_row)
127 {
128 i_row = j_Equa * num_subblocks_Equa + index_offset + l_row;
129 if (i_row < numMyRows + index_offset)
130 {
131 for (k = col_coupleBlock_ptr[i_col - numMyCols] - index_offset;
132 k < col_coupleBlock_ptr[i_col - numMyCols + 1] - index_offset; ++k)
133 {
134 if (col_coupleBlock_index[k] == i_row)
135 {
136 for (ic = 0; ic < col_block_size; ++ic)
137 {
138 i_Sol = ic + col_block_size * l_col;
139 for (ir = 0; ir < row_block_size; ++ir)
140 {
141 i_Equa = ir + row_block_size * l_row;
142 col_coupleBlock_val[k * block_size + ir + row_block_size * ic] +=
143 array[INDEX4
144 (i_Equa, i_Sol, k_Equa, k_Sol, num_Equa, num_Sol, NN_Equa)];
145 }
146 }
147 break;
148 }
149 }
150 }
151 }
152 }
153 }
154 }
155 }
156 }
157 else
158 {
159 mainBlock_ptr = in->mainBlock->pattern->ptr;
160 mainBlock_index = in->mainBlock->pattern->index;
161 mainBlock_val = in->mainBlock->val;
162 col_coupleBlock_ptr = in->col_coupleBlock->pattern->ptr;
163 col_coupleBlock_index = in->col_coupleBlock->pattern->index;
164 col_coupleBlock_val = in->col_coupleBlock->val;
165 row_coupleBlock_ptr = in->row_coupleBlock->pattern->ptr;
166 row_coupleBlock_index = in->row_coupleBlock->pattern->index;
167 row_coupleBlock_val = in->row_coupleBlock->val;
168
169 for (k_Equa = 0; k_Equa < NN_Equa; ++k_Equa)
170 { /* Down columns of array */
171 j_Equa = Nodes_Equa[k_Equa];
172 for (l_row = 0; l_row < num_subblocks_Equa; ++l_row)
173 {
174 i_row = j_Equa * num_subblocks_Equa + l_row;
175 /* only look at the matrix rows stored on this processor */
176 if (i_row < numMyRows)
177 {
178 for (k_Sol = 0; k_Sol < NN_Sol; ++k_Sol)
179 { /* Across rows of array */
180 j_Sol = Nodes_Sol[k_Sol];
181 for (l_col = 0; l_col < num_subblocks_Sol; ++l_col)
182 {
183 /* only look at the matrix rows stored on this processor */
184 i_col = j_Sol * num_subblocks_Sol + index_offset + l_col;
185 if (i_col < numMyCols + index_offset)
186 {
187 for (k = mainBlock_ptr[i_row] - index_offset;
188 k < mainBlock_ptr[i_row + 1] - index_offset; ++k)
189 {
190 if (mainBlock_index[k] == i_col)
191 {
192 /* Entry array(k_Sol, j_Equa) is a block (row_block_size x col_block_size) */
193 for (ic = 0; ic < col_block_size; ++ic)
194 {
195 i_Sol = ic + col_block_size * l_col;
196 for (ir = 0; ir < row_block_size; ++ir)
197 {
198 i_Equa = ir + row_block_size * l_row;
199 mainBlock_val[k * block_size + ir + row_block_size * ic] +=
200 array[INDEX4
201 (i_Equa, i_Sol, k_Equa, k_Sol, num_Equa, num_Sol, NN_Equa)];
202 }
203 }
204 break;
205 }
206 }
207 }
208 else
209 {
210 for (k = col_coupleBlock_ptr[i_row] - index_offset;
211 k < col_coupleBlock_ptr[i_row + 1] - index_offset; ++k)
212 {
213 if (col_coupleBlock_index[k] == i_col - numMyCols)
214 {
215 /* Entry array(k_Sol, j_Equa) is a block (row_block_size x col_block_size) */
216 for (ic = 0; ic < col_block_size; ++ic)
217 {
218 i_Sol = ic + col_block_size * l_col;
219 for (ir = 0; ir < row_block_size; ++ir)
220 {
221 i_Equa = ir + row_block_size * l_row;
222 col_coupleBlock_val[k * block_size + ir + row_block_size * ic] +=
223 array[INDEX4
224 (i_Equa, i_Sol, k_Equa, k_Sol, num_Equa, num_Sol, NN_Equa)];
225 }
226 }
227 break;
228 }
229 }
230 }
231 }
232 }
233 }
234 else
235 {
236 for (k_Sol = 0; k_Sol < NN_Sol; ++k_Sol)
237 { /* Across rows of array */
238 j_Sol = Nodes_Sol[k_Sol];
239 for (l_col = 0; l_col < num_subblocks_Sol; ++l_col)
240 {
241 i_col = j_Sol * num_subblocks_Sol + index_offset + l_col;
242 if (i_col < numMyCols + index_offset)
243 {
244 for (k = row_coupleBlock_ptr[i_row - numMyRows] - index_offset;
245 k < row_coupleBlock_ptr[i_row - numMyRows + 1] - index_offset; ++k)
246 {
247 if (row_coupleBlock_index[k] == i_col)
248 {
249 /* Entry array(k_Sol, j_Equa) is a block (row_block_size x col_block_size) */
250 for (ic = 0; ic < col_block_size; ++ic)
251 {
252 i_Sol = ic + col_block_size * l_col;
253 for (ir = 0; ir < row_block_size; ++ir)
254 {
255 i_Equa = ir + row_block_size * l_row;
256 row_coupleBlock_val[k * block_size + ir + row_block_size * ic] +=
257 array[INDEX4
258 (i_Equa, i_Sol, k_Equa, k_Sol, num_Equa, num_Sol, NN_Equa)];
259 }
260 }
261 break;
262 }
263 }
264 }
265 }
266 }
267 }
268 }
269 }
270 }
271 }
272
273 } // namespace dudley
274

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26