revision 82 by jgs, Tue Oct 26 06:53:54 2004 UTC revision 97 by jgs, Tue Dec 14 05:39:33 2004 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,int NN_Equa,maybelong* Nodes_Equa, int num_Equa,
29                                                        int NN_col,maybelong* Nodes_col,int col_block_size, double* array) {                                                        int NN_Sol,maybelong* Nodes_Sol, int num_Sol, double* array) {
30    int kr,jr,jc,kc,ir,ic;    int 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;    int row_block_size=in->row_block_size;
32    switch(in->type) {    int col_block_size=in->col_block_size;
33       case CSR:    int block_size=in->block_size;
34            for (kr=0;kr<NN_row;kr++) {    int num_subblocks_Equa=num_Equa/row_block_size;
35              jr=Nodes_row[kr];    int 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+num_subblocks_Sol*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+num_subblocks_Equa*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.1  2004/10/26 06:53:57  jgs   * Revision 1.2  2004/12/14 05:39:30  jgs
94   * Initial revision   * *** empty log message ***
95     *
96     * Revision 1.1.1.1.2.1  2004/11/12 06:58:19  gross
97     * 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
98     *
99     * Revision 1.1.1.1  2004/10/26 06:53:57  jgs
100     * initial import of project esys2
101   *   *
102   * Revision 1.1  2004/07/02 04:21:13  gross   * Revision 1.1  2004/07/02 04:21:13  gross
103   * Finley C code has been included   * Finley C code has been included

Legend:
 Removed from v.82 changed lines Added in v.97