/[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 1552 by gross, Thu May 8 08:52:41 2008 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  /**************************************************************/  /* $Id$ */
3    
4  /* Finley: Converting an element list into a matrix shape     */  /*******************************************************
5     *
6     *           Copyright 2003-2007 by ACceSS MNRF
7     *       Copyright 2007 by University of Queensland
8     *
9     *                http://esscc.uq.edu.au
10     *        Primary Business: Queensland, Australia
11     *  Licensed under the Open Software License version 3.0
12     *     http://www.opensource.org/licenses/osl-3.0.php
13     *
14     *******************************************************/
15    
16  /**************************************************************/  /**************************************************************/
17    
18  /*  Author: gross@access.edu.au */  /* Finley: Converting an element list into a matrix shape     */
 /*  Version: $Id$ */  
19    
20  /**************************************************************/  /**************************************************************/
21    
# Line 25  Line 23 
23    
24  /* Translate from distributed/local array indices to global indices */  /* Translate from distributed/local array indices to global indices */
25    
 #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_map,
33                                         bool_t reduce_col_order, index_t* col_Label) {                                         bool_t reduce_col_order, index_t* col_map) {
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_list is an array of linked lists. Each entry is a row (DOF) and contains the indices to the non-zero columns */
35    index_t color, num_CPUs = 1, my_CPU = 0;    index_t color, *id=NULL;
36    dim_t e,kr,kc,NN_row,NN_col,i,icol,irow;    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 ); */  
   
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;      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;
48      if (reduce_row_order) {         }
49         row_node=elements->ReferenceElement->Type->linearNodes;         if (reduce_row_order) {
50         NN_row=elements->LinearReferenceElement->Type->numNodes;            row_node=elements->ReferenceElement->Type->linearNodes;
51      } else {            NN_row=elements->LinearReferenceElement->Type->numNodes;
52         row_node=id;         } else {
53         NN_row=elements->ReferenceElement->Type->numNodes;            row_node=id;
54              NN_row=elements->ReferenceElement->Type->numNodes;
55           }
56           for (color=elements->minColor;color<=elements->maxColor;color++) {
57               #pragma omp for private(e,irow,kr,kc,icol) schedule(static)
58               for (e=0;e<elements->numElements;e++) {
59                   if (elements->Color[e]==color) {
60                       for (kr=0;kr<NN_row;kr++) {
61                         irow=row_map[elements->Nodes[INDEX2(row_node[kr],e,NN)]];
62                         for (kc=0;kc<NN_col;kc++) {
63                              icol=col_map[elements->Nodes[INDEX2(col_node[kc],e,NN)]];
64                              Finley_IndexList_insertIndex(&(index_list[irow]),icol);
65                         }
66                       }
67                   }
68               }
69           }
70           TMPMEMFREE(id);
71      }      }
72      if (num_CPUs == 1) {    }
73      return;
74    }
75    
76    void Finley_IndexList_insertElementsWithRowRange(Finley_IndexList* index_list, index_t firstRow, index_t lastRow,
77                                                     Finley_ElementFile* elements, index_t* row_map, index_t* col_map)
78    {
79      index_t color;
80      dim_t e,kr,kc,i,icol,irow, NN;
81      if (elements!=NULL) {
82        NN=elements->numNodes;
83      for (color=elements->minColor;color<=elements->maxColor;color++) {      for (color=elements->minColor;color<=elements->maxColor;color++) {
84          #pragma omp for private(e,irow,kr,kc,icol) schedule(static)             #pragma omp for private(e,irow,kr,kc,icol) schedule(static)
85          for (e=0;e<elements->numElements;e++) {             for (e=0;e<elements->numElements;e++) {
86              if (elements->Color[e]==color) {                 if (elements->Color[e]==color) {
87                  for (kr=0;kr<NN_row;kr++) {                     for (kr=0;kr<NN;kr++) {
88                    irow=row_Label[elements->Nodes[INDEX2(row_node[kr],e,NN)]];                       irow=row_map[elements->Nodes[INDEX2(kr,e,NN)]];
89                    for (kc=0;kc<NN_col;kc++) {                       if ((firstRow<=irow) && (irow < lastRow)) {
90                         icol=col_Label[elements->Nodes[INDEX2(col_node[kc],e,NN)]];                            irow-=firstRow;
91                         Finley_IndexList_insertIndex(&(index_list[irow]),icol);                            for (kc=0;kc<NN;kc++) {
92                                  icol=col_map[elements->Nodes[INDEX2(kc,e,NN)]];
93                                  Finley_IndexList_insertIndex(&(index_list[irow]),icol);
94                              }
95                          }
96                    }                    }
97                  }                 }
98              }             }
         }  
       }  
     }  
     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;  
99      }      }
     else {  
       col_degreeOfFreedomDistribution = mesh->Nodes->degreeOfFreedomDistribution;  
     }  
     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, 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 */  
100    }    }
   return;  
101  }  }
102    
103  /* inserts row index row into the Finley_IndexList in if it does not exist */  /* inserts row index row into the Finley_IndexList in if it does not exist */
# Line 161  void Finley_IndexList_insertIndex(Finley Line 127  void Finley_IndexList_insertIndex(Finley
127    
128  /* counts the number of row indices in the Finley_IndexList in */  /* counts the number of row indices in the Finley_IndexList in */
129    
130  dim_t Finley_IndexList_count(Finley_IndexList* in) {  dim_t Finley_IndexList_count(Finley_IndexList* in, index_t range_min,index_t range_max) {
131      dim_t i;
132      dim_t out=0;
133      register index_t itmp;
134    if (in==NULL) {    if (in==NULL) {
135       return 0;       return 0;
136    } else {    } else {
137       return (in->n)+Finley_IndexList_count(in->extension);      for (i=0;i<in->n;i++) {
138              itmp=in->index[i];
139              if ((itmp>=range_min) && (range_max>itmp)) ++out;
140        }
141         return out+Finley_IndexList_count(in->extension, range_min,range_max);
142    }    }
143  }  }
144    
145  /* count the number of row indices in the Finley_IndexList in */  /* count the number of row indices in the Finley_IndexList in */
146    
147  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, index_t index_offset) {
148    dim_t i;    dim_t i, ptr;
149      register index_t itmp;
150    if (in!=NULL) {    if (in!=NULL) {
151      for (i=0;i<in->n;i++) array[i]=in->index[i];      ptr=0;
152      Finley_IndexList_toArray(in->extension,&(array[in->n]));      for (i=0;i<in->n;i++) {
153              itmp=in->index[i];
154              if ((itmp>=range_min) && (range_max>itmp)) {
155                 array[ptr]=itmp+index_offset;
156                 ptr++;
157              }
158    
159        }
160        Finley_IndexList_toArray(in->extension,&(array[ptr]), range_min, range_max, index_offset);
161    }    }
162  }  }
163    
# Line 188  void Finley_IndexList_free(Finley_IndexL Line 170  void Finley_IndexList_free(Finley_IndexL
170    }    }
171  }  }
172    
173  /*  /* creates a Paso_pattern from a range of indices */
174   * $Log$  Paso_Pattern* Finley_IndexList_createPattern(dim_t n0, dim_t n,Finley_IndexList* index_list,index_t range_min,index_t range_max,index_t index_offset)
175   * Revision 1.6  2005/09/15 03:44:22  jgs  {
176   * Merge of development branch dev-02 back to main trunk on 2005-09-15     dim_t *ptr=NULL;
177   *     register dim_t s,i,itmp;
178   * Revision 1.5.2.1  2005/09/07 06:26:18  gross     index_t *index=NULL;
179   * the solver from finley are put into the standalone package paso now     Paso_Pattern* out=NULL;
180   *  
181   * Revision 1.5  2005/07/08 04:07:51  jgs     ptr=MEMALLOC(n+1-n0,index_t);
182   * Merge of development branch back to main trunk on 2005-07-08     if (! Finley_checkPtr(ptr) ) {
183   *         /* get the number of connections per row */
184   * Revision 1.4  2004/12/15 07:08:32  jgs         #pragma omp parallel for schedule(static) private(i)
185   * *** empty log message ***         for(i=n0;i<n;++i) {
186   * Revision 1.1.1.1.2.3  2005/06/29 02:34:50  gross                ptr[i-n0]=Finley_IndexList_count(&index_list[i],range_min,range_max);
187   * some changes towards 64 integers in finley         }
188   *         /* accumulate ptr */
189   * Revision 1.1.1.1.2.2  2004/11/24 01:37:13  gross         s=0;
190   * some changes dealing with the integer overflow in memory allocation. Finley solves 4M unknowns now         for(i=n0;i<n;++i) {
191   *                 itmp=ptr[i-n0];
192   *                 ptr[i-n0]=s;
193   *                 s+=itmp;
194   */         }
195           ptr[n-n0]=s;
196           /* fill index */
197           index=MEMALLOC(ptr[n-n0],index_t);
198           if (! Finley_checkPtr(index)) {
199                  #pragma omp parallel for schedule(static)
200                  for(i=n0;i<n;++i) {
201                      Finley_IndexList_toArray(&index_list[i],&index[ptr[i-n0]],range_min,range_max,index_offset);
202                  }
203                  out=Paso_Pattern_alloc(PATTERN_FORMAT_DEFAULT,1,1,n-n0,ptr,index);
204           }
205      }
206      if (! Finley_noError()) {
207            MEMFREE(ptr);
208            MEMFREE(index);
209            Paso_Pattern_free(out);
210      }
211      return out;
212    }

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

  ViewVC Help
Powered by ViewVC 1.1.26