/[escript]/trunk/esys2/finley/src/finleyC/System_add.c
ViewVC logotype

Diff of /trunk/esys2/finley/src/finleyC/System_add.c

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

revision 99 by jgs, Tue Dec 14 05:39:33 2004 UTC revision 100 by jgs, Wed Dec 15 03:48:48 2004 UTC
# Line 4  Line 4 
4    
5  /* Finley: SystemMatrix and SystemVector */  /* Finley: SystemMatrix and SystemVector */
6    
7  /*  adds the matrix array[Equa,Sol,NN,NN] onto the matrix in. */  /*  adds the matrix array[row_block_size,col_block_size,NN,NN] onto the matrix in. */
8  /* the rows/columns are given by */  /* the rows/columns are given by */
9  /*  i_Equa+Equa*Nodes_Equa[Nodes[j_Equa]] (i_Equa=0:Equa; j_Equa=0:NN_Equa). */  /*  ir+row_block_size*Nodes_row[Nodes[jr]] (ir=0:row_block_size; jr=0:NN_row). */
10  /*  the routine has to be called from a parallel region                        */  /*  the routine has to be called from a parallel region                        */
11    
12  /*  This routine assumes that in->Equa=in->Sol=1, i.e. */  /*  This routine assumes that in->row_block_size=in->col_block_size=1, i.e. */
13  /*  array is fully packed. */  /*  array is fully packed. */
14  /* TODO: the case in->Equa!=1  */  /* TODO: the case in->row_block_size!=1  */
15    
16  /**************************************************************/  /**************************************************************/
17    
# Line 25  Line 25 
25    
26  /**************************************************************/  /**************************************************************/
27    
28  void  Finley_SystemMatrix_add(Finley_SystemMatrix* in,int NN_Equa,maybelong* Nodes_Equa, int num_Equa,  void  Finley_SystemMatrix_add(Finley_SystemMatrix* in,int NN_row,maybelong* Nodes_row, int row_block_size,
29                                                        int NN_Sol,maybelong* Nodes_Sol, int num_Sol, double* array) {                                                        int NN_col,maybelong* Nodes_col,int col_block_size, double* array) {
30    int k_Equa,j_Equa,j_Sol,k_Sol,i_Equa,i_Sol,l_col,l_row,ic,ir,index,k,iptr;    int kr,jr,jc,kc,ir,ic;
31    int row_block_size=in->row_block_size;    maybelong irow,icol,k;
32    int col_block_size=in->col_block_size;    switch(in->type) {
33    int block_size=in->block_size;       case CSR:
34    int num_subblocks_Equa=num_Equa/row_block_size;            for (kr=0;kr<NN_row;kr++) {
35    int num_subblocks_Sol=num_Sol/col_block_size;              jr=Nodes_row[kr];
36                for (ir=0;ir<row_block_size;ir++) {
37    if (in->type==CSR) {                 irow=ir+row_block_size*jr;
38            for (k_Equa=0;k_Equa<NN_Equa;++k_Equa) {                 for (kc=0;kc<NN_col;kc++) {
39              j_Equa=Nodes_Equa[k_Equa];               jc=Nodes_col[kc];
40              for (l_row=0;l_row<num_subblocks_Equa;++l_row) {               for (ic=0;ic<col_block_size;ic++) {
41                 iptr=j_Equa*num_subblocks_Equa+l_row;                  icol=ic+col_block_size*jc+INDEX_OFFSET;
42                 for (k_Sol=0;k_Sol<NN_Sol;++k_Sol) {                  for(k=in->ptr[irow]-PTR_OFFSET;k<in->ptr[irow+1]-PTR_OFFSET;k++) {
43                   j_Sol=Nodes_Sol[k_Sol];                     if (in->index[k]==icol) {
44                   for (l_col=0;l_col<num_subblocks_Sol;++l_col) {                            #pragma omp atomic
45                      index=j_Sol*num_subblocks_Sol+INDEX_OFFSET+l_col;                    in->val[k]+=array[INDEX4(ir,ic,kr,kc,row_block_size,col_block_size,NN_row)];
                 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)];  
                                 }  
                           }  
46                            break;                            break;
47                          }                         }
48                      }                  }
49                   }               }
50                 }                }
              }  
51            }            }
52       } else {        }
53           for (k_Sol=0;k_Sol<NN_Sol;k_Sol++) {        break;
54              j_Sol=Nodes_Sol[k_Sol];     case CSC:
55              for (l_col=0;l_col<num_subblocks_Sol;++l_col) {          for (kr=0;kr<NN_row;kr++) {
56                 iptr=j_Sol*num_subblocks_Sol+l_col;            jr=Nodes_row[kr];
57                 for (k_Equa=0;k_Equa<NN_Equa;k_Equa++) {            for (ir=0;ir<row_block_size;ir++) {
58                   j_Equa=Nodes_Equa[k_Equa];               irow=ir+row_block_size*jr;
59                   for (l_row=0;l_row<num_subblocks_Equa;++l_row) {               for (kc=0;kc<NN_col;kc++) {
60                      index=j_Equa*num_subblocks_Equa+INDEX_OFFSET+l_row;              jc=Nodes_col[kc];
61                  for (k=in->pattern->ptr[iptr]-PTR_OFFSET;k<in->pattern->ptr[iptr+1]-PTR_OFFSET;++k) {              for (ic=0;ic<col_block_size;ic++) {
62                      if (in->pattern->index[k]==index) {                 icol=ic+col_block_size*jc+INDEX_OFFSET;
63                            for (ic=0;ic<col_block_size;++ic) {                 for(k=in->ptr[icol]-PTR_OFFSET;k<in->ptr[icol+1]-PTR_OFFSET;k++) {
64                                  i_Sol=ic+num_subblocks_Sol*l_col;                    if (in->index[k]==irow) {
65                                  for (ir=0;ir<row_block_size;++ir) {                           #pragma omp atomic
66                                     i_Equa=ir+num_subblocks_Equa*l_row;                   in->val[k]+=array[INDEX4(ir,ic,kr,kc,row_block_size,col_block_size,NN_row)];
67                             in->val[k*block_size+ir+row_block_size*ic]+=                           break;
68                                             array[INDEX4(i_Equa,i_Sol,k_Equa,k_Sol,num_Equa,num_Sol,NN_Equa)];                        }
69                                  }                     }
70                            }                  }
71                            break;           }
72                          }        }
73                      }          }
74                   }      break;
75                 }     default:
76              }         Finley_ErrorCode=TYPE_ERROR;
77           }         sprintf(Finley_ErrorMsg,"Unknown matrix type.");
78     }    } /* switch in->type */
79  }  }
80  /*  /*
81   * $Log$   * $Log$
82   * Revision 1.2  2004/12/14 05:39:30  jgs   * Revision 1.3  2004/12/15 03:48:46  jgs
83   * *** empty log message ***   * *** empty log message ***
84   *   *
  * 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  
  *  
85   * Revision 1.1.1.1  2004/10/26 06:53:57  jgs   * Revision 1.1.1.1  2004/10/26 06:53:57  jgs
86   * initial import of project esys2   * initial import of project esys2
87   *   *

Legend:
Removed from v.99  
changed lines
  Added in v.100

  ViewVC Help
Powered by ViewVC 1.1.26