/[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 100 by jgs, Wed Dec 15 03:48:48 2004 UTC revision 123 by jgs, Fri Jul 8 04:08:13 2005 UTC
# Line 4  Line 4 
4    
5  /* Finley: SystemMatrix and SystemVector */  /* Finley: SystemMatrix and SystemVector */
6    
7  /*  adds the matrix array[row_block_size,col_block_size,NN,NN] onto the matrix in. */  /*  adds the matrix array[Equa,Sol,NN,NN] onto the matrix in. */
8  /* the rows/columns are given by */  /* the rows/columns are given by */
9  /*  ir+row_block_size*Nodes_row[Nodes[jr]] (ir=0:row_block_size; jr=0:NN_row). */  /*  i_Equa+Equa*Nodes_Equa[Nodes[j_Equa]] (i_Equa=0:Equa; j_Equa=0:NN_Equa). */
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->row_block_size=in->col_block_size=1, i.e. */  /*  This routine assumes that in->Equa=in->Sol=1, i.e. */
13  /*  array is fully packed. */  /*  array is fully packed. */
14  /* TODO: the case in->row_block_size!=1  */  /* TODO: the case in->Equa!=1  */
15    
16  /**************************************************************/  /**************************************************************/
17    
# Line 25  Line 25 
25    
26  /**************************************************************/  /**************************************************************/
27    
28  void  Finley_SystemMatrix_add(Finley_SystemMatrix* in,int NN_row,maybelong* Nodes_row, int row_block_size,  void  Finley_SystemMatrix_add(Finley_SystemMatrix* in,dim_t NN_Equa,index_t* Nodes_Equa, dim_t num_Equa,
29                                                        int NN_col,maybelong* Nodes_col,int col_block_size, double* array) {                                                        dim_t NN_Sol,index_t* Nodes_Sol, dim_t num_Sol, double* array) {
30    int kr,jr,jc,kc,ir,ic;    dim_t k_Equa,j_Equa,j_Sol,k_Sol,i_Equa,i_Sol,l_col,l_row,ic,ir,index,k,iptr;
31    maybelong irow,icol,k;    dim_t row_block_size=in->row_block_size;
32    switch(in->type) {    dim_t col_block_size=in->col_block_size;
33       case CSR:    dim_t block_size=in->block_size;
34            for (kr=0;kr<NN_row;kr++) {    dim_t num_subblocks_Equa=num_Equa/row_block_size;
35              jr=Nodes_row[kr];    dim_t num_subblocks_Sol=num_Sol/col_block_size;
36              for (ir=0;ir<row_block_size;ir++) {  
37                 irow=ir+row_block_size*jr;    if (in->type==CSR) {
38                 for (kc=0;kc<NN_col;kc++) {            for (k_Equa=0;k_Equa<NN_Equa;++k_Equa) {
39               jc=Nodes_col[kc];              j_Equa=Nodes_Equa[k_Equa];
40               for (ic=0;ic<col_block_size;ic++) {              for (l_row=0;l_row<num_subblocks_Equa;++l_row) {
41                  icol=ic+col_block_size*jc+INDEX_OFFSET;                 iptr=j_Equa*num_subblocks_Equa+l_row;
42                  for(k=in->ptr[irow]-PTR_OFFSET;k<in->ptr[irow+1]-PTR_OFFSET;k++) {                 for (k_Sol=0;k_Sol<NN_Sol;++k_Sol) {
43                     if (in->index[k]==icol) {                   j_Sol=Nodes_Sol[k_Sol];
44                            #pragma omp atomic                   for (l_col=0;l_col<num_subblocks_Sol;++l_col) {
45                    in->val[k]+=array[INDEX4(ir,ic,kr,kc,row_block_size,col_block_size,NN_row)];                      index=j_Sol*num_subblocks_Sol+INDEX_OFFSET+l_col;
46                    for (k=in->pattern->ptr[iptr]-PTR_OFFSET;k<in->pattern->ptr[iptr+1]-PTR_OFFSET;++k) {
47                        if (in->pattern->index[k]==index) {
48                              for (ic=0;ic<col_block_size;++ic) {
49                                    i_Sol=ic+col_block_size*l_col;
50                                    for (ir=0;ir<row_block_size;++ir) {
51                                       i_Equa=ir+row_block_size*l_row;
52                               in->val[k*block_size+ir+row_block_size*ic]+=
53                                               array[INDEX4(i_Equa,i_Sol,k_Equa,k_Sol,num_Equa,num_Sol,NN_Equa)];
54                                    }
55                              }
56                            break;                            break;
57                         }                          }
58                  }                      }
59               }                   }
60                }                 }
61                 }
62            }            }
63        }       } else {
64        break;           for (k_Sol=0;k_Sol<NN_Sol;k_Sol++) {
65     case CSC:              j_Sol=Nodes_Sol[k_Sol];
66          for (kr=0;kr<NN_row;kr++) {              for (l_col=0;l_col<num_subblocks_Sol;++l_col) {
67            jr=Nodes_row[kr];                 iptr=j_Sol*num_subblocks_Sol+l_col;
68            for (ir=0;ir<row_block_size;ir++) {                 for (k_Equa=0;k_Equa<NN_Equa;k_Equa++) {
69               irow=ir+row_block_size*jr;                   j_Equa=Nodes_Equa[k_Equa];
70               for (kc=0;kc<NN_col;kc++) {                   for (l_row=0;l_row<num_subblocks_Equa;++l_row) {
71              jc=Nodes_col[kc];                      index=j_Equa*num_subblocks_Equa+INDEX_OFFSET+l_row;
72              for (ic=0;ic<col_block_size;ic++) {                  for (k=in->pattern->ptr[iptr]-PTR_OFFSET;k<in->pattern->ptr[iptr+1]-PTR_OFFSET;++k) {
73                 icol=ic+col_block_size*jc+INDEX_OFFSET;                      if (in->pattern->index[k]==index) {
74                 for(k=in->ptr[icol]-PTR_OFFSET;k<in->ptr[icol+1]-PTR_OFFSET;k++) {                            for (ic=0;ic<col_block_size;++ic) {
75                    if (in->index[k]==irow) {                                  i_Sol=ic+col_block_size*l_col;
76                           #pragma omp atomic                                  for (ir=0;ir<row_block_size;++ir) {
77                   in->val[k]+=array[INDEX4(ir,ic,kr,kc,row_block_size,col_block_size,NN_row)];                                     i_Equa=ir+row_block_size*l_row;
78                           break;                             in->val[k*block_size+ir+row_block_size*ic]+=
79                        }                                             array[INDEX4(i_Equa,i_Sol,k_Equa,k_Sol,num_Equa,num_Sol,NN_Equa)];
80                     }                                  }
81                  }                            }
82           }                            break;
83        }                          }
84          }                      }
85      break;                   }
86     default:                 }
87         Finley_ErrorCode=TYPE_ERROR;              }
88         sprintf(Finley_ErrorMsg,"Unknown matrix type.");           }
89    } /* switch in->type */     }
90  }  }
91  /*  /*
92   * $Log$   * $Log$
93   * Revision 1.3  2004/12/15 03:48:46  jgs   * Revision 1.6  2005/07/08 04:07:58  jgs
94   * *** empty log message ***   * Merge of development branch back to main trunk on 2005-07-08
95     *
96     * Revision 1.1.1.1.2.3  2005/06/29 02:34:56  gross
97     * some changes towards 64 integers in finley
98     *
99     * Revision 1.1.1.1.2.2  2005/03/15 07:23:55  gross
100     * 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.
101     *
102     * Revision 1.1.1.1.2.1  2004/11/12 06:58:19  gross
103     * 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
104   *   *
105   * Revision 1.1.1.1  2004/10/26 06:53:57  jgs   * Revision 1.1.1.1  2004/10/26 06:53:57  jgs
106   * initial import of project esys2   * initial import of project esys2

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

  ViewVC Help
Powered by ViewVC 1.1.26