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