/[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 969 - (hide annotations)
Tue Feb 13 23:02:23 2007 UTC (12 years, 9 months ago) by ksteube
File MIME type: text/plain
File size: 7004 byte(s)
Parallelization using MPI for solution of implicit problems.

Parallelization for explicit problems has already been accomplished in
the main SVN branch.

This is incomplete and is not ready for use.


1 jgs 150 /*
2 elspeth 616 ************************************************************
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 jgs 150 */
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    
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 gross 416 index_t index_offset=(in->type & MATRIX_FORMAT_OFFSET1 ? 1:0);
39 jgs 150 dim_t k_Equa,j_Equa,j_Sol,k_Sol,i_Equa,i_Sol,l_col,l_row,ic,ir,index,k,iptr;
40     dim_t row_block_size=in->row_block_size;
41     dim_t col_block_size=in->col_block_size;
42     dim_t block_size=in->block_size;
43     dim_t num_subblocks_Equa=num_Equa/row_block_size;
44     dim_t num_subblocks_Sol=num_Sol/col_block_size;
45    
46 ksteube 969 printf("ksteube addToSystemMatrix NN_Sol=%d num_subblocks_Sol=%d NN_Equa=%d num_subblocks_Equa=%d\n", NN_Sol, num_subblocks_Sol, NN_Equa, num_subblocks_Equa);
47 gross 416 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 jgs 150 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 gross 416 }
72     }
73 ksteube 969 } else if (in->type & MATRIX_FORMAT_TRILINOS_CRS) {
74     #ifdef PASO_MPI
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 < row_degreeOfFreedomDistribution->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(col_degreeOfFreedomDistribution, my_CPU, 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_blocksize;++l_col) {
84     icol=j_Sol*col_blocksize+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 gross 416 } else {
94 ksteube 969 for (k_Equa=0;k_Equa<NN_Equa;++k_Equa) { /* Down columns of array */
95 gross 416 j_Equa=Nodes_Equa[k_Equa];
96     for (l_row=0;l_row<num_subblocks_Equa;++l_row) {
97     iptr=j_Equa*num_subblocks_Equa+l_row;
98 ksteube 969 for (k_Sol=0;k_Sol<NN_Sol;++k_Sol) { /* Across rows of array */
99 gross 416 j_Sol=Nodes_Sol[k_Sol];
100     for (l_col=0;l_col<num_subblocks_Sol;++l_col) {
101     index=j_Sol*num_subblocks_Sol+index_offset+l_col;
102     for (k=in->pattern->ptr[iptr]-index_offset;k<in->pattern->ptr[iptr+1]-index_offset;++k) {
103 jgs 150 if (in->pattern->index[k]==index) {
104 ksteube 969 for (ic=0;ic<col_block_size;++ic) { /* Entry array(k_Sol, j_Equa) is a block (row_block_size x col_block_size) */
105 jgs 150 i_Sol=ic+col_block_size*l_col;
106     for (ir=0;ir<row_block_size;++ir) {
107     i_Equa=ir+row_block_size*l_row;
108     in->val[k*block_size+ir+row_block_size*ic]+=
109     array[INDEX4(i_Equa,i_Sol,k_Equa,k_Sol,num_Equa,num_Sol,NN_Equa)];
110 ksteube 969 /* printf("ksteube assigning val[ %d (%d,%d) ] += %f \n", k*block_size+ir+row_block_size*ic, k_Equa, k_Sol, array[INDEX4(i_Equa,i_Sol,k_Equa,k_Sol,num_Equa,num_Sol,NN_Equa)]); */
111 jgs 150 }
112     }
113     break;
114     }
115     }
116     }
117     }
118 gross 416 }
119     }
120 jgs 150 }
121     }
122     /*
123     * $Log$
124     * Revision 1.2 2005/09/15 03:44:21 jgs
125     * Merge of development branch dev-02 back to main trunk on 2005-09-15
126     *
127     * Revision 1.1.2.1 2005/09/07 09:47:29 gross
128     * missing file
129     *
130     * Revision 1.6 2005/07/08 04:07:58 jgs
131     * Merge of development branch back to main trunk on 2005-07-08
132     *
133     * Revision 1.1.1.1.2.3 2005/06/29 02:34:56 gross
134     * some changes towards 64 integers in finley
135     *
136     * Revision 1.1.1.1.2.2 2005/03/15 07:23:55 gross
137     * 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.
138     *
139     * Revision 1.1.1.1.2.1 2004/11/12 06:58:19 gross
140     * 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
141     *
142     * Revision 1.1.1.1 2004/10/26 06:53:57 jgs
143     * initial import of project esys2
144     *
145     * Revision 1.1 2004/07/02 04:21:13 gross
146     * Finley C code has been included
147     *
148     *
149     */

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26