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