/[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

trunk/finley/src/IndexList.c revision 969 by ksteube, Tue Feb 13 23:02:23 2007 UTC branches/domexper/dudley/src/IndexList.c revision 3086 by jfenwick, Thu Aug 5 05:07:58 2010 UTC
# Line 1  Line 1 
 /*  
  ************************************************************  
  *          Copyright 2006 by ACcESS MNRF                   *  
  *                                                          *  
  *              http://www.access.edu.au                    *  
  *       Primary Business: Queensland, Australia            *  
  *  Licensed under the Open Software License version 3.0    *  
  *     http://www.opensource.org/licenses/osl-3.0.php       *  
  *                                                          *  
  ************************************************************  
 */  
1    
2  /**************************************************************/  /*******************************************************
3    *
4    * Copyright (c) 2003-2010 by University of Queensland
5    * Earth Systems Science Computational Center (ESSCC)
6    * http://www.uq.edu.au/esscc
7    *
8    * Primary Business: Queensland, Australia
9    * Licensed under the Open Software License version 3.0
10    * http://www.opensource.org/licenses/osl-3.0.php
11    *
12    *******************************************************/
13    
 /* Finley: Converting an element list into a matrix shape     */  
14    
15  /**************************************************************/  /**************************************************************/
16    
17  /*  Author: gross@access.edu.au */  /* Dudley: Converting an element list into a matrix shape     */
 /*  Version: $Id$ */  
18    
19  /**************************************************************/  /**************************************************************/
20    
# Line 25  Line 22 
22    
23  /* Translate from distributed/local array indices to global indices */  /* Translate from distributed/local array indices to global indices */
24    
 #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  
   
25  /**************************************************************/  /**************************************************************/
26  /* inserts the contributions from the element matrices of elements  /* inserts the contributions from the element matrices of elements
27     into the row index col. If symmetric is set, only the upper     into the row index col. If symmetric is set, only the upper
28     triangle of the matrix is stored. */     triangle of the matrix is stored. */
29    
30  void Finley_IndexList_insertElements(Finley_IndexList* index_list, Finley_Mesh* mesh, Finley_ElementFile* elements,  void Dudley_IndexList_insertElements(Dudley_IndexList* index_list, Dudley_ElementFile* elements,
31                                         bool_t reduce_row_order, index_t* row_Label,                                         bool_t reduce_row_order, index_t* row_map,
32                                         bool_t reduce_col_order, index_t* col_Label) {                                         bool_t reduce_col_order, index_t* col_map) {
33    /* 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 */
34    index_t color, num_CPUs = 1, my_CPU = 0;    index_t color;
35    dim_t e,kr,kc,NN_row,NN_col,i,icol,irow;    Dudley_ReferenceElement*refElement;
36  #ifdef PASO_MPI    dim_t e, kr, kc, NN_row, NN_col, icol, irow, NN, *row_node=NULL, *col_node=NULL, isub, numSub;
     num_CPUs = mesh->MPIInfo->size;  
     my_CPU = mesh->MPIInfo->rank;  
 #endif  
   /* print_mesh_statistics( mesh, TRUE ); */  
   
37    if (elements!=NULL) {    if (elements!=NULL) {
38      dim_t NN=elements->ReferenceElement->Type->numNodes;      NN=elements->numNodes;
39      index_t id[NN],*row_node,*col_node;      refElement= Dudley_ReferenceElementSet_borrowReferenceElement(elements->referenceElementSet, FALSE);
40      for (i=0;i<NN;i++) id[i]=i;        
41      if (reduce_col_order) {      if (reduce_col_order) {
42         col_node=elements->ReferenceElement->Type->linearNodes;            numSub=1;
43         NN_col=elements->LinearReferenceElement->Type->numNodes;            col_node=refElement->Type->linearNodes;
44              NN_col=(refElement->LinearBasisFunctions->Type->numShapes) * (refElement->Type->numSides) ;
45      } else {      } else {
46         col_node=id;            numSub=refElement->Type->numSubElements;
47         NN_col=elements->ReferenceElement->Type->numNodes;            col_node=refElement->Type->subElementNodes;
48              NN_col=(refElement->BasisFunctions->Type->numShapes) * (refElement->Type->numSides) ;
49      }      }
50    
51      if (reduce_row_order) {      if (reduce_row_order) {
52         row_node=elements->ReferenceElement->Type->linearNodes;            numSub=1;
53         NN_row=elements->LinearReferenceElement->Type->numNodes;            row_node=refElement->Type->linearNodes;
54              NN_row=(refElement->LinearBasisFunctions->Type->numShapes) * (refElement->Type->numSides) ;
55      } else {      } else {
56         row_node=id;            numSub=refElement->Type->numSubElements;
57         NN_row=elements->ReferenceElement->Type->numNodes;            row_node=refElement->Type->subElementNodes;
58      }            NN_row=(refElement->BasisFunctions->Type->numShapes) * (refElement->Type->numSides) ;
59      if (num_CPUs == 1) {      }
60    
61        for (color=elements->minColor;color<=elements->maxColor;color++) {
62               #pragma omp for private(e,irow,kr,kc,icol,isub) schedule(static)
63               for (e=0;e<elements->numElements;e++) {
64                   if (elements->Color[e]==color) {
65                       for (isub=0;isub<numSub; isub++) {
66                           for (kr=0;kr<NN_row;kr++) {
67                               irow=row_map[elements->Nodes[INDEX2(row_node[INDEX2(kr,isub,NN_row)],e,NN)]];
68                               for (kc=0;kc<NN_col;kc++) {
69                                   icol=col_map[elements->Nodes[INDEX2(col_node[INDEX2(kc,isub,NN_col)],e,NN)]];
70                                   Dudley_IndexList_insertIndex(&(index_list[irow]),icol);
71                               }
72                           }
73                       }
74                   }
75               }
76           }
77      }
78      return;
79    }
80    
81    
82    void Dudley_IndexList_insertElementsWithRowRange(Dudley_IndexList* index_list, index_t firstRow, index_t lastRow,
83                                                     Dudley_ElementFile* elements, index_t* row_map, index_t* col_map)
84    {
85    /* this does not resolve macro elements */
86        index_t color;
87      dim_t e,kr,kc,icol,irow, NN;
88      if (elements!=NULL) {
89        NN=elements->numNodes;
90      for (color=elements->minColor;color<=elements->maxColor;color++) {      for (color=elements->minColor;color<=elements->maxColor;color++) {
91          #pragma omp for private(e,irow,kr,kc,icol) schedule(static)             #pragma omp for private(e,irow,kr,kc,icol) schedule(static)
92          for (e=0;e<elements->numElements;e++) {             for (e=0;e<elements->numElements;e++) {
93              if (elements->Color[e]==color) {                 if (elements->Color[e]==color) {
94                  for (kr=0;kr<NN_row;kr++) {                     for (kr=0;kr<NN;kr++) {
95                    irow=row_Label[elements->Nodes[INDEX2(row_node[kr],e,NN)]];                       irow=row_map[elements->Nodes[INDEX2(kr,e,NN)]];
96                    for (kc=0;kc<NN_col;kc++) {                       if ((firstRow<=irow) && (irow < lastRow)) {
97                         icol=col_Label[elements->Nodes[INDEX2(col_node[kc],e,NN)]];                            irow-=firstRow;
98                         Finley_IndexList_insertIndex(&(index_list[irow]),icol);                            for (kc=0;kc<NN;kc++) {
99                                  icol=col_map[elements->Nodes[INDEX2(kc,e,NN)]];
100                                  Dudley_IndexList_insertIndex(&(index_list[irow]),icol);
101                              }
102                          }
103                    }                    }
104                  }                 }
105              }             }
         }  
       }  
     }  
     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;  
106      }      }
107      else {    }
108        col_degreeOfFreedomDistribution = mesh->Nodes->degreeOfFreedomDistribution;  }
109      }  void Dudley_IndexList_insertElementsWithRowRangeNoMainDiagonal(Dudley_IndexList* index_list, index_t firstRow, index_t lastRow,
110      if (reduce_row_order) {                                                                Dudley_ElementFile* elements, index_t* row_map, index_t* col_map)
111        row_degreeOfFreedomDistribution = mesh->Nodes->reducedDegreeOfFreedomDistribution;  {
112      }    /* this does not resolve macro elements */
113      else {    index_t color;
114        row_degreeOfFreedomDistribution = mesh->Nodes->degreeOfFreedomDistribution;    dim_t e,kr,kc,icol,irow, NN,irow_loc;
115      if (elements!=NULL) {
116        NN=elements->numNodes;
117        for (color=elements->minColor;color<=elements->maxColor;color++) {
118               #pragma omp for private(e,irow,kr,kc,icol,irow_loc) schedule(static)
119               for (e=0;e<elements->numElements;e++) {
120                   if (elements->Color[e]==color) {
121                       for (kr=0;kr<NN;kr++) {
122                         irow=row_map[elements->Nodes[INDEX2(kr,e,NN)]];
123                         if ((firstRow<=irow) && (irow < lastRow)) {
124                              irow_loc=irow-firstRow;
125                              for (kc=0;kc<NN;kc++) {
126                                  icol=col_map[elements->Nodes[INDEX2(kc,e,NN)]];
127                                  if (icol != irow) Dudley_IndexList_insertIndex(&(index_list[irow_loc]),icol);
128                              }
129                          }
130                      }
131                   }
132               }
133      }      }
     /* 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);  
                     }  
           }  
                 }  
         }  
       }  
 #endif  
     }   /* More than one CPU */  
134    }    }
   return;  
135  }  }
136    
137  /* inserts row index row into the Finley_IndexList in if it does not exist */  /* inserts row index row into the Dudley_IndexList in if it does not exist */
138    
139  void Finley_IndexList_insertIndex(Finley_IndexList* in, index_t index) {  void Dudley_IndexList_insertIndex(Dudley_IndexList* in, index_t index) {
140    dim_t i;    dim_t i;
141    /* is index in in? */    /* is index in in? */
142    for (i=0;i<in->n;i++) {    for (i=0;i<in->n;i++) {
# Line 146  void Finley_IndexList_insertIndex(Finley Line 146  void Finley_IndexList_insertIndex(Finley
146    if (in->n==INDEXLIST_LENGTH) {    if (in->n==INDEXLIST_LENGTH) {
147       /* if in->index is full check the extension */       /* if in->index is full check the extension */
148       if (in->extension==NULL) {       if (in->extension==NULL) {
149          in->extension=TMPMEMALLOC(1,Finley_IndexList);          in->extension=TMPMEMALLOC(1,Dudley_IndexList);
150          if (Finley_checkPtr(in->extension)) return;          if (Dudley_checkPtr(in->extension)) return;
151          in->extension->n=0;          in->extension->n=0;
152          in->extension->extension=NULL;          in->extension->extension=NULL;
153       }       }
154       Finley_IndexList_insertIndex(in->extension,index);       Dudley_IndexList_insertIndex(in->extension,index);
155    } else {    } else {
156       /* insert index into in->index*/       /* insert index into in->index*/
157       in->index[in->n]=index;       in->index[in->n]=index;
# Line 159  void Finley_IndexList_insertIndex(Finley Line 159  void Finley_IndexList_insertIndex(Finley
159    }    }
160  }  }
161    
162  /* counts the number of row indices in the Finley_IndexList in */  /* counts the number of row indices in the Dudley_IndexList in */
163    
164  dim_t Finley_IndexList_count(Finley_IndexList* in) {  dim_t Dudley_IndexList_count(Dudley_IndexList* in, index_t range_min,index_t range_max) {
165      dim_t i;
166      dim_t out=0;
167      register index_t itmp;
168    if (in==NULL) {    if (in==NULL) {
169       return 0;       return 0;
170    } else {    } else {
171       return (in->n)+Finley_IndexList_count(in->extension);      for (i=0;i<in->n;i++) {
172              itmp=in->index[i];
173              if ((itmp>=range_min) && (range_max>itmp)) ++out;
174        }
175         return out+Dudley_IndexList_count(in->extension, range_min,range_max);
176    }    }
177  }  }
178    
179  /* count the number of row indices in the Finley_IndexList in */  /* count the number of row indices in the Dudley_IndexList in */
180    
181  void Finley_IndexList_toArray(Finley_IndexList* in, index_t* array) {  void Dudley_IndexList_toArray(Dudley_IndexList* in, index_t* array, index_t range_min,index_t range_max, index_t index_offset) {
182    dim_t i;    dim_t i, ptr;
183      register index_t itmp;
184    if (in!=NULL) {    if (in!=NULL) {
185      for (i=0;i<in->n;i++) array[i]=in->index[i];      ptr=0;
186      Finley_IndexList_toArray(in->extension,&(array[in->n]));      for (i=0;i<in->n;i++) {
187              itmp=in->index[i];
188              if ((itmp>=range_min) && (range_max>itmp)) {
189                 array[ptr]=itmp+index_offset;
190                 ptr++;
191              }
192    
193        }
194        Dudley_IndexList_toArray(in->extension,&(array[ptr]), range_min, range_max, index_offset);
195    }    }
196  }  }
197    
198  /* deallocates the Finley_IndexList in by recursive calls */  /* deallocates the Dudley_IndexList in by recursive calls */
199    
200  void Finley_IndexList_free(Finley_IndexList* in) {  void Dudley_IndexList_free(Dudley_IndexList* in) {
201    if (in!=NULL) {    if (in!=NULL) {
202      Finley_IndexList_free(in->extension);      Dudley_IndexList_free(in->extension);
203      TMPMEMFREE(in);      TMPMEMFREE(in);
204    }    }
205  }  }
206    
207  /*  /* creates a Paso_pattern from a range of indices */
208   * $Log$  Paso_Pattern* Dudley_IndexList_createPattern(dim_t n0, dim_t n,Dudley_IndexList* index_list,index_t range_min,index_t range_max,index_t index_offset)
209   * Revision 1.6  2005/09/15 03:44:22  jgs  {
210   * Merge of development branch dev-02 back to main trunk on 2005-09-15     dim_t *ptr=NULL;
211   *     register dim_t s,i,itmp;
212   * Revision 1.5.2.1  2005/09/07 06:26:18  gross     index_t *index=NULL;
213   * the solver from finley are put into the standalone package paso now     Paso_Pattern* out=NULL;
214   *  
215   * Revision 1.5  2005/07/08 04:07:51  jgs     ptr=MEMALLOC(n+1-n0,index_t);
216   * Merge of development branch back to main trunk on 2005-07-08     if (! Dudley_checkPtr(ptr) ) {
217   *         /* get the number of connections per row */
218   * Revision 1.4  2004/12/15 07:08:32  jgs         #pragma omp parallel for schedule(static) private(i)
219   * *** empty log message ***         for(i=n0;i<n;++i) {
220   * Revision 1.1.1.1.2.3  2005/06/29 02:34:50  gross                ptr[i-n0]=Dudley_IndexList_count(&index_list[i],range_min,range_max);
221   * some changes towards 64 integers in finley         }
222   *         /* accumulate ptr */
223   * Revision 1.1.1.1.2.2  2004/11/24 01:37:13  gross         s=0;
224   * some changes dealing with the integer overflow in memory allocation. Finley solves 4M unknowns now         for(i=n0;i<n;++i) {
225   *                 itmp=ptr[i-n0];
226   *                 ptr[i-n0]=s;
227   *                 s+=itmp;
228   */         }
229           ptr[n-n0]=s;
230           /* fill index */
231           index=MEMALLOC(ptr[n-n0],index_t);
232           if (! Dudley_checkPtr(index)) {
233                  #pragma omp parallel for schedule(static)
234                  for(i=n0;i<n;++i) {
235                      Dudley_IndexList_toArray(&index_list[i],&index[ptr[i-n0]],range_min,range_max,index_offset);
236                  }
237                  out=Paso_Pattern_alloc(PATTERN_FORMAT_DEFAULT,n-n0,range_max+index_offset,ptr,index);
238           }
239      }
240      if (! Dudley_noError()) {
241            MEMFREE(ptr);
242            MEMFREE(index);
243            Paso_Pattern_free(out);
244      }
245      return out;
246    }

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

  ViewVC Help
Powered by ViewVC 1.1.26