revision 616 by elspeth, Wed Mar 22 02:46:56 2006 UTC revision 2548 by jfenwick, Mon Jul 20 06:20:06 2009 UTC
# Line 1  Line 1
1  /*
2   ************************************************************  /*******************************************************
3   *          Copyright 2006 by ACcESS MNRF                   *  *
4   *                                                          *  * Copyright (c) 2003-2009 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    if (in->type & MATRIX_FORMAT_CSC) {    if (in->type & MATRIX_FORMAT_CSC) {
50             /* MATRIX_FORMAT_CSC does not support MPI !!!!! */
51             mainBlock_ptr=in->mainBlock->pattern->ptr;
52             mainBlock_index=in->mainBlock->pattern->index;
53             mainBlock_val=in->mainBlock->val;
54           for (k_Sol=0;k_Sol<NN_Sol;k_Sol++) {           for (k_Sol=0;k_Sol<NN_Sol;k_Sol++) {
55              j_Sol=Nodes_Sol[k_Sol];              j_Sol=Nodes_Sol[k_Sol];
56              for (l_col=0;l_col<num_subblocks_Sol;++l_col) {              for (l_col=0;l_col<num_subblocks_Sol;++l_col) {
57                 iptr=j_Sol*num_subblocks_Sol+l_col;                 i_col=j_Sol*num_subblocks_Sol+l_col;
58                 for (k_Equa=0;k_Equa<NN_Equa;k_Equa++) {                 for (k_Equa=0;k_Equa<NN_Equa;k_Equa++) {
59                   j_Equa=Nodes_Equa[k_Equa];                   j_Equa=Nodes_Equa[k_Equa];
60                   for (l_row=0;l_row<num_subblocks_Equa;++l_row) {                   for (l_row=0;l_row<num_subblocks_Equa;++l_row) {
61                      index=j_Equa*num_subblocks_Equa+index_offset+l_row;                      i_row=j_Equa*num_subblocks_Equa+index_offset+l_row;
62                  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) {
63                      if (in->pattern->index[k]==index) {                      if (mainBlock_index[k] == i_row) {
64                            for (ic=0;ic<col_block_size;++ic) {                            for (ic=0;ic<col_block_size;++ic) {
65                                  i_Sol=ic+col_block_size*l_col;                                  i_Sol=ic+col_block_size*l_col;
66                                  for (ir=0;ir<row_block_size;++ir) {                                  for (ir=0;ir<row_block_size;++ir) {
67                                     i_Equa=ir+row_block_size*l_row;                                     i_Equa=ir+row_block_size*l_row;
68                             in->val[k*block_size+ir+row_block_size*ic]+=                             mainBlock_val[k*block_size+ir+row_block_size*ic]+=
69                                             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)];
70                                  }                                  }
71                            }                            }
76                 }                 }
77              }              }
78           }           }
79     } else {     } else if (in->type & MATRIX_FORMAT_TRILINOS_CRS) {
80            for (k_Equa=0;k_Equa<NN_Equa;++k_Equa) {         /* this needs to be modified */
81           #ifdef TRILINOS
82              for (k_Equa=0;k_Equa<NN_Equa;++k_Equa) { /* Down columns of array */
83              j_Equa=Nodes_Equa[k_Equa];              j_Equa=Nodes_Equa[k_Equa];
84              for (l_row=0;l_row<num_subblocks_Equa;++l_row) {          if (j_Equa < in->mainBlock->pattern->output_node_distribution->numLocal) {
85                 iptr=j_Equa*num_subblocks_Equa+l_row;                for (k_Sol=0;k_Sol<NN_Sol;++k_Sol) { /* Across rows of array */
86                 for (k_Sol=0;k_Sol<NN_Sol;++k_Sol) {                  j_Sol=Nodes_Sol[k_Sol];
87                   j_Sol=Nodes_Sol[k_Sol];                  for (l_row=0;l_row<num_subblocks_Equa;++l_row) {
88                   for (l_col=0;l_col<num_subblocks_Sol;++l_col) {                    irow=j_Equa*row_block_size+l_row;
89                      index=j_Sol*num_subblocks_Sol+index_offset+l_col;                    for (l_col=0;l_col<col_block_size;++l_col) {
90                  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;
91                      if (in->pattern->index[k]==index) {               /* irow is local and icol is global */
92                            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)]);
93                                  i_Sol=ic+col_block_size*l_col;                }
94                                  for (ir=0;ir<row_block_size;++ir) {                  }
95                                     i_Equa=ir+row_block_size*l_row;                }
96                             in->val[k*block_size+ir+row_block_size*ic]+=              }
97                                             array[INDEX4(i_Equa,i_Sol,k_Equa,k_Sol,num_Equa,num_Sol,NN_Equa)];            }
98                                  }         #endif
99                            }     } else {
100                            break;           mainBlock_ptr=in->mainBlock->pattern->ptr;
101             mainBlock_index=in->mainBlock->pattern->index;
102             mainBlock_val=in->mainBlock->val;
103             col_coupleBlock_ptr=in->col_coupleBlock->pattern->ptr;
104             col_coupleBlock_index=in->col_coupleBlock->pattern->index;
105             col_coupleBlock_val=in->col_coupleBlock->val;
106             row_coupleBlock_ptr=in->row_coupleBlock->pattern->ptr;
107             row_coupleBlock_index=in->row_coupleBlock->pattern->index;
108             row_coupleBlock_val=in->row_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=col_coupleBlock_ptr[i_row]-index_offset;k<col_coupleBlock_ptr[i_row+1]-index_offset;++k) {
138                                   if (col_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                                          col_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                        } else {
155                          for (k_Sol=0;k_Sol<NN_Sol;++k_Sol) { /* Across rows of array */
156                            j_Sol=Nodes_Sol[k_Sol];
157                            for (l_col=0;l_col<num_subblocks_Sol;++l_col) {
158                               i_col=j_Sol*num_subblocks_Sol+index_offset+l_col;
159                               if (i_col < numMyCols + index_offset ) {
160                               for (k=row_coupleBlock_ptr[i_row-numMyRows]-index_offset;k<row_coupleBlock_ptr[i_row-numMyRows+1]-index_offset;++k) {
161                                   if (row_coupleBlock_index[k] == i_col) {
162                                         /* Entry array(k_Sol, j_Equa) is a block (row_block_size x col_block_size) */
163                                         for (ic=0;ic<col_block_size;++ic) {
164                                               i_Sol=ic+col_block_size*l_col;
165                                               for (ir=0;ir<row_block_size;++ir) {
166                                                  i_Equa=ir+row_block_size*l_row;
167                                          row_coupleBlock_val[k*block_size+ir+row_block_size*ic]+=
168                                                          array[INDEX4(i_Equa,i_Sol,k_Equa,k_Sol,num_Equa,num_Sol,NN_Equa)];
169                                               }
170                                         }
171                                         break;
172                                       }
173                                   }
174                               }
175                          }                          }
176                          }
177                      }                      }
178                   }                }
179                 }          }
}
}
180     }     }
181  }  }
/*
* \$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.2548