/[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 3259 - (show annotations)
Mon Oct 11 01:48:14 2010 UTC (9 years, 3 months ago) by jfenwick
File MIME type: text/plain
File size: 9924 byte(s)
Merging dudley and scons updates from branches

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26