/[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

trunk/finley/src/finleyC/Assemble_addToSystemMatrix.c revision 155 by jgs, Wed Nov 9 02:02:19 2005 UTC trunk/finley/src/Assemble_addToSystemMatrix.c revision 1322 by ksteube, Thu Sep 27 04:40:43 2007 UTC
# Line 1  Line 1 
1  /*  
2   ******************************************************************************  /* $Id$ */
3   *                                                                            *  
4   *       COPYRIGHT  ACcESS 2003,2004,2005 -  All Rights Reserved              *  /*******************************************************
5   *                                                                            *   *
6   * This software is the property of ACcESS. No part of this code              *   *           Copyright 2003-2007 by ACceSS MNRF
7   * may be copied in any form or by any means without the expressed written    *   *       Copyright 2007 by University of Queensland
8   * consent of ACcESS.  Copying, use or modification of this software          *   *
9   * by any unauthorised person is illegal unless that person has a software    *   *                http://esscc.uq.edu.au
10   * license agreement with ACcESS.                                             *   *        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 26  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    dim_t k_Equa,j_Equa,j_Sol,k_Sol,i_Equa,i_Sol,l_col,l_row,ic,ir,index,k,iptr;    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,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    if (in->type==CSR) {  
51            for (k_Equa=0;k_Equa<NN_Equa;++k_Equa) {    if (in->type & MATRIX_FORMAT_CSC) {
52              j_Equa=Nodes_Equa[k_Equa];           /* MATRIX_FORMAT_CSC does not support MPI !!!!! */
53              for (l_row=0;l_row<num_subblocks_Equa;++l_row) {           mainBlock_ptr=in->mainBlock->pattern->ptr;
54                 iptr=j_Equa*num_subblocks_Equa+l_row;           mainBlock_index=in->mainBlock->pattern->index;
55                 for (k_Sol=0;k_Sol<NN_Sol;++k_Sol) {           mainBlock_val=in->mainBlock->val;
                  j_Sol=Nodes_Sol[k_Sol];  
                  for (l_col=0;l_col<num_subblocks_Sol;++l_col) {  
                     index=j_Sol*num_subblocks_Sol+INDEX_OFFSET+l_col;  
                 for (k=in->pattern->ptr[iptr]-PTR_OFFSET;k<in->pattern->ptr[iptr+1]-PTR_OFFSET;++k) {  
                     if (in->pattern->index[k]==index) {  
                           for (ic=0;ic<col_block_size;++ic) {  
                                 i_Sol=ic+col_block_size*l_col;  
                                 for (ir=0;ir<row_block_size;++ir) {  
                                    i_Equa=ir+row_block_size*l_row;  
                            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;  
                         }  
                     }  
                  }  
                }  
              }  
           }  
      } else {  
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]-PTR_OFFSET;k<in->pattern->ptr[iptr+1]-PTR_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 96  void  Finley_Assemble_addToSystemMatrix( Line 78  void  Finley_Assemble_addToSystemMatrix(
78                 }                 }
79              }              }
80           }           }
81       } else if (in->type & MATRIX_FORMAT_TRILINOS_CRS) {
82           /* 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];
86            if (j_Equa < in->mainBlock->pattern->output_node_distribution->numLocal) {
87                  for (k_Sol=0;k_Sol<NN_Sol;++k_Sol) { /* Across rows of array */
88                    j_Sol=Nodes_Sol[k_Sol];
89                    for (l_row=0;l_row<num_subblocks_Equa;++l_row) {
90                      irow=j_Equa*row_block_size+l_row;
91                      for (l_col=0;l_col<col_block_size;++l_col) {
92                         icol=j_Sol*col_block_size+index_offset+l_col;
93                 /* irow is local and icol is global */
94                 Trilinos_SumIntoMyValues(in->trilinos_data, irow, icol, array[INDEX4(l_row,l_col,k_Equa,k_Sol,num_Equa,num_Sol,NN_Equa)]);
95                  }
96                    }
97                  }
98                }
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.155  
changed lines
  Added in v.1322

  ViewVC Help
Powered by ViewVC 1.1.26