/[escript]/trunk/dudley/src/IndexList.cpp
ViewVC logotype

Diff of /trunk/dudley/src/IndexList.cpp

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

revision 969 by ksteube, Tue Feb 13 23:02:23 2007 UTC revision 1028 by gross, Wed Mar 14 00:15:24 2007 UTC
# Line 23  Line 23 
23    
24  #include "IndexList.h"  #include "IndexList.h"
25    
 /* Translate from distributed/local array indices to global indices */  
   
 #ifdef PASO_MPI  
 int Finley_IndexList_localToGlobal(Finley_NodeDistribution *dofDistribution, int my_CPU, int localIndex) {  
   /*  
     get global id of icol  
     if icol is internal node (on this CPU): use icol+vtxdist[my_CPU]  
     else use indexExternal[icol-numLocal] to get global index of node  
     (actually DOF...the NodeDistribution structure should have been called DofDistribution)  
   */  
   if (localIndex < dofDistribution->numLocal) {  
     localIndex = localIndex + dofDistribution->vtxdist[my_CPU];  
   }  
   else {  
     localIndex = dofDistribution->indexExternal[localIndex-dofDistribution->numLocal];  
   }  
   return(localIndex);  
 }  
 #endif  
   
26  /**************************************************************/  /**************************************************************/
27  /* inserts the contributions from the element matrices of elements  /* inserts the contributions from the element matrices of elements
28     into the row index col. If symmetric is set, only the upper     into the row index col. If symmetric is set, only the upper
29     triangle of the matrix is stored. */     triangle of the matrix is stored. */
30    
31  void Finley_IndexList_insertElements(Finley_IndexList* index_list, Finley_Mesh* mesh, Finley_ElementFile* elements,  void Finley_IndexList_insertElements(Finley_IndexList* index_list, Finley_ElementFile* elements,
32                                         bool_t reduce_row_order, index_t* row_Label,                                         bool_t reduce_row_order, index_t* row_Label,
33                                         bool_t reduce_col_order, index_t* col_Label) {                                         bool_t reduce_col_order, index_t* col_Label) {
34    /* index_list is an array of linked lists. Each entry is a row (DOF) and contains the indices to the non-zero columns */    index_t color, *id=NULL;
35    index_t color, num_CPUs = 1, my_CPU = 0;    dim_t e,kr,kc,NN_row,NN_col,i,icol,irow, NN, *row_node=NULL,*col_node=NULL;
   dim_t e,kr,kc,NN_row,NN_col,i,icol,irow;  
 #ifdef PASO_MPI  
     num_CPUs = mesh->MPIInfo->size;  
     my_CPU = mesh->MPIInfo->rank;  
 #endif  
   /* print_mesh_statistics( mesh, TRUE ); */  
36    
37    if (elements!=NULL) {    if (elements!=NULL) {
38      dim_t NN=elements->ReferenceElement->Type->numNodes;      NN=elements->ReferenceElement->Type->numNodes;
39      index_t id[NN],*row_node,*col_node;      id=TMPMEMALLOC(NN, index_t);
40      for (i=0;i<NN;i++) id[i]=i;      if (! Finley_checkPtr(id) ) {
41      if (reduce_col_order) {          for (i=0;i<NN;i++) id[i]=i;
42         col_node=elements->ReferenceElement->Type->linearNodes;          if (reduce_col_order) {
43         NN_col=elements->LinearReferenceElement->Type->numNodes;             col_node=elements->ReferenceElement->Type->linearNodes;
44      } else {             NN_col=elements->LinearReferenceElement->Type->numNodes;
45         col_node=id;          } else {
46         NN_col=elements->ReferenceElement->Type->numNodes;             col_node=id;
47      }             NN_col=elements->ReferenceElement->Type->numNodes;
     if (reduce_row_order) {  
        row_node=elements->ReferenceElement->Type->linearNodes;  
        NN_row=elements->LinearReferenceElement->Type->numNodes;  
     } else {  
        row_node=id;  
        NN_row=elements->ReferenceElement->Type->numNodes;  
     }  
     if (num_CPUs == 1) {  
     for (color=elements->minColor;color<=elements->maxColor;color++) {  
         #pragma omp for private(e,irow,kr,kc,icol) schedule(static)  
         for (e=0;e<elements->numElements;e++) {  
             if (elements->Color[e]==color) {  
                 for (kr=0;kr<NN_row;kr++) {  
                   irow=row_Label[elements->Nodes[INDEX2(row_node[kr],e,NN)]];  
                   for (kc=0;kc<NN_col;kc++) {  
                        icol=col_Label[elements->Nodes[INDEX2(col_node[kc],e,NN)]];  
                        Finley_IndexList_insertIndex(&(index_list[irow]),icol);  
                   }  
                 }  
             }  
48          }          }
49        }          if (reduce_row_order) {
50      }             row_node=elements->ReferenceElement->Type->linearNodes;
51      else {  /* More than one CPU (what's below should also work for one CPU, but let's build confidence in it first) */             NN_row=elements->LinearReferenceElement->Type->numNodes;
52  #ifdef PASO_MPI          } else {
53      Finley_NodeDistribution *row_degreeOfFreedomDistribution;             row_node=id;
54      Finley_NodeDistribution *col_degreeOfFreedomDistribution;             NN_row=elements->ReferenceElement->Type->numNodes;
55      if (reduce_col_order) {          }
56        col_degreeOfFreedomDistribution = mesh->Nodes->reducedDegreeOfFreedomDistribution;          for (color=elements->minColor;color<=elements->maxColor;color++) {
57      }              #pragma omp for private(e,irow,kr,kc,icol) schedule(static)
58      else {              for (e=0;e<elements->numElements;e++) {
59        col_degreeOfFreedomDistribution = mesh->Nodes->degreeOfFreedomDistribution;                  if (elements->Color[e]==color) {
60      }                      for (kr=0;kr<NN_row;kr++) {
61      if (reduce_row_order) {                        irow=row_Label[elements->Nodes[INDEX2(row_node[kr],e,NN)]];
62        row_degreeOfFreedomDistribution = mesh->Nodes->reducedDegreeOfFreedomDistribution;                        for (kc=0;kc<NN_col;kc++) {
63      }                             icol=col_Label[elements->Nodes[INDEX2(col_node[kc],e,NN)]];
64      else {                             Finley_IndexList_insertIndex(&(index_list[irow]),icol);
65        row_degreeOfFreedomDistribution = mesh->Nodes->degreeOfFreedomDistribution;                        }
     }  
     /* Not using loop over colors as above */ {  
         #pragma omp for private(e,irow,kr,kc,icol) schedule(static)  
     for (e=0;e<elements->numElements;e++) {  
                 for (kr=0;kr<NN_row;kr++) {  
                   irow=row_Label[elements->Nodes[INDEX2(row_node[kr],e,NN)]];  
           if (irow < row_degreeOfFreedomDistribution->numLocal) {  
                     for (kc=0;kc<NN_col;kc++) {  
                /* Get the local col ID */  
                        icol=col_Label[elements->Nodes[INDEX2(col_node[kc],e,NN)]];  
                /* Convert to global col ID (row ID is saved as local value) */  
                icol = Finley_IndexList_localToGlobal(col_degreeOfFreedomDistribution, my_CPU, icol);  
                        Finley_IndexList_insertIndex(&(index_list[irow]),icol);  
                printf("ksteube Finley_IndexList_insertIndex cpu= %d irow= %d icol= %d\n", my_CPU, irow, icol);  
66                      }                      }
           }  
67                  }                  }
68                }
69          }          }
70            TMPMEMFREE(id);
71        }        }
 #endif  
     }   /* More than one CPU */  
72    }    }
73    return;    return;
74  }  }

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

  ViewVC Help
Powered by ViewVC 1.1.26