/[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 trunk/dudley/src/IndexList.c revision 3312 by gross, Tue Oct 26 07:54: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  /* Finley: Converting an element list into a matrix shape     */  * 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    
14  /**************************************************************/  /**************************************************************/
15    
16  /*  Author: gross@access.edu.au */  /* Dudley: Converting an element list into a matrix shape     */
 /*  Version: $Id$ */  
17    
18  /**************************************************************/  /**************************************************************/
19    
# Line 25  Line 21 
21    
22  /* Translate from distributed/local array indices to global indices */  /* Translate from distributed/local array indices to global indices */
23    
 #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  
   
24  /**************************************************************/  /**************************************************************/
25  /* inserts the contributions from the element matrices of elements  /* inserts the contributions from the element matrices of elements
26     into the row index col. If symmetric is set, only the upper     into the row index col. If symmetric is set, only the upper
27     triangle of the matrix is stored. */     triangle of the matrix is stored. */
28    
29  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,
30                                         bool_t reduce_row_order, index_t* row_Label,                       bool_t reduce_row_order, index_t * row_map,
31                                         bool_t reduce_col_order, index_t* col_Label) {                       bool_t reduce_col_order, index_t * col_map)
32    /* index_list is an array of linked lists. Each entry is a row (DOF) and contains the indices to the non-zero columns */  {
33    index_t color, num_CPUs = 1, my_CPU = 0;      /* index_list is an array of linked lists. Each entry is a row (DOF) and contains the indices to the non-zero columns */
34    dim_t e,kr,kc,NN_row,NN_col,i,icol,irow;      index_t color;
35  #ifdef PASO_MPI      dim_t e, kr, kc, NN_row, NN_col, icol, irow, NN;
36      num_CPUs = mesh->MPIInfo->size;      if (elements != NULL)
37      my_CPU = mesh->MPIInfo->rank;      {
38  #endif      NN = elements->numNodes;
39    /* print_mesh_statistics( mesh, TRUE ); */      NN_col = (elements->numShapes);
40        NN_row = (elements->numShapes);
41    if (elements!=NULL) {  
42      dim_t NN=elements->ReferenceElement->Type->numNodes;      for (color = elements->minColor; color <= elements->maxColor; color++)
43      index_t id[NN],*row_node,*col_node;      {
44      for (i=0;i<NN;i++) id[i]=i;  #pragma omp for private(e,irow,kr,kc,icol) schedule(static)
45      if (reduce_col_order) {          for (e = 0; e < elements->numElements; e++)
46         col_node=elements->ReferenceElement->Type->linearNodes;          {
47         NN_col=elements->LinearReferenceElement->Type->numNodes;          if (elements->Color[e] == color)
48      } else {          {
49         col_node=id;              for (kr = 0; kr < NN_row; kr++)
50         NN_col=elements->ReferenceElement->Type->numNodes;              {
51      }              irow = row_map[elements->Nodes[INDEX2(kr, e, NN)]];
52      if (reduce_row_order) {              for (kc = 0; kc < NN_col; kc++)
53         row_node=elements->ReferenceElement->Type->linearNodes;              {
54         NN_row=elements->LinearReferenceElement->Type->numNodes;                  icol = col_map[elements->Nodes[INDEX2(kc, e, NN)]];
55      } else {                  Dudley_IndexList_insertIndex(&(index_list[irow]), icol);
56         row_node=id;              }
57         NN_row=elements->ReferenceElement->Type->numNodes;              }
58      }          }
59      if (num_CPUs == 1) {          }
60      for (color=elements->minColor;color<=elements->maxColor;color++) {      }
61          #pragma omp for private(e,irow,kr,kc,icol) schedule(static)      }
62          for (e=0;e<elements->numElements;e++) {      return;
63              if (elements->Color[e]==color) {  }
64                  for (kr=0;kr<NN_row;kr++) {  
65                    irow=row_Label[elements->Nodes[INDEX2(row_node[kr],e,NN)]];  void Dudley_IndexList_insertElementsWithRowRange(Dudley_IndexList * index_list, index_t firstRow, index_t lastRow,
66                    for (kc=0;kc<NN_col;kc++) {                           Dudley_ElementFile * elements, index_t * row_map, index_t * col_map)
67                         icol=col_Label[elements->Nodes[INDEX2(col_node[kc],e,NN)]];  {
68                         Finley_IndexList_insertIndex(&(index_list[irow]),icol);  /* this does not resolve macro elements */
69                    }      index_t color;
70                  }      dim_t e, kr, kc, icol, irow, NN;
71              }      if (elements != NULL)
72          }      {
73        }      NN = elements->numNodes;
74      }      for (color = elements->minColor; color <= elements->maxColor; color++)
75      else {  /* More than one CPU (what's below should also work for one CPU, but let's build confidence in it first) */      {
76  #ifdef PASO_MPI  #pragma omp for private(e,irow,kr,kc,icol) schedule(static)
77      Finley_NodeDistribution *row_degreeOfFreedomDistribution;          for (e = 0; e < elements->numElements; e++)
78      Finley_NodeDistribution *col_degreeOfFreedomDistribution;          {
79      if (reduce_col_order) {          if (elements->Color[e] == color)
80        col_degreeOfFreedomDistribution = mesh->Nodes->reducedDegreeOfFreedomDistribution;          {
81      }              for (kr = 0; kr < NN; kr++)
82      else {              {
83        col_degreeOfFreedomDistribution = mesh->Nodes->degreeOfFreedomDistribution;              irow = row_map[elements->Nodes[INDEX2(kr, e, NN)]];
84      }              if ((firstRow <= irow) && (irow < lastRow))
85      if (reduce_row_order) {              {
86        row_degreeOfFreedomDistribution = mesh->Nodes->reducedDegreeOfFreedomDistribution;                  irow -= firstRow;
87      }                  for (kc = 0; kc < NN; kc++)
88      else {                  {
89        row_degreeOfFreedomDistribution = mesh->Nodes->degreeOfFreedomDistribution;                  icol = col_map[elements->Nodes[INDEX2(kc, e, NN)]];
90      }                  Dudley_IndexList_insertIndex(&(index_list[irow]), icol);
91      /* Not using loop over colors as above */ {                  }
92          #pragma omp for private(e,irow,kr,kc,icol) schedule(static)              }
93      for (e=0;e<elements->numElements;e++) {              }
94                  for (kr=0;kr<NN_row;kr++) {          }
95                    irow=row_Label[elements->Nodes[INDEX2(row_node[kr],e,NN)]];          }
96            if (irow < row_degreeOfFreedomDistribution->numLocal) {      }
97                      for (kc=0;kc<NN_col;kc++) {      }
98                 /* Get the local col ID */  }
99                         icol=col_Label[elements->Nodes[INDEX2(col_node[kc],e,NN)]];  
100                 /* Convert to global col ID (row ID is saved as local value) */  void Dudley_IndexList_insertElementsWithRowRangeNoMainDiagonal(Dudley_IndexList * index_list, index_t firstRow,
101                 icol = Finley_IndexList_localToGlobal(col_degreeOfFreedomDistribution, my_CPU, icol);                                     index_t lastRow, Dudley_ElementFile * elements,
102                         Finley_IndexList_insertIndex(&(index_list[irow]),icol);                                     index_t * row_map, index_t * col_map)
103                 printf("ksteube Finley_IndexList_insertIndex cpu= %d irow= %d icol= %d\n", my_CPU, irow, icol);  {
104                      }      /* this does not resolve macro elements */
105            }      index_t color;
106                  }      dim_t e, kr, kc, icol, irow, NN, irow_loc;
107          }      if (elements != NULL)
108        }      {
109  #endif      NN = elements->numNodes;
110      }   /* More than one CPU */      for (color = elements->minColor; color <= elements->maxColor; color++)
111    }      {
112    return;  #pragma omp for private(e,irow,kr,kc,icol,irow_loc) schedule(static)
113  }          for (e = 0; e < elements->numElements; e++)
114            {
115  /* inserts row index row into the Finley_IndexList in if it does not exist */          if (elements->Color[e] == color)
116            {
117  void Finley_IndexList_insertIndex(Finley_IndexList* in, index_t index) {              for (kr = 0; kr < NN; kr++)
118    dim_t i;              {
119    /* is index in in? */              irow = row_map[elements->Nodes[INDEX2(kr, e, NN)]];
120    for (i=0;i<in->n;i++) {              if ((firstRow <= irow) && (irow < lastRow))
121      if (in->index[i]==index)  return;              {
122    }                  irow_loc = irow - firstRow;
123    /* index could not be found */                  for (kc = 0; kc < NN; kc++)
124    if (in->n==INDEXLIST_LENGTH) {                  {
125       /* if in->index is full check the extension */                  icol = col_map[elements->Nodes[INDEX2(kc, e, NN)]];
126       if (in->extension==NULL) {                  if (icol != irow)
127          in->extension=TMPMEMALLOC(1,Finley_IndexList);                      Dudley_IndexList_insertIndex(&(index_list[irow_loc]), icol);
128          if (Finley_checkPtr(in->extension)) return;                  }
129          in->extension->n=0;              }
130          in->extension->extension=NULL;              }
131       }          }
132       Finley_IndexList_insertIndex(in->extension,index);          }
133    } else {      }
134       /* insert index into in->index*/      }
135       in->index[in->n]=index;  }
136       in->n++;  
137    }  /* inserts row index row into the Dudley_IndexList in if it does not exist */
138  }  
139    void Dudley_IndexList_insertIndex(Dudley_IndexList * in, index_t index)
140  /* counts the number of row indices in the Finley_IndexList in */  {
141        dim_t i;
142  dim_t Finley_IndexList_count(Finley_IndexList* in) {      /* is index in in? */
143    if (in==NULL) {      for (i = 0; i < in->n; i++)
144       return 0;      {
145    } else {      if (in->index[i] == index)
146       return (in->n)+Finley_IndexList_count(in->extension);          return;
147    }      }
148  }      /* index could not be found */
149        if (in->n == INDEXLIST_LENGTH)
150  /* count the number of row indices in the Finley_IndexList in */      {
151        /* if in->index is full check the extension */
152  void Finley_IndexList_toArray(Finley_IndexList* in, index_t* array) {      if (in->extension == NULL)
153    dim_t i;      {
154    if (in!=NULL) {          in->extension = TMPMEMALLOC(1, Dudley_IndexList);
155      for (i=0;i<in->n;i++) array[i]=in->index[i];          if (Dudley_checkPtr(in->extension))
156      Finley_IndexList_toArray(in->extension,&(array[in->n]));          return;
157    }          in->extension->n = 0;
158  }          in->extension->extension = NULL;
159        }
160  /* deallocates the Finley_IndexList in by recursive calls */      Dudley_IndexList_insertIndex(in->extension, index);
161        }
162  void Finley_IndexList_free(Finley_IndexList* in) {      else
163    if (in!=NULL) {      {
164      Finley_IndexList_free(in->extension);      /* insert index into in->index */
165      TMPMEMFREE(in);      in->index[in->n] = index;
166    }      in->n++;
167  }      }
168    }
169  /*  
170   * $Log$  /* counts the number of row indices in the Dudley_IndexList in */
171   * Revision 1.6  2005/09/15 03:44:22  jgs  
172   * Merge of development branch dev-02 back to main trunk on 2005-09-15  dim_t Dudley_IndexList_count(Dudley_IndexList * in, index_t range_min, index_t range_max)
173   *  {
174   * Revision 1.5.2.1  2005/09/07 06:26:18  gross      dim_t i;
175   * the solver from finley are put into the standalone package paso now      dim_t out = 0;
176   *      register index_t itmp;
177   * Revision 1.5  2005/07/08 04:07:51  jgs      if (in == NULL)
178   * Merge of development branch back to main trunk on 2005-07-08      {
179   *      return 0;
180   * Revision 1.4  2004/12/15 07:08:32  jgs      }
181   * *** empty log message ***      else
182   * Revision 1.1.1.1.2.3  2005/06/29 02:34:50  gross      {
183   * some changes towards 64 integers in finley      for (i = 0; i < in->n; i++)
184   *      {
185   * Revision 1.1.1.1.2.2  2004/11/24 01:37:13  gross          itmp = in->index[i];
186   * some changes dealing with the integer overflow in memory allocation. Finley solves 4M unknowns now          if ((itmp >= range_min) && (range_max > itmp))
187   *          ++out;
188   *      }
189   *      return out + Dudley_IndexList_count(in->extension, range_min, range_max);
190   */      }
191    }
192    
193    /* count the number of row indices in the Dudley_IndexList in */
194    
195    void Dudley_IndexList_toArray(Dudley_IndexList * in, index_t * array, index_t range_min, index_t range_max,
196                      index_t index_offset)
197    {
198        dim_t i, ptr;
199        register index_t itmp;
200        if (in != NULL)
201        {
202        ptr = 0;
203        for (i = 0; i < in->n; i++)
204        {
205            itmp = in->index[i];
206            if ((itmp >= range_min) && (range_max > itmp))
207            {
208            array[ptr] = itmp + index_offset;
209            ptr++;
210            }
211    
212        }
213        Dudley_IndexList_toArray(in->extension, &(array[ptr]), range_min, range_max, index_offset);
214        }
215    }
216    
217    /* deallocates the Dudley_IndexList in by recursive calls */
218    
219    void Dudley_IndexList_free(Dudley_IndexList * in)
220    {
221        if (in != NULL)
222        {
223        Dudley_IndexList_free(in->extension);
224        TMPMEMFREE(in);
225        }
226    }
227    
228    /* creates a Paso_pattern from a range of indices */
229    Paso_Pattern *Dudley_IndexList_createPattern(dim_t n0, dim_t n, Dudley_IndexList * index_list, index_t range_min,
230                             index_t range_max, index_t index_offset)
231    {
232        dim_t *ptr = NULL;
233        register dim_t s, i, itmp;
234        index_t *index = NULL;
235        Paso_Pattern *out = NULL;
236    
237        ptr = MEMALLOC(n + 1 - n0, index_t);
238        if (!Dudley_checkPtr(ptr))
239        {
240        /* get the number of connections per row */
241    #pragma omp parallel for schedule(static) private(i)
242        for (i = n0; i < n; ++i)
243        {
244            ptr[i - n0] = Dudley_IndexList_count(&index_list[i], range_min, range_max);
245        }
246        /* accumulate ptr */
247        s = 0;
248        for (i = n0; i < n; ++i)
249        {
250            itmp = ptr[i - n0];
251            ptr[i - n0] = s;
252            s += itmp;
253        }
254        ptr[n - n0] = s;
255        /* fill index */
256        index = MEMALLOC(ptr[n - n0], index_t);
257        if (!Dudley_checkPtr(index))
258        {
259    #pragma omp parallel for schedule(static)
260            for (i = n0; i < n; ++i)
261            {
262            Dudley_IndexList_toArray(&index_list[i], &index[ptr[i - n0]], range_min, range_max, index_offset);
263            }
264            out = Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT, n - n0, range_max + index_offset, ptr, index);
265        }
266        }
267        if (!Dudley_noError())
268        {
269        MEMFREE(ptr);
270        MEMFREE(index);
271        Paso_Pattern_free(out);
272        }
273        return out;
274    }

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

  ViewVC Help
Powered by ViewVC 1.1.26