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

Annotation of /trunk/finley/src/Assemble_addToSystemMatrix.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3071 - (hide annotations)
Wed Jul 21 05:37:30 2010 UTC (11 years, 3 months ago) by gross
File MIME type: text/plain
File size: 11013 byte(s)
new darcy flux
1 ksteube 1312
2     /*******************************************************
3 ksteube 1811 *
4 jfenwick 2881 * Copyright (c) 2003-2010 by University of Queensland
5 ksteube 1811 * 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 ksteube 1312
14 ksteube 1811
15 jgs 150 /**************************************************************/
16    
17     /* Finley: 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 ksteube 1312 #include "IndexList.h"
32 jgs 150
33     /**************************************************************/
34    
35     void Finley_Assemble_addToSystemMatrix(Paso_SystemMatrix* in,dim_t NN_Equa,index_t* Nodes_Equa, dim_t num_Equa,
36 ksteube 1312 dim_t NN_Sol,index_t* Nodes_Sol, dim_t num_Sol, double* array) {
37 gross 416 index_t index_offset=(in->type & MATRIX_FORMAT_OFFSET1 ? 1:0);
38 phornby 1628 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 gross 1552 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 jgs 150 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 ksteube 1312 dim_t numMyCols=in->pattern->mainPattern->numInput;
47     dim_t numMyRows=in->pattern->mainPattern->numOutput;
48 jgs 150
49 gross 416 if (in->type & MATRIX_FORMAT_CSC) {
50 ksteube 1312 /* MATRIX_FORMAT_CSC does not support MPI !!!!! */
51 gross 3071 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 ksteube 1312 } 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 ksteube 1322 /* irow is local and icol is global */
140 ksteube 1312 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 gross 416 } else {
148 ksteube 1312 mainBlock_ptr=in->mainBlock->pattern->ptr;
149     mainBlock_index=in->mainBlock->pattern->index;
150     mainBlock_val=in->mainBlock->val;
151 gross 1552 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 ksteube 1312
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 gross 1552 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 ksteube 1312 /* 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 gross 1552 col_coupleBlock_val[k*block_size+ir+row_block_size*ic]+=
193 ksteube 1312 array[INDEX4(i_Equa,i_Sol,k_Equa,k_Sol,num_Equa,num_Sol,NN_Equa)];
194     }
195     }
196     break;
197     }
198     }
199     }
200 jgs 150 }
201 ksteube 1312 }
202 gross 1552 } 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 ksteube 1312 }
227     }
228 jgs 150 }
229     }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26