/[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 1223 - (show annotations)
Fri Aug 3 02:40:39 2007 UTC (12 years, 9 months ago) by gross
File MIME type: text/plain
File size: 7943 byte(s)
first attemt towards an improved MPI version.  

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,i_row, i_col;
41 index_t *mainBlock_ptr, *mainBlock_index, *coupleBlock_ptr, *coupleBlock_index;
42 double *mainBlock_val, *coupleBlock_val;
43 dim_t row_block_size=in->row_block_size;
44 dim_t col_block_size=in->col_block_size;
45 dim_t block_size=in->block_size;
46 dim_t num_subblocks_Equa=num_Equa/row_block_size;
47 dim_t num_subblocks_Sol=num_Sol/col_block_size;
48 dim_t numMyCols=in->pattern->mainPattern->numInput;
49 dim_t numMyRows=in->pattern->mainPattern->numOutput;
50
51
52 if (in->type & MATRIX_FORMAT_CSC) {
53 /* MATRIX_FORMAT_CSC does not support MPI !!!!! */
54 mainBlock_ptr=in->mainBlock->pattern->ptr;
55 mainBlock_index=in->mainBlock->pattern->index;
56 mainBlock_val=in->mainBlock->val;
57 for (k_Sol=0;k_Sol<NN_Sol;k_Sol++) {
58 j_Sol=Nodes_Sol[k_Sol];
59 for (l_col=0;l_col<num_subblocks_Sol;++l_col) {
60 i_col=j_Sol*num_subblocks_Sol+l_col;
61 for (k_Equa=0;k_Equa<NN_Equa;k_Equa++) {
62 j_Equa=Nodes_Equa[k_Equa];
63 for (l_row=0;l_row<num_subblocks_Equa;++l_row) {
64 i_row=j_Equa*num_subblocks_Equa+index_offset+l_row;
65 for (k=mainBlock_ptr[i_col]-index_offset; k<mainBlock_ptr[i_col+1]-index_offset;++k) {
66 if (mainBlock_index[k] == i_row) {
67 for (ic=0;ic<col_block_size;++ic) {
68 i_Sol=ic+col_block_size*l_col;
69 for (ir=0;ir<row_block_size;++ir) {
70 i_Equa=ir+row_block_size*l_row;
71 mainBlock_val[k*block_size+ir+row_block_size*ic]+=
72 array[INDEX4(i_Equa,i_Sol,k_Equa,k_Sol,num_Equa,num_Sol,NN_Equa)];
73 }
74 }
75 break;
76 }
77 }
78 }
79 }
80 }
81 }
82 } else if (in->type & MATRIX_FORMAT_TRILINOS_CRS) {
83 /* this needs to be modified */
84 #ifdef TRILINOS
85 for (k_Equa=0;k_Equa<NN_Equa;++k_Equa) { /* Down columns of array */
86 j_Equa=Nodes_Equa[k_Equa];
87 if (j_Equa < in->mainBlock->pattern->output_node_distribution->numLocal) {
88 for (k_Sol=0;k_Sol<NN_Sol;++k_Sol) { /* Across rows of array */
89 j_Sol=Nodes_Sol[k_Sol];
90 for (l_row=0;l_row<num_subblocks_Equa;++l_row) {
91 irow=j_Equa*row_block_size+l_row;
92 for (l_col=0;l_col<col_block_size;++l_col) {
93 icol=j_Sol*col_block_size+index_offset+l_col;
94 // irow is local and icol is global
95 Trilinos_SumIntoMyValues(in->trilinos_data, irow, icol, array[INDEX4(l_row,l_col,k_Equa,k_Sol,num_Equa,num_Sol,NN_Equa)]);
96 }
97 }
98 }
99 }
100 }
101 #endif
102 } else {
103
104 mainBlock_ptr=in->mainBlock->pattern->ptr;
105 mainBlock_index=in->mainBlock->pattern->index;
106 mainBlock_val=in->mainBlock->val;
107 coupleBlock_ptr=in->coupleBlock->pattern->ptr;
108 coupleBlock_index=in->coupleBlock->pattern->index;
109 coupleBlock_val=in->coupleBlock->val;
110
111 for (k_Equa=0;k_Equa<NN_Equa;++k_Equa) { /* Down columns of array */
112 j_Equa=Nodes_Equa[k_Equa];
113 for (l_row=0;l_row<num_subblocks_Equa;++l_row) {
114 i_row=j_Equa*num_subblocks_Equa+l_row;
115 /* only look at the matrix rows stored on this processor */
116 if (i_row < numMyRows) {
117 for (k_Sol=0;k_Sol<NN_Sol;++k_Sol) { /* Across rows of array */
118 j_Sol=Nodes_Sol[k_Sol];
119 for (l_col=0;l_col<num_subblocks_Sol;++l_col) {
120 /* only look at the matrix rows stored on this processor */
121 i_col=j_Sol*num_subblocks_Sol+index_offset+l_col;
122 if (i_col < numMyCols + index_offset ) {
123 for (k=mainBlock_ptr[i_row]-index_offset;k<mainBlock_ptr[i_row+1]-index_offset;++k) {
124 if (mainBlock_index[k]==i_col) {
125 /* Entry array(k_Sol, j_Equa) is a block (row_block_size x col_block_size) */
126 for (ic=0;ic<col_block_size;++ic) {
127 i_Sol=ic+col_block_size*l_col;
128 for (ir=0;ir<row_block_size;++ir) {
129 i_Equa=ir+row_block_size*l_row;
130 mainBlock_val[k*block_size+ir+row_block_size*ic]+=
131 array[INDEX4(i_Equa,i_Sol,k_Equa,k_Sol,num_Equa,num_Sol,NN_Equa)];
132 }
133 }
134 break;
135 }
136 }
137 } else {
138 for (k=coupleBlock_ptr[i_row]-index_offset;k<coupleBlock_ptr[i_row+1]-index_offset;++k) {
139 if (coupleBlock_index[k] == i_col-numMyCols) {
140 /* Entry array(k_Sol, j_Equa) is a block (row_block_size x col_block_size) */
141 for (ic=0;ic<col_block_size;++ic) {
142 i_Sol=ic+col_block_size*l_col;
143 for (ir=0;ir<row_block_size;++ir) {
144 i_Equa=ir+row_block_size*l_row;
145 coupleBlock_val[k*block_size+ir+row_block_size*ic]+=
146 array[INDEX4(i_Equa,i_Sol,k_Equa,k_Sol,num_Equa,num_Sol,NN_Equa)];
147 }
148 }
149 break;
150 }
151 }
152 }
153 }
154 }
155 } /* end i_row check */
156 }
157 }
158 }
159 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26