/[escript]/trunk-mpi-branch/finley/src/IndexList.c
ViewVC logotype

Diff of /trunk-mpi-branch/finley/src/IndexList.c

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

revision 1222 by ksteube, Tue May 15 03:23:17 2007 UTC revision 1223 by gross, Fri Aug 3 02:40:39 2007 UTC
# Line 25  Line 25 
25    
26  /* Translate from distributed/local array indices to global indices */  /* Translate from distributed/local array indices to global indices */
27    
 int Finley_IndexList_localToGlobal(Finley_NodeDistribution *dofDistribution, 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)  
   */  
   index_t my_CPU=dofDistribution->MPIInfo->rank;  
   if (localIndex < dofDistribution->numLocal) {  
     localIndex = localIndex + dofDistribution->vtxdist[my_CPU];  
   }  
   else {  
     localIndex = dofDistribution->indexExternal[localIndex-dofDistribution->numLocal];  
   }  
   return(localIndex);  
 }  
   
28  /**************************************************************/  /**************************************************************/
29  /* inserts the contributions from the element matrices of elements  /* inserts the contributions from the element matrices of elements
30     into the row index col. If symmetric is set, only the upper     into the row index col. If symmetric is set, only the upper
31     triangle of the matrix is stored. */     triangle of the matrix is stored. */
32    
33  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,
34                                         bool_t reduce_row_order, index_t* row_Label,                                         bool_t reduce_row_order, index_t* row_map,
35                                         bool_t reduce_col_order, index_t* col_Label) {                                         bool_t reduce_col_order, index_t* col_map) {
36    /* index_list is an array of linked lists. Each entry is a row (DOF) and contains the indices to the non-zero columns */    /* index_list is an array of linked lists. Each entry is a row (DOF) and contains the indices to the non-zero columns */
37    index_t color, *id=NULL, num_CPUs = 1, my_CPU = 0;    index_t color, *id=NULL;
38    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, NN, *row_node=NULL,*col_node=NULL;
 #ifdef PASO_MPI  
     num_CPUs = mesh->MPIInfo->size;  
     my_CPU = mesh->MPIInfo->rank;  
 #endif  
   /* print_mesh_statistics( mesh, TRUE ); */  
   
39    if (elements!=NULL) {    if (elements!=NULL) {
40      dim_t NN=elements->ReferenceElement->Type->numNodes;      NN=elements->numNodes;
     index_t *row_node,*col_node;  
41      id=TMPMEMALLOC(NN, index_t);      id=TMPMEMALLOC(NN, index_t);
42      if (! Finley_checkPtr(id) ) {      if (! Finley_checkPtr(id) ) {
43      for (i=0;i<NN;i++) id[i]=i;         for (i=0;i<NN;i++) id[i]=i;
44      if (reduce_col_order) {         if (reduce_col_order) {
45         col_node=elements->ReferenceElement->Type->linearNodes;            col_node=elements->ReferenceElement->Type->linearNodes;
46         NN_col=elements->LinearReferenceElement->Type->numNodes;            NN_col=elements->LinearReferenceElement->Type->numNodes;
47      } else {         } else {
48         col_node=id;            col_node=id;
49         NN_col=elements->ReferenceElement->Type->numNodes;            NN_col=elements->ReferenceElement->Type->numNodes;
50      }         }
51      if (reduce_row_order) {         if (reduce_row_order) {
52         row_node=elements->ReferenceElement->Type->linearNodes;            row_node=elements->ReferenceElement->Type->linearNodes;
53         NN_row=elements->LinearReferenceElement->Type->numNodes;            NN_row=elements->LinearReferenceElement->Type->numNodes;
54      } else {         } else {
55         row_node=id;            row_node=id;
56         NN_row=elements->ReferenceElement->Type->numNodes;            NN_row=elements->ReferenceElement->Type->numNodes;
57      }         }
58      if (num_CPUs == 1) {         for (color=elements->minColor;color<=elements->maxColor;color++) {
59      for (color=elements->minColor;color<=elements->maxColor;color++) {             #pragma omp for private(e,irow,kr,kc,icol) schedule(static)
60          #pragma omp for private(e,irow,kr,kc,icol) schedule(static)             for (e=0;e<elements->numElements;e++) {
61          for (e=0;e<elements->numElements;e++) {                 if (elements->Color[e]==color) {
62              if (elements->Color[e]==color) {                     for (kr=0;kr<NN_row;kr++) {
63                  for (kr=0;kr<NN_row;kr++) {                       irow=row_map[elements->Nodes[INDEX2(row_node[kr],e,NN)]];
64                    irow=row_Label[elements->Nodes[INDEX2(row_node[kr],e,NN)]];                       for (kc=0;kc<NN_col;kc++) {
65                    for (kc=0;kc<NN_col;kc++) {                            icol=col_map[elements->Nodes[INDEX2(col_node[kc],e,NN)]];
66                         icol=col_Label[elements->Nodes[INDEX2(col_node[kc],e,NN)]];                            Finley_IndexList_insertIndex(&(index_list[irow]),icol);
67                         Finley_IndexList_insertIndex(&(index_list[irow]),icol);                       }
68                    }                     }
69                  }                 }
70              }             }
71          }         }
72        }         TMPMEMFREE(id);
     } else {    /* More than one CPU (what's below should also work for one CPU, but let's build confidence in it first) */  
 #ifdef PASO_MPI  
     Finley_NodeDistribution *row_degreeOfFreedomDistribution;  
     Finley_NodeDistribution *col_degreeOfFreedomDistribution;  
     if (reduce_col_order) {  
       col_degreeOfFreedomDistribution = mesh->Nodes->reducedDegreeOfFreedomDistribution;  
     }  
     else {  
       col_degreeOfFreedomDistribution = mesh->Nodes->degreeOfFreedomDistribution;  
73      }      }
     if (reduce_row_order) {  
       row_degreeOfFreedomDistribution = mesh->Nodes->reducedDegreeOfFreedomDistribution;  
     }  
     else {  
       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, icol);  
                        Finley_IndexList_insertIndex(&(index_list[irow]),icol);  
                     }  
           }  
                 }  
         }  
       }  
 #endif  
     }   /* More than one CPU */  
     TMPMEMFREE(id);  
   }  
