/[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 969 by ksteube, Tue Feb 13 23:02:23 2007 UTC
# Line 43  void  Finley_Assemble_addToSystemMatrix( Line 43  void  Finley_Assemble_addToSystemMatrix(
43    dim_t num_subblocks_Equa=num_Equa/row_block_size;    dim_t num_subblocks_Equa=num_Equa/row_block_size;
44    dim_t num_subblocks_Sol=num_Sol/col_block_size;    dim_t num_subblocks_Sol=num_Sol/col_block_size;
45    
46      printf("ksteube addToSystemMatrix NN_Sol=%d num_subblocks_Sol=%d NN_Equa=%d num_subblocks_Equa=%d\n", NN_Sol, num_subblocks_Sol, NN_Equa, num_subblocks_Equa);
47    if (in->type & MATRIX_FORMAT_CSC) {    if (in->type & MATRIX_FORMAT_CSC) {
48           for (k_Sol=0;k_Sol<NN_Sol;k_Sol++) {           for (k_Sol=0;k_Sol<NN_Sol;k_Sol++) {
49              j_Sol=Nodes_Sol[k_Sol];              j_Sol=Nodes_Sol[k_Sol];
# Line 69  void  Finley_Assemble_addToSystemMatrix( Line 70  void  Finley_Assemble_addToSystemMatrix(
70                 }                 }
71              }              }
72           }           }
73       } else if (in->type & MATRIX_FORMAT_TRILINOS_CRS) {
74    #ifdef PASO_MPI
75              for (k_Equa=0;k_Equa<NN_Equa;++k_Equa) { /* Down columns of array */
76                j_Equa=Nodes_Equa[k_Equa];
77            if (j_Equa < row_degreeOfFreedomDistribution->numLocal) {
78                  for (k_Sol=0;k_Sol<NN_Sol;++k_Sol) { /* Across rows of array */
79                    j_Sol=Nodes_Sol[k_Sol];
80                j_Sol = Finley_IndexList_localToGlobal(col_degreeOfFreedomDistribution, my_CPU, j_Sol);
81                    for (l_row=0;l_row<num_subblocks_Equa;++l_row) {
82                      irow=j_Equa*row_block_size+l_row;
83                      for (l_col=0;l_col<col_blocksize;++l_col) {
84                         icol=j_Sol*col_blocksize+index_offset+l_col;
85                 // irow is local and icol is global
86                 Trilinos_SumIntoMyValues(in->trilinos_data, irow, icol, array[INDEX4(l_row,l_col,k_Equa,k_Sol,num_Equa,num_Sol,NN_Equa)]);
87                  }
88                    }
89                  }
90                }
91              }
92    #endif
93     } else {     } else {
94            for (k_Equa=0;k_Equa<NN_Equa;++k_Equa) {            for (k_Equa=0;k_Equa<NN_Equa;++k_Equa) { /* Down columns of array */
95              j_Equa=Nodes_Equa[k_Equa];              j_Equa=Nodes_Equa[k_Equa];
96              for (l_row=0;l_row<num_subblocks_Equa;++l_row) {              for (l_row=0;l_row<num_subblocks_Equa;++l_row) {
97                 iptr=j_Equa*num_subblocks_Equa+l_row;                 iptr=j_Equa*num_subblocks_Equa+l_row;
98                 for (k_Sol=0;k_Sol<NN_Sol;++k_Sol) {                 for (k_Sol=0;k_Sol<NN_Sol;++k_Sol) { /* Across rows of array */
99                   j_Sol=Nodes_Sol[k_Sol];                   j_Sol=Nodes_Sol[k_Sol];
100                   for (l_col=0;l_col<num_subblocks_Sol;++l_col) {                   for (l_col=0;l_col<num_subblocks_Sol;++l_col) {
101                      index=j_Sol*num_subblocks_Sol+index_offset+l_col;                      index=j_Sol*num_subblocks_Sol+index_offset+l_col;
102                  for (k=in->pattern->ptr[iptr]-index_offset;k<in->pattern->ptr[iptr+1]-index_offset;++k) {                  for (k=in->pattern->ptr[iptr]-index_offset;k<in->pattern->ptr[iptr+1]-index_offset;++k) {
103                      if (in->pattern->index[k]==index) {                      if (in->pattern->index[k]==index) {
104                            for (ic=0;ic<col_block_size;++ic) {                            for (ic=0;ic<col_block_size;++ic) { /* Entry array(k_Sol, j_Equa) is a block (row_block_size x col_block_size) */
105                                  i_Sol=ic+col_block_size*l_col;                                  i_Sol=ic+col_block_size*l_col;
106                                  for (ir=0;ir<row_block_size;++ir) {                                  for (ir=0;ir<row_block_size;++ir) {
107                                     i_Equa=ir+row_block_size*l_row;                                     i_Equa=ir+row_block_size*l_row;
108                             in->val[k*block_size+ir+row_block_size*ic]+=                             in->val[k*block_size+ir+row_block_size*ic]+=
109                                             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)];
110                       /* printf("ksteube assigning val[ %d (%d,%d) ] += %f \n", k*block_size+ir+row_block_size*ic, k_Equa, k_Sol, array[INDEX4(i_Equa,i_Sol,k_Equa,k_Sol,num_Equa,num_Sol,NN_Equa)]); */
111                                  }                                  }
112                            }                            }
113                            break;                            break;

Legend:
Removed from v.616  
changed lines
  Added in v.969

  ViewVC Help
Powered by ViewVC 1.1.26