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