revision 616 by elspeth, Wed Mar 22 02:46:56 2006 UTC revision 1811 by ksteube, Thu Sep 25 23:11:13 2008 UTC
# Line 1  Line 1
1  /*
2   ************************************************************  /*******************************************************
3   *          Copyright 2006 by ACcESS MNRF                   *  *
4   *                                                          *  * Copyright (c) 2003-2008 by University of Queensland
5   *              http://www.access.edu.au                    *  * Earth Systems Science Computational Center (ESSCC)
6   *       Primary Business: Queensland, Australia            *  * http://www.uq.edu.au/esscc
7   *  Licensed under the Open Software License version 3.0    *  *
9   *                                                          *  * Licensed under the Open Software License version 3.0
11  */  *
12    *******************************************************/
13
14
15  /**************************************************************/  /**************************************************************/
16
17  /* Finley: SystemMatrix and SystemVector */  /* Finley: SystemMatrix and SystemVector */
# Line 24  Line 27
27
28  /**************************************************************/  /**************************************************************/
29
/*  Author: gross@access.edu.au */
/*  Version: \$Id\$ */

/**************************************************************/

30  #include "Assemble.h"  #include "Assemble.h"
31    #include "IndexList.h"
32
33  /**************************************************************/  /**************************************************************/
34
35  void  Finley_Assemble_addToSystemMatrix(Paso_SystemMatrix* in,dim_t NN_Equa,index_t* Nodes_Equa, dim_t num_Equa,  void  Finley_Assemble_addToSystemMatrix(Paso_SystemMatrix* in,dim_t NN_Equa,index_t* Nodes_Equa, dim_t num_Equa,
36                                                        dim_t NN_Sol,index_t* Nodes_Sol, dim_t num_Sol, double* array) {                                          dim_t NN_Sol,index_t* Nodes_Sol, dim_t num_Sol, double* array) {
37    index_t index_offset=(in->type & MATRIX_FORMAT_OFFSET1 ? 1:0);    index_t index_offset=(in->type & MATRIX_FORMAT_OFFSET1 ? 1:0);
38    dim_t k_Equa,j_Equa,j_Sol,k_Sol,i_Equa,i_Sol,l_col,l_row,ic,ir,index,k,iptr;    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      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    dim_t row_block_size=in->row_block_size;    dim_t row_block_size=in->row_block_size;
42    dim_t col_block_size=in->col_block_size;    dim_t col_block_size=in->col_block_size;
43    dim_t block_size=in->block_size;    dim_t block_size=in->block_size;
44    dim_t num_subblocks_Equa=num_Equa/row_block_size;    dim_t num_subblocks_Equa=num_Equa/row_block_size;
45    dim_t num_subblocks_Sol=num_Sol/col_block_size;    dim_t num_subblocks_Sol=num_Sol/col_block_size;
46      dim_t numMyCols=in->pattern->mainPattern->numInput;
47      dim_t numMyRows=in->pattern->mainPattern->numOutput;
48
49
50    if (in->type & MATRIX_FORMAT_CSC) {    if (in->type & MATRIX_FORMAT_CSC) {
51             /* MATRIX_FORMAT_CSC does not support MPI !!!!! */
52             mainBlock_ptr=in->mainBlock->pattern->ptr;
53             mainBlock_index=in->mainBlock->pattern->index;
54             mainBlock_val=in->mainBlock->val;
55           for (k_Sol=0;k_Sol<NN_Sol;k_Sol++) {           for (k_Sol=0;k_Sol<NN_Sol;k_Sol++) {
56              j_Sol=Nodes_Sol[k_Sol];              j_Sol=Nodes_Sol[k_Sol];
57              for (l_col=0;l_col<num_subblocks_Sol;++l_col) {              for (l_col=0;l_col<num_subblocks_Sol;++l_col) {
58                 iptr=j_Sol*num_subblocks_Sol+l_col;                 i_col=j_Sol*num_subblocks_Sol+l_col;
59                 for (k_Equa=0;k_Equa<NN_Equa;k_Equa++) {                 for (k_Equa=0;k_Equa<NN_Equa;k_Equa++) {
60                   j_Equa=Nodes_Equa[k_Equa];                   j_Equa=Nodes_Equa[k_Equa];
61                   for (l_row=0;l_row<num_subblocks_Equa;++l_row) {                   for (l_row=0;l_row<num_subblocks_Equa;++l_row) {
62                      index=j_Equa*num_subblocks_Equa+index_offset+l_row;                      i_row=j_Equa*num_subblocks_Equa+index_offset+l_row;
63                  for (k=in->pattern->ptr[iptr]-index_offset;k<in->pattern->ptr[iptr+1]-index_offset;++k) {                  for (k=mainBlock_ptr[i_col]-index_offset; k<mainBlock_ptr[i_col+1]-index_offset;++k) {
64                      if (in->pattern->index[k]==index) {                      if (mainBlock_index[k] == i_row) {
65                            for (ic=0;ic<col_block_size;++ic) {                            for (ic=0;ic<col_block_size;++ic) {
66                                  i_Sol=ic+col_block_size*l_col;                                  i_Sol=ic+col_block_size*l_col;
67                                  for (ir=0;ir<row_block_size;++ir) {                                  for (ir=0;ir<row_block_size;++ir) {
68                                     i_Equa=ir+row_block_size*l_row;                                     i_Equa=ir+row_block_size*l_row;
69                             in->val[k*block_size+ir+row_block_size*ic]+=                             mainBlock_val[k*block_size+ir+row_block_size*ic]+=
70                                             array[INDEX4(i_Equa,i_Sol,k_Equa,k_Sol,num_Equa,num_Sol,NN_Equa)];                                             array[INDEX4(i_Equa,i_Sol,k_Equa,k_Sol,num_Equa,num_Sol,NN_Equa)];
71                                  }                                  }
72                            }                            }
77                 }                 }
78              }              }
79           }           }
80     } else {     } else if (in->type & MATRIX_FORMAT_TRILINOS_CRS) {
81            for (k_Equa=0;k_Equa<NN_Equa;++k_Equa) {         /* this needs to be modified */
82           #ifdef TRILINOS
83              for (k_Equa=0;k_Equa<NN_Equa;++k_Equa) { /* Down columns of array */
84              j_Equa=Nodes_Equa[k_Equa];              j_Equa=Nodes_Equa[k_Equa];
85              for (l_row=0;l_row<num_subblocks_Equa;++l_row) {          if (j_Equa < in->mainBlock->pattern->output_node_distribution->numLocal) {
86                 iptr=j_Equa*num_subblocks_Equa+l_row;                for (k_Sol=0;k_Sol<NN_Sol;++k_Sol) { /* Across rows of array */
87                 for (k_Sol=0;k_Sol<NN_Sol;++k_Sol) {                  j_Sol=Nodes_Sol[k_Sol];
88                   j_Sol=Nodes_Sol[k_Sol];                  for (l_row=0;l_row<num_subblocks_Equa;++l_row) {
89                   for (l_col=0;l_col<num_subblocks_Sol;++l_col) {                    irow=j_Equa*row_block_size+l_row;
90                      index=j_Sol*num_subblocks_Sol+index_offset+l_col;                    for (l_col=0;l_col<col_block_size;++l_col) {
91                  for (k=in->pattern->ptr[iptr]-index_offset;k<in->pattern->ptr[iptr+1]-index_offset;++k) {                       icol=j_Sol*col_block_size+index_offset+l_col;
92                      if (in->pattern->index[k]==index) {               /* irow is local and icol is global */
93                            for (ic=0;ic<col_block_size;++ic) {               Trilinos_SumIntoMyValues(in->trilinos_data, irow, icol, array[INDEX4(l_row,l_col,k_Equa,k_Sol,num_Equa,num_Sol,NN_Equa)]);
94                                  i_Sol=ic+col_block_size*l_col;                }
95                                  for (ir=0;ir<row_block_size;++ir) {                  }
96                                     i_Equa=ir+row_block_size*l_row;                }
97                             in->val[k*block_size+ir+row_block_size*ic]+=              }
98                                             array[INDEX4(i_Equa,i_Sol,k_Equa,k_Sol,num_Equa,num_Sol,NN_Equa)];            }
99                                  }         #endif
100                            }     } else {
101                            break;
102             mainBlock_ptr=in->mainBlock->pattern->ptr;
103             mainBlock_index=in->mainBlock->pattern->index;
104             mainBlock_val=in->mainBlock->val;
105             col_coupleBlock_ptr=in->col_coupleBlock->pattern->ptr;
106             col_coupleBlock_index=in->col_coupleBlock->pattern->index;
107             col_coupleBlock_val=in->col_coupleBlock->val;
108             row_coupleBlock_ptr=in->row_coupleBlock->pattern->ptr;
109             row_coupleBlock_index=in->row_coupleBlock->pattern->index;
110             row_coupleBlock_val=in->row_coupleBlock->val;
111
112             for (k_Equa=0;k_Equa<NN_Equa;++k_Equa) { /* Down columns of array */
113                    j_Equa=Nodes_Equa[k_Equa];
114                    for (l_row=0;l_row<num_subblocks_Equa;++l_row) {
115                       i_row=j_Equa*num_subblocks_Equa+l_row;
116                       /* only look at the matrix rows stored on this processor */
117                       if (i_row < numMyRows) {
118                          for (k_Sol=0;k_Sol<NN_Sol;++k_Sol) { /* Across rows of array */
119                            j_Sol=Nodes_Sol[k_Sol];
120                            for (l_col=0;l_col<num_subblocks_Sol;++l_col) {
121                               /* only look at the matrix rows stored on this processor */
122                               i_col=j_Sol*num_subblocks_Sol+index_offset+l_col;
123                               if (i_col < numMyCols + index_offset ) {
124                               for (k=mainBlock_ptr[i_row]-index_offset;k<mainBlock_ptr[i_row+1]-index_offset;++k) {
125                                   if (mainBlock_index[k]==i_col) {
126                                         /* Entry array(k_Sol, j_Equa) is a block (row_block_size x col_block_size) */
127                                         for (ic=0;ic<col_block_size;++ic) {
128                                               i_Sol=ic+col_block_size*l_col;
129                                               for (ir=0;ir<row_block_size;++ir) {
130                                                  i_Equa=ir+row_block_size*l_row;
131                                          mainBlock_val[k*block_size+ir+row_block_size*ic]+=
132                                                          array[INDEX4(i_Equa,i_Sol,k_Equa,k_Sol,num_Equa,num_Sol,NN_Equa)];
133                                               }
134                                         }
135                                         break;
136                                       }
137                                   }
138                               } else {
139                               for (k=col_coupleBlock_ptr[i_row]-index_offset;k<col_coupleBlock_ptr[i_row+1]-index_offset;++k) {
140                                   if (col_coupleBlock_index[k] == i_col-numMyCols) {
141                                         /* Entry array(k_Sol, j_Equa) is a block (row_block_size x col_block_size) */
142                                         for (ic=0;ic<col_block_size;++ic) {
143                                               i_Sol=ic+col_block_size*l_col;
144                                               for (ir=0;ir<row_block_size;++ir) {
145                                                  i_Equa=ir+row_block_size*l_row;
146                                          col_coupleBlock_val[k*block_size+ir+row_block_size*ic]+=
147                                                          array[INDEX4(i_Equa,i_Sol,k_Equa,k_Sol,num_Equa,num_Sol,NN_Equa)];
148                                               }
149                                         }
150                                         break;
151                                       }
152                                   }
153                               }
154                            }
155                          }
156                        } else {
157                          for (k_Sol=0;k_Sol<NN_Sol;++k_Sol) { /* Across rows of array */
158                            j_Sol=Nodes_Sol[k_Sol];
159                            for (l_col=0;l_col<num_subblocks_Sol;++l_col) {
160                               i_col=j_Sol*num_subblocks_Sol+index_offset+l_col;
161                               if (i_col < numMyCols + index_offset ) {
162                               for (k=row_coupleBlock_ptr[i_row-numMyRows]-index_offset;k<row_coupleBlock_ptr[i_row-numMyRows+1]-index_offset;++k) {
163                                   if (row_coupleBlock_index[k] == i_col) {
164                                         /* Entry array(k_Sol, j_Equa) is a block (row_block_size x col_block_size) */
165                                         for (ic=0;ic<col_block_size;++ic) {
166                                               i_Sol=ic+col_block_size*l_col;
167                                               for (ir=0;ir<row_block_size;++ir) {
168                                                  i_Equa=ir+row_block_size*l_row;
169                                          row_coupleBlock_val[k*block_size+ir+row_block_size*ic]+=
170                                                          array[INDEX4(i_Equa,i_Sol,k_Equa,k_Sol,num_Equa,num_Sol,NN_Equa)];
171                                               }
172                                         }
173                                         break;
174                                       }
175                                   }
176                               }
177                          }                          }
178                          }
179                      }                      }
180                   }                }
181                 }          }
}
}
182     }     }
183  }  }
/*
* \$Log\$
* Revision 1.2  2005/09/15 03:44:21  jgs
* Merge of development branch dev-02 back to main trunk on 2005-09-15
*
* Revision 1.1.2.1  2005/09/07 09:47:29  gross
* missing file
*
* Revision 1.6  2005/07/08 04:07:58  jgs
* Merge of development branch back to main trunk on 2005-07-08
*
* Revision 1.1.1.1.2.3  2005/06/29 02:34:56  gross
* some changes towards 64 integers in finley
*
* Revision 1.1.1.1.2.2  2005/03/15 07:23:55  gross
* 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.
*
* Revision 1.1.1.1.2.1  2004/11/12 06:58:19  gross
* 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
*
* Revision 1.1.1.1  2004/10/26 06:53:57  jgs
* initial import of project esys2
*
* Revision 1.1  2004/07/02 04:21:13  gross
* Finley C code has been included
*
*
*/

Legend:
 Removed from v.616 changed lines Added in v.1811