74    }    }
75    return;    return;
76  }  }
# Line 162  void Finley_IndexList_insertIndex(Finley Line 102  void Finley_IndexList_insertIndex(Finley
102    
103  /* counts the number of row indices in the Finley_IndexList in */  /* counts the number of row indices in the Finley_IndexList in */
104    
105  dim_t Finley_IndexList_count(Finley_IndexList* in) {  dim_t Finley_IndexList_count(Finley_IndexList* in, index_t range_min,index_t range_max) {
106      dim_t i;
107      dim_t out=0;
108      register index_t itmp;
109    if (in==NULL) {    if (in==NULL) {
110       return 0;       return 0;
111    } else {    } else {
112       return (in->n)+Finley_IndexList_count(in->extension);      for (i=0;i<in->n;i++) {
113              itmp=in->index[i];
114              if ((itmp>=range_min) && (range_max>itmp)) ++out;
115        }
116         return out+Finley_IndexList_count(in->extension, range_min,range_max);
117    }    }
118  }  }
119    
120  /* count the number of row indices in the Finley_IndexList in */  /* count the number of row indices in the Finley_IndexList in */
121    
122  void Finley_IndexList_toArray(Finley_IndexList* in, index_t* array) {  void Finley_IndexList_toArray(Finley_IndexList* in, index_t* array, index_t range_min,index_t range_max) {
123    dim_t i;    dim_t i;
124      register index_t itmp;
125    if (in!=NULL) {    if (in!=NULL) {
126      for (i=0;i<in->n;i++) array[i]=in->index[i];      for (i=0;i<in->n;i++) {
127      Finley_IndexList_toArray(in->extension,&(array[in->n]));            itmp=in->index[i];
128              if ((itmp>=range_min) && (range_max>itmp)) array[i]=itmp;
129        }
130        Finley_IndexList_toArray(in->extension,&(array[in->n]), range_min, range_max);
131    }    }
132  }  }
133    
# Line 189  void Finley_IndexList_free(Finley_IndexL Line 140  void Finley_IndexList_free(Finley_IndexL
140    }    }
141  }  }
142    
143    /* creates a Paso_pattern from a range of indices */
144    Paso_Pattern* Finley_IndexList_createPattern(dim_t n,Finley_IndexList* index_list,index_t range_min,index_t range_max)
145    {
146       dim_t *ptr=NULL;
147       register dim_t s,i,itmp;
148       index_t *index=NULL;
149       Paso_Pattern* out=NULL;
150    
151       ptr=MEMALLOC(n+1,index_t);
152       if (! Finley_checkPtr(ptr) ) {
153           /* get the number of connections per row */
154           #pragma omp parallel for schedule(static) private(i)
155           for(i=0;i<n;++i) {
156                  ptr[i]=Finley_IndexList_count(&index_list[i],range_min,range_max);
157           }
158           /* accumulate ptr */
159           s=0;
160           for(i=0;i<n;++i) {
161                   itmp=ptr[i];
162                   ptr[i]=s;
163                   s+=itmp;
164           }
165           ptr[n]=s;
166           /* fill index */
167           index=MEMALLOC(ptr[n],index_t);
168           if (! Finley_checkPtr(index)) {
169                  #pragma omp parallel for schedule(static)
170                  for(i=0;i<n;++i) Finley_IndexList_toArray(&index_list[i],&index[ptr[i]],range_min,range_max);
171    
172                  out=Paso_Pattern_alloc(PATTERN_FORMAT_DEFAULT,n,ptr,index);
173           }
174      }
175      if (! Finley_noError()) {
176            MEMFREE(ptr);
177            MEMFREE(index);
178            Paso_Pattern_free(out);
179      }
180      return out;
181    }
182    
183    Paso_Coupler* Finley_IndexList_createCoupler(Finley_IndexList *index_list)
184    {
185        Paso_Coupler* out=NULL;
186        index_t *sharedInput=NULL, *offsetInSharedInput=NULL, *offsetInRemoteInput=NULL;
187  /*  /*
188   * $Log$  
189   * Revision 1.6  2005/09/15 03:44:22  jgs      if (! (rowMap->MPIInfo->comm == colMap->MPIInfo->comm ) ) {
190   * Merge of development branch dev-02 back to main trunk on 2005-09-15          Finley_setError(SYSTEM_ERROR,"Finley_IndexList_createCoupler: communicator for row and column DOFMap must be identical.");
191   *          return NULL;
192   * Revision 1.5.2.1  2005/09/07 06:26:18  gross      }
193   * the solver from finley are put into the standalone package paso now      if (! (rowMap->numNeighbours == colMap->numNeighbours ) ) {
194   *          Finley_setError(SYSTEM_ERROR,"Finley_IndexList_createCoupler: number of neighbours for row and column DOFMap must be identical.");
195   * Revision 1.5  2005/07/08 04:07:51  jgs          return NULL;
196   * Merge of development branch back to main trunk on 2005-07-08      }
197   *      
198   * Revision 1.4  2004/12/15 07:08:32  jgs  
199   * *** empty log message ***      out=Paso_Coupler_alloc(rowMap->numNeighbours,
200   * Revision 1.1.1.1.2.3  2005/06/29 02:34:50  gross                             rowMap->neighbours,
201   * some changes towards 64 integers in finley                             sharedInput,
202   *                             offsetInSharedInput,
203   * Revision 1.1.1.1.2.2  2004/11/24 01:37:13  gross                             offsetInRemoteInput,
204   * some changes dealing with the integer overflow in memory allocation. Finley solves 4M unknowns now                             rowMap->MPIInfo);
205   *  */
206   *     return out;
207   *  }
  */  

Legend:
Removed from v.1222  
changed lines
  Added in v.1223

  ViewVC Help
Powered by ViewVC 1.1.26