revision 616 by elspeth, Wed Mar 22 02:46:56 2006 UTC revision 1388 by trankine, Fri Jan 11 07:45:58 2008 UTC
# Line 1  Line 1
1  /*
2   ************************************************************  /* \$Id\$ */
3   *          Copyright 2006 by ACcESS MNRF                   *
4   *                                                          *  /*******************************************************
5   *              http://www.access.edu.au                    *   *
6   *       Primary Business: Queensland, Australia            *   *           Copyright 2003-2007 by ACceSS MNRF
7   *  Licensed under the Open Software License version 3.0    *   *       Copyright 2007 by University of Queensland
9   *                                                          *   *                http://esscc.uq.edu.au
10   ************************************************************   *        Primary Business: Queensland, Australia
13     *
14     *******************************************************/
15
16  /**************************************************************/  /**************************************************************/
17
18  /* Finley: SystemMatrix and SystemVector */  /* Finley: SystemMatrix and SystemVector */
# Line 24  Line 28
28
29  /**************************************************************/  /**************************************************************/
30
/*  Author: gross@access.edu.au */
/*  Version: \$Id\$ */

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

31  #include "Assemble.h"  #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,  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) {                                          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);    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,iptr;    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;    dim_t row_block_size=in->row_block_size;
43    dim_t col_block_size=in->col_block_size;    dim_t col_block_size=in->col_block_size;
44    dim_t block_size=in->block_size;    dim_t block_size=in->block_size;
45    dim_t num_subblocks_Equa=num_Equa/row_block_size;    dim_t num_subblocks_Equa=num_Equa/row_block_size;
46    dim_t num_subblocks_Sol=num_Sol/col_block_size;    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) {    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++) {           for (k_Sol=0;k_Sol<NN_Sol;k_Sol++) {
57              j_Sol=Nodes_Sol[k_Sol];              j_Sol=Nodes_Sol[k_Sol];
58              for (l_col=0;l_col<num_subblocks_Sol;++l_col) {              for (l_col=0;l_col<num_subblocks_Sol;++l_col) {
59                 iptr=j_Sol*num_subblocks_Sol+l_col;                 i_col=j_Sol*num_subblocks_Sol+l_col;
60                 for (k_Equa=0;k_Equa<NN_Equa;k_Equa++) {                 for (k_Equa=0;k_Equa<NN_Equa;k_Equa++) {
61                   j_Equa=Nodes_Equa[k_Equa];                   j_Equa=Nodes_Equa[k_Equa];
62                   for (l_row=0;l_row<num_subblocks_Equa;++l_row) {                   for (l_row=0;l_row<num_subblocks_Equa;++l_row) {
63                      index=j_Equa*num_subblocks_Equa+index_offset+l_row;                      i_row=j_Equa*num_subblocks_Equa+index_offset+l_row;
64                  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) {
65                      if (in->pattern->index[k]==index) {                      if (mainBlock_index[k] == i_row) {
66                            for (ic=0;ic<col_block_size;++ic) {                            for (ic=0;ic<col_block_size;++ic) {
67                                  i_Sol=ic+col_block_size*l_col;                                  i_Sol=ic+col_block_size*l_col;
68                                  for (ir=0;ir<row_block_size;++ir) {                                  for (ir=0;ir<row_block_size;++ir) {
69                                     i_Equa=ir+row_block_size*l_row;                                     i_Equa=ir+row_block_size*l_row;
70                             in->val[k*block_size+ir+row_block_size*ic]+=                             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)];                                             array[INDEX4(i_Equa,i_Sol,k_Equa,k_Sol,num_Equa,num_Sol,NN_Equa)];
72                                  }                                  }
73                            }                            }
78                 }                 }
79              }              }
80           }           }
81     } else {     } else if (in->type & MATRIX_FORMAT_TRILINOS_CRS) {
82            for (k_Equa=0;k_Equa<NN_Equa;++k_Equa) {         /* 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];              j_Equa=Nodes_Equa[k_Equa];
86              for (l_row=0;l_row<num_subblocks_Equa;++l_row) {          if (j_Equa < in->mainBlock->pattern->output_node_distribution->numLocal) {
87                 iptr=j_Equa*num_subblocks_Equa+l_row;                for (k_Sol=0;k_Sol<NN_Sol;++k_Sol) { /* Across rows of array */
88                 for (k_Sol=0;k_Sol<NN_Sol;++k_Sol) {                  j_Sol=Nodes_Sol[k_Sol];
89                   j_Sol=Nodes_Sol[k_Sol];                  for (l_row=0;l_row<num_subblocks_Equa;++l_row) {
90                   for (l_col=0;l_col<num_subblocks_Sol;++l_col) {                    irow=j_Equa*row_block_size+l_row;
91                      index=j_Sol*num_subblocks_Sol+index_offset+l_col;                    for (l_col=0;l_col<col_block_size;++l_col) {
92                  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;
93                      if (in->pattern->index[k]==index) {               /* irow is local and icol is global */
94                            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)]);
95                                  i_Sol=ic+col_block_size*l_col;                }
96                                  for (ir=0;ir<row_block_size;++ir) {                  }
97                                     i_Equa=ir+row_block_size*l_row;                }
98                             in->val[k*block_size+ir+row_block_size*ic]+=              }
array[INDEX4(i_Equa,i_Sol,k_Equa,k_Sol,num_Equa,num_Sol,NN_Equa)];
}
}
break;
}
}
}
}
}
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  }  }
/*
* \$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.1388