/[escript]/trunk/finley/src/Assemble_addToSystemMatrix.c
ViewVC logotype

Diff of /trunk/finley/src/Assemble_addToSystemMatrix.c

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 616 by elspeth, Wed Mar 22 02:46:56 2006 UTC revision 1552 by gross, Thu May 8 08:52:41 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
8   *     http://www.opensource.org/licenses/osl-3.0.php       *   *
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 */  /* 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, *col_coupleBlock_ptr, *col_coupleBlock_index, *row_coupleBlock_ptr, *row_coupleBlock_index;
41      double *mainBlock_val, *row_coupleBlock_val, *col_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                            }                            }
# Line 69  void  Finley_Assemble_addToSystemMatrix( Line 78  void  Finley_Assemble_addToSystemMatrix(
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]+=              }
99                                             array[INDEX4(i_Equa,i_Sol,k_Equa,k_Sol,num_Equa,num_Sol,NN_Equa)];            }
100                                  }         #endif
101                            }     } else {
102                            break;  
103             mainBlock_ptr=in->mainBlock->pattern->ptr;
104             mainBlock_index=in->mainBlock->pattern->index;
105             mainBlock_val=in->mainBlock->val;
106             col_coupleBlock_ptr=in->col_coupleBlock->pattern->ptr;
107             col_coupleBlock_index=in->col_coupleBlock->pattern->index;
108             col_coupleBlock_val=in->col_coupleBlock->val;
109             row_coupleBlock_ptr=in->row_coupleBlock->pattern->ptr;
110             row_coupleBlock_index=in->row_coupleBlock->pattern->index;
111             row_coupleBlock_val=in->row_coupleBlock->val;
112    
113             for (k_Equa=0;k_Equa<NN_Equa;++k_Equa) { /* Down columns of array */
114                    j_Equa=Nodes_Equa[k_Equa];
115                    for (l_row=0;l_row<num_subblocks_Equa;++l_row) {
116                       i_row=j_Equa*num_subblocks_Equa+l_row;
117                       /* only look at the matrix rows stored on this processor */
118                       if (i_row < numMyRows) {
119                          for (k_Sol=0;k_Sol<NN_Sol;++k_Sol) { /* Across rows of array */
120                            j_Sol=Nodes_Sol[k_Sol];
121                            for (l_col=0;l_col<num_subblocks_Sol;++l_col) {
122                               /* only look at the matrix rows stored on this processor */
123                               i_col=j_Sol*num_subblocks_Sol+index_offset+l_col;
124                               if (i_col < numMyCols + index_offset ) {
125                               for (k=mainBlock_ptr[i_row]-index_offset;k<mainBlock_ptr[i_row+1]-index_offset;++k) {
126                                   if (mainBlock_index[k]==i_col) {
127                                         /* Entry array(k_Sol, j_Equa) is a block (row_block_size x col_block_size) */
128                                         for (ic=0;ic<col_block_size;++ic) {
129                                               i_Sol=ic+col_block_size*l_col;
130                                               for (ir=0;ir<row_block_size;++ir) {
131                                                  i_Equa=ir+row_block_size*l_row;
132                                          mainBlock_val[k*block_size+ir+row_block_size*ic]+=
133                                                          array[INDEX4(i_Equa,i_Sol,k_Equa,k_Sol,num_Equa,num_Sol,NN_Equa)];
134                                               }
135                                         }
136                                         break;
137                                       }
138                                   }
139                               } else {
140                               for (k=col_coupleBlock_ptr[i_row]-index_offset;k<col_coupleBlock_ptr[i_row+1]-index_offset;++k) {
141                                   if (col_coupleBlock_index[k] == i_col-numMyCols) {
142                                         /* Entry array(k_Sol, j_Equa) is a block (row_block_size x col_block_size) */
143                                         for (ic=0;ic<col_block_size;++ic) {
144                                               i_Sol=ic+col_block_size*l_col;
145                                               for (ir=0;ir<row_block_size;++ir) {
146                                                  i_Equa=ir+row_block_size*l_row;
147                                          col_coupleBlock_val[k*block_size+ir+row_block_size*ic]+=
148                                                          array[INDEX4(i_Equa,i_Sol,k_Equa,k_Sol,num_Equa,num_Sol,NN_Equa)];
149                                               }
150                                         }
151                                         break;
152                                       }
153                                   }
154                               }
155                          }                          }
156                          }
157                        } else {
158                          for (k_Sol=0;k_Sol<NN_Sol;++k_Sol) { /* Across rows of array */
159                            j_Sol=Nodes_Sol[k_Sol];
160                            for (l_col=0;l_col<num_subblocks_Sol;++l_col) {
161                               i_col=j_Sol*num_subblocks_Sol+index_offset+l_col;
162                               if (i_col < numMyCols + index_offset ) {
163                               for (k=row_coupleBlock_ptr[i_row-numMyRows]-index_offset;k<row_coupleBlock_ptr[i_row-numMyRows+1]-index_offset;++k) {
164                                   if (row_coupleBlock_index[k] == i_col) {
165                                         /* Entry array(k_Sol, j_Equa) is a block (row_block_size x col_block_size) */
166                                         for (ic=0;ic<col_block_size;++ic) {
167                                               i_Sol=ic+col_block_size*l_col;
168                                               for (ir=0;ir<row_block_size;++ir) {
169                                                  i_Equa=ir+row_block_size*l_row;
170                                          row_coupleBlock_val[k*block_size+ir+row_block_size*ic]+=
171                                                          array[INDEX4(i_Equa,i_Sol,k_Equa,k_Sol,num_Equa,num_Sol,NN_Equa)];
172                                               }
173                                         }
174                                         break;
175                                       }
176                                   }
177                               }
178                            }
179                          }
180                      }                      }
181                   }                }
182                 }          }
              }  
           }  
183     }     }
184  }  }
 /*  
  * $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.1552

  ViewVC Help
Powered by ViewVC 1.1.26