/[escript]/branches/domexper/dudley/src/Assemble_addToSystemMatrix.c
ViewVC logotype

Contents of /branches/domexper/dudley/src/Assemble_addToSystemMatrix.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3086 - (show annotations)
Thu Aug 5 05:07:58 2010 UTC (9 years, 5 months ago) by jfenwick
File MIME type: text/plain
File size: 11013 byte(s)
Another pass at removing finley

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
17 /* Dudley: SystemMatrix and SystemVector */
18
19 /* adds the matrix array[Equa,Sol,NN,NN] onto the matrix in. */
20 /* the rows/columns are given by */
21 /* i_Equa+Equa*Nodes_Equa[Nodes[j_Equa]] (i_Equa=0:Equa; j_Equa=0:NN_Equa). */
22 /* the routine has to be called from a parallel region */
23
24 /* This routine assumes that in->Equa=in->Sol=1, i.e. */
25 /* array is fully packed. */
26 /* TODO: the case in->Equa!=1 */
27
28 /**************************************************************/
29
30 #include "Assemble.h"
31 #include "IndexList.h"
32
33 /**************************************************************/
34
35 void Dudley_Assemble_addToSystemMatrix(Paso_SystemMatrix* in,dim_t NN_Equa,index_t* Nodes_Equa, dim_t num_Equa,
36 dim_t NN_Sol,index_t* Nodes_Sol, dim_t num_Sol, double* array) {
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, *row_coupleBlock_index;
40 double *mainBlock_val, *row_coupleBlock_val, *col_coupleBlock_val;
41 dim_t row_block_size=in->row_block_size;
42 dim_t col_block_size=in->col_block_size;
43 dim_t block_size=in->block_size;
44 dim_t num_subblocks_Equa=num_Equa/row_block_size;
45 dim_t num_subblocks_Sol=num_Sol/col_block_size;
46 dim_t numMyCols=in->pattern->mainPattern->numInput;
47 dim_t numMyRows=in->pattern->mainPattern->numOutput;
48
49 if (in->type & MATRIX_FORMAT_CSC) {
50 /* MATRIX_FORMAT_CSC does not support MPI !!!!! */
51 mainBlock_ptr=in->mainBlock->pattern->ptr;
52 mainBlock_index=in->mainBlock->pattern->index;
53 mainBlock_val=in->mainBlock->val;
54 col_coupleBlock_ptr=in->col_coupleBlock->pattern->ptr;
55 col_coupleBlock_index=in->col_coupleBlock->pattern->index;
56 col_coupleBlock_val=in->col_coupleBlock->val;
57 row_coupleBlock_ptr=in->row_coupleBlock->pattern->ptr;
58 row_coupleBlock_index=in->row_coupleBlock->pattern->index;
59 row_coupleBlock_val=in->row_coupleBlock->val;
60
61 for (k_Sol=0;k_Sol<NN_Sol;++k_Sol) { /* Down columns of array */
62 j_Sol=Nodes_Sol[k_Sol];
63 for (l_col=0;l_col<num_subblocks_Sol;++l_col) {
64 i_col=j_Sol*num_subblocks_Sol+l_col;
65 if (i_col < numMyCols) {
66 for (k_Equa=0;k_Equa<NN_Equa;++k_Equa) { /* Across cols of array */
67 j_Equa=Nodes_Equa[k_Equa];
68 for (l_row=0;l_row<num_subblocks_Equa;++l_row) {
69 i_row=j_Equa*num_subblocks_Equa+index_offset+l_row;
70 if (i_row < numMyRows + index_offset ) {
71 for (k=mainBlock_ptr[i_col]-index_offset;k<mainBlock_ptr[i_col+1]-index_offset;++k) {
72 if (mainBlock_index[k]==i_row) {
73 /* Entry array(k_Equa, j_Sol) is a block (col_block_size x col_block_size) */
74 for (ic=0;ic<col_block_size;++ic) {
75 i_Sol=ic+col_block_size*l_col;;
76 for (ir=0;ir<row_block_size;++ir) {
77 i_Equa=ir+row_block_size*l_row;
78 mainBlock_val[k*block_size+ir+row_block_size*ic]+=
79 array[INDEX4(i_Equa,i_Sol,k_Equa,k_Sol,num_Equa,num_Sol,NN_Equa)];
80 }
81 }
82 break;
83 }
84 }
85 } else {
86 for (k=col_coupleBlock_ptr[i_col]-index_offset;k<col_coupleBlock_ptr[i_col+1]-index_offset;++k) {
87 if (row_coupleBlock_index[k] == i_row-numMyRows) {
88 for (ic=0;ic<col_block_size;++ic) {
89 i_Sol=ic+col_block_size*l_col;
90 for (ir=0;ir<row_block_size;++ir) {
91 i_Equa=ir+row_block_size*l_row;
92 row_coupleBlock_val[k*block_size+ir+row_block_size*ic]+=
93 array[INDEX4(i_Equa,i_Sol,k_Equa,k_Sol,num_Equa,num_Sol,NN_Equa)];;
94 }
95 }
96 break;
97 }
98 }
99 }
100 }
101 }
102 } else {
103 for (k_Equa=0;k_Equa<NN_Equa;++k_Equa) { /* Across rows of array */
104 j_Equa=Nodes_Equa[k_Equa];
105 for (l_row=0;l_row<num_subblocks_Equa;++l_row) {
106 i_row=j_Equa*num_subblocks_Equa+index_offset+l_row;
107 if (i_row < numMyRows + index_offset ) {
108 for (k=col_coupleBlock_ptr[i_col-numMyCols]-index_offset;k<col_coupleBlock_ptr[i_col-numMyCols+1]-index_offset;++k) {
109 if (col_coupleBlock_index[k] == i_row) {
110 for (ic=0;ic<col_block_size;++ic) {
111 i_Sol=ic+col_block_size*l_col;
112 for (ir=0;ir<row_block_size;++ir) {
113 i_Equa=ir+row_block_size*l_row;
114 col_coupleBlock_val[k*block_size+ir+row_block_size*ic]+=
115 array[INDEX4(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 }
126 }
127 } else if (in->type & MATRIX_FORMAT_TRILINOS_CRS) {
128 /* this needs to be modified */
129 #ifdef TRILINOS
130 for (k_Equa=0;k_Equa<NN_Equa;++k_Equa) { /* Down columns of array */
131 j_Equa=Nodes_Equa[k_Equa];
132 if (j_Equa < in->mainBlock->pattern->output_node_distribution->numLocal) {
133 for (k_Sol=0;k_Sol<NN_Sol;++k_Sol) { /* Across rows of array */
134 j_Sol=Nodes_Sol[k_Sol];
135 for (l_row=0;l_row<num_subblocks_Equa;++l_row) {
136 irow=j_Equa*row_block_size+l_row;
137 for (l_col=0;l_col<col_block_size;++l_col) {
138 icol=j_Sol*col_block_size+index_offset+l_col;
139 /* irow is local and icol is global */
140 Trilinos_SumIntoMyValues(in->trilinos_data, irow, icol, array[INDEX4(l_row,l_col,k_Equa,k_Sol,num_Equa,num_Sol,NN_Equa)]);
141 }
142 }
143 }
144 }
145 }
146 #endif
147 } else {
148 mainBlock_ptr=in->mainBlock->pattern->ptr;
149 mainBlock_index=in->mainBlock->pattern->index;
150 mainBlock_val=in->mainBlock->val;
151 col_coupleBlock_ptr=in->col_coupleBlock->pattern->ptr;
152 col_coupleBlock_index=in->col_coupleBlock->pattern->index;
153 col_coupleBlock_val=in->col_coupleBlock->val;
154 row_coupleBlock_ptr=in->row_coupleBlock->pattern->ptr;
155 row_coupleBlock_index=in->row_coupleBlock->pattern->index;
156 row_coupleBlock_val=in->row_coupleBlock->val;
157
158 for (k_Equa=0;k_Equa<NN_Equa;++k_Equa) { /* Down columns of array */
159 j_Equa=Nodes_Equa[k_Equa];
160 for (l_row=0;l_row<num_subblocks_Equa;++l_row) {
161 i_row=j_Equa*num_subblocks_Equa+l_row;
162 /* only look at the matrix rows stored on this processor */
163 if (i_row < numMyRows) {
164 for (k_Sol=0;k_Sol<NN_Sol;++k_Sol) { /* Across rows of array */
165 j_Sol=Nodes_Sol[k_Sol];
166 for (l_col=0;l_col<num_subblocks_Sol;++l_col) {
167 /* only look at the matrix rows stored on this processor */
168 i_col=j_Sol*num_subblocks_Sol+index_offset+l_col;
169 if (i_col < numMyCols + index_offset ) {
170 for (k=mainBlock_ptr[i_row]-index_offset;k<mainBlock_ptr[i_row+1]-index_offset;++k) {
171 if (mainBlock_index[k]==i_col) {
172 /* Entry array(k_Sol, j_Equa) is a block (row_block_size x col_block_size) */
173 for (ic=0;ic<col_block_size;++ic) {
174 i_Sol=ic+col_block_size*l_col;
175 for (ir=0;ir<row_block_size;++ir) {
176 i_Equa=ir+row_block_size*l_row;
177 mainBlock_val[k*block_size+ir+row_block_size*ic]+=
178 array[INDEX4(i_Equa,i_Sol,k_Equa,k_Sol,num_Equa,num_Sol,NN_Equa)];
179 }
180 }
181 break;
182 }
183 }
184 } else {
185 for (k=col_coupleBlock_ptr[i_row]-index_offset;k<col_coupleBlock_ptr[i_row+1]-index_offset;++k) {
186 if (col_coupleBlock_index[k] == i_col-numMyCols) {
187 /* Entry array(k_Sol, j_Equa) is a block (row_block_size x col_block_size) */
188 for (ic=0;ic<col_block_size;++ic) {
189 i_Sol=ic+col_block_size*l_col;
190 for (ir=0;ir<row_block_size;++ir) {
191 i_Equa=ir+row_block_size*l_row;
192 col_coupleBlock_val[k*block_size+ir+row_block_size*ic]+=
193 array[INDEX4(i_Equa,i_Sol,k_Equa,k_Sol,num_Equa,num_Sol,NN_Equa)];
194 }
195 }
196 break;
197 }
198 }
199 }
200 }
201 }
202 } else {
203 for (k_Sol=0;k_Sol<NN_Sol;++k_Sol) { /* Across rows of array */
204 j_Sol=Nodes_Sol[k_Sol];
205 for (l_col=0;l_col<num_subblocks_Sol;++l_col) {
206 i_col=j_Sol*num_subblocks_Sol+index_offset+l_col;
207 if (i_col < numMyCols + index_offset ) {
208 for (k=row_coupleBlock_ptr[i_row-numMyRows]-index_offset;k<row_coupleBlock_ptr[i_row-numMyRows+1]-index_offset;++k) {
209 if (row_coupleBlock_index[k] == i_col) {
210 /* Entry array(k_Sol, j_Equa) is a block (row_block_size x col_block_size) */
211 for (ic=0;ic<col_block_size;++ic) {
212 i_Sol=ic+col_block_size*l_col;
213 for (ir=0;ir<row_block_size;++ir) {
214 i_Equa=ir+row_block_size*l_row;
215 row_coupleBlock_val[k*block_size+ir+row_block_size*ic]+=
216 array[INDEX4(i_Equa,i_Sol,k_Equa,k_Sol,num_Equa,num_Sol,NN_Equa)];
217 }
218 }
219 break;
220 }
221 }
222 }
223 }
224 }
225 }
226 }
227 }
228 }
229 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26