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

Contents of /trunk-mpi-branch/finley/src/Assemble_addToSystemMatrix.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1087 - (show annotations)
Thu Apr 12 10:01:47 2007 UTC (13 years, 1 month ago) by gross
File MIME type: text/plain
File size: 8691 byte(s)
the MPI version of PASO.PCG is running. 
There is a bug in the rectangular mesh generators but they need to be
revised in any case to clean up the code.


1 /*
2 ************************************************************
3 * Copyright 2006 by ACcESS MNRF *
4 * *
5 * http://www.access.edu.au *
6 * Primary Business: Queensland, Australia *
7 * Licensed under the Open Software License version 3.0 *
8 * http://www.opensource.org/licenses/osl-3.0.php *
9 * *
10 ************************************************************
11 */
12 /**************************************************************/
13
14 /* Finley: SystemMatrix and SystemVector */
15
16 /* adds the matrix array[Equa,Sol,NN,NN] onto the matrix in. */
17 /* the rows/columns are given by */
18 /* i_Equa+Equa*Nodes_Equa[Nodes[j_Equa]] (i_Equa=0:Equa; j_Equa=0:NN_Equa). */
19 /* the routine has to be called from a parallel region */
20
21 /* This routine assumes that in->Equa=in->Sol=1, i.e. */
22 /* array is fully packed. */
23 /* TODO: the case in->Equa!=1 */
24
25 /**************************************************************/
26
27 /* Author: gross@access.edu.au */
28 /* Version: $Id$ */
29
30 /**************************************************************/
31
32 #include "Assemble.h"
33 #include "IndexList.h"
34
35 /**************************************************************/
36
37 void Finley_Assemble_addToSystemMatrix(Paso_SystemMatrix* in,dim_t NN_Equa,index_t* Nodes_Equa, dim_t num_Equa,
38 dim_t NN_Sol,index_t* Nodes_Sol, dim_t num_Sol, double* array) {
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,index,k,iptr, irow, icol;
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
47 if (in->type & MATRIX_FORMAT_CSC) {
48 for (k_Sol=0;k_Sol<NN_Sol;k_Sol++) {
49 j_Sol=Nodes_Sol[k_Sol];
50 for (l_col=0;l_col<num_subblocks_Sol;++l_col) {
51 iptr=j_Sol*num_subblocks_Sol+l_col;
52 for (k_Equa=0;k_Equa<NN_Equa;k_Equa++) {
53 j_Equa=Nodes_Equa[k_Equa];
54 for (l_row=0;l_row<num_subblocks_Equa;++l_row) {
55 index=j_Equa*num_subblocks_Equa+index_offset+l_row;
56 for (k=in->pattern->ptr[iptr]-index_offset;k<in->pattern->ptr[iptr+1]-index_offset;++k) {
57 if (in->pattern->index[k]==index) {
58 for (ic=0;ic<col_block_size;++ic) {
59 i_Sol=ic+col_block_size*l_col;
60 for (ir=0;ir<row_block_size;++ir) {
61 i_Equa=ir+row_block_size*l_row;
62 in->val[k*block_size+ir+row_block_size*ic]+=
63 array[INDEX4(i_Equa,i_Sol,k_Equa,k_Sol,num_Equa,num_Sol,NN_Equa)];
64 }
65 }
66 break;
67 }
68 }
69 }
70 }
71 }
72 }
73 } else if (in->type & MATRIX_FORMAT_TRILINOS_CRS) {
74 #ifdef TRILINOS
75 for (k_Equa=0;k_Equa<NN_Equa;++k_Equa) { /* Down columns of array */
76 j_Equa=Nodes_Equa[k_Equa];
77 if (j_Equa < in->pattern->output_node_distribution->numLocal) {
78 for (k_Sol=0;k_Sol<NN_Sol;++k_Sol) { /* Across rows of array */
79 j_Sol=Nodes_Sol[k_Sol];
80 j_Sol = Finley_IndexList_localToGlobal(in->pattern->input_node_distribution, j_Sol);
81 for (l_row=0;l_row<num_subblocks_Equa;++l_row) {
82 irow=j_Equa*row_block_size+l_row;
83 for (l_col=0;l_col<col_block_size;++l_col) {
84 icol=j_Sol*col_block_size+index_offset+l_col;
85 // irow is local and icol is global
86 Trilinos_SumIntoMyValues(in->trilinos_data, irow, icol, array[INDEX4(l_row,l_col,k_Equa,k_Sol,num_Equa,num_Sol,NN_Equa)]);
87 }
88 }
89 }
90 }
91 }
92 #endif
93 } else {
94 if (in->mpi_info->size == 1) {
95 for (k_Equa=0;k_Equa<NN_Equa;++k_Equa) { /* Down columns of array */
96 j_Equa=Nodes_Equa[k_Equa];
97 for (l_row=0;l_row<num_subblocks_Equa;++l_row) {
98 iptr=j_Equa*num_subblocks_Equa+l_row;
99 for (k_Sol=0;k_Sol<NN_Sol;++k_Sol) { /* Across rows of array */
100 j_Sol=Nodes_Sol[k_Sol];
101 for (l_col=0;l_col<num_subblocks_Sol;++l_col) {
102 index=j_Sol*num_subblocks_Sol+index_offset+l_col;
103 for (k=in->pattern->ptr[iptr]-index_offset;k<in->pattern->ptr[iptr+1]-index_offset;++k) {
104 if (in->pattern->index[k]==index) {
105 for (ic=0;ic<col_block_size;++ic) { /* Entry array(k_Sol, j_Equa) is a block (row_block_size x col_block_size) */
106 i_Sol=ic+col_block_size*l_col;
107 for (ir=0;ir<row_block_size;++ir) {
108 i_Equa=ir+row_block_size*l_row;
109 in->val[k*block_size+ir+row_block_size*ic]+=
110 array[INDEX4(i_Equa,i_Sol,k_Equa,k_Sol,num_Equa,num_Sol,NN_Equa)];
111 }
112 }
113 break;
114 }
115 }
116 }
117 }
118 }
119 }
120 } else {
121 for (k_Equa=0;k_Equa<NN_Equa;++k_Equa) { /* Down columns of array */
122 j_Equa=Nodes_Equa[k_Equa];
123 if (j_Equa < in->pattern->output_node_distribution->numLocal) {
124 for (l_row=0;l_row<num_subblocks_Equa;++l_row) {
125 iptr=j_Equa*num_subblocks_Equa+l_row;
126 for (k_Sol=0;k_Sol<NN_Sol;++k_Sol) { /* Across rows of array */
127 j_Sol=Nodes_Sol[k_Sol];
128 j_Sol = Finley_IndexList_localToGlobal(in->pattern->input_node_distribution, j_Sol);
129 for (l_col=0;l_col<num_subblocks_Sol;++l_col) {
130 index=j_Sol*num_subblocks_Sol+index_offset+l_col;
131 for (k=in->pattern->ptr[iptr]-index_offset;k<in->pattern->ptr[iptr+1]-index_offset;++k) {
132 if (in->pattern->index[k]==index) {
133 for (ic=0;ic<col_block_size;++ic) { /* Entry array(k_Sol, j_Equa) is a block (row_block_size x col_block_size) */
134 i_Sol=ic+col_block_size*l_col;
135 for (ir=0;ir<row_block_size;++ir) {
136 i_Equa=ir+row_block_size*l_row;
137 in->val[k*block_size+ir+row_block_size*ic]+=
138 array[INDEX4(i_Equa,i_Sol,k_Equa,k_Sol,num_Equa,num_Sol,NN_Equa)];
139 }
140 }
141 break;
142 }
143 }
144 }
145 }
146 }
147 }
148 }
149 }
150 }
151 }
152 /*
153 * $Log$
154 * Revision 1.2 2005/09/15 03:44:21 jgs
155 * Merge of development branch dev-02 back to main trunk on 2005-09-15
156 *
157 * Revision 1.1.2.1 2005/09/07 09:47:29 gross
158 * missing file
159 *
160 * Revision 1.6 2005/07/08 04:07:58 jgs
161 * Merge of development branch back to main trunk on 2005-07-08
162 *
163 * Revision 1.1.1.1.2.3 2005/06/29 02:34:56 gross
164 * some changes towards 64 integers in finley
165 *
166 * Revision 1.1.1.1.2.2 2005/03/15 07:23:55 gross
167 * Finley's interface to the SCSL library can deal with systems of PDEs now. tests shows that the SCSL library cannot deal with problems with more then 200000 unknowns. problem has been reported to SGI.
168 *
169 * Revision 1.1.1.1.2.1 2004/11/12 06:58:19 gross
170 * a lot of changes to get the linearPDE class running: most important change is that there is no matrix format exposed to the user anymore. the format is chosen by the Domain according to the solver and symmetry
171 *
172 * Revision 1.1.1.1 2004/10/26 06:53:57 jgs
173 * initial import of project esys2
174 *
175 * Revision 1.1 2004/07/02 04:21:13 gross
176 * Finley C code has been included
177 *
178 *
179 */

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26