/[escript]/trunk/dudley/src/Assemble_addToSystemMatrix.c
ViewVC logotype

Contents of /trunk/dudley/src/Assemble_addToSystemMatrix.c

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26