/[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/finley/IndexList.c revision 201 by jgs, Wed Nov 23 04:10:21 2005 UTC trunk/dudley/src/IndexList.c revision 3981 by jfenwick, Fri Sep 21 02:47:54 2012 UTC
# Line 1  Line 1 
 /*  
  ******************************************************************************  
  *                                                                            *  
  *       COPYRIGHT  ACcESS 2003,2004,2005 -  All Rights Reserved              *  
  *                                                                            *  
  * This software is the property of ACcESS. No part of this code              *  
  * may be copied in any form or by any means without the expressed written    *  
  * consent of ACcESS.  Copying, use or modification of this software          *  
  * by any unauthorised person is illegal unless that person has a software    *  
  * license agreement with ACcESS.                                             *  
  *                                                                            *  
  ******************************************************************************  
 */  
   
 /**************************************************************/  
1    
2  /* Finley: Converting an element list into a matrix shape     */  /*****************************************************************************
3    *
4    * Copyright (c) 2003-2012 by University of Queensland
5    * http://www.uq.edu.au
6    *
7    * Primary Business: Queensland, Australia
8    * Licensed under the Open Software License version 3.0
9    * http://www.opensource.org/licenses/osl-3.0.php
10    *
11    * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12    * Development since 2012 by School of Earth Sciences
13    *
14    *****************************************************************************/
15    
16  /**************************************************************/  /************************************************************************************/
17    
18  /*  Author: gross@access.edu.au */  /* Dudley: Converting an element list into a matrix shape     */
 /*  Version: $Id$ */  
19    
20  /**************************************************************/  /************************************************************************************/
21    
22  #include "IndexList.h"  #include "IndexList.h"
23    
24  /**************************************************************/  /* Translate from distributed/local array indices to global indices */
25    
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_ElementFile* elements,  void Dudley_IndexList_insertElements(Dudley_IndexList * index_list, Dudley_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_t color;  {
35    dim_t e,kr,kc,NN_row,NN_col,i,icol,irow;      /* index_list is an array of linked lists. Each entry is a row (DOF) and contains the indices to the non-zero columns */
36        index_t color;
37    if (elements!=NULL) {      dim_t e, kr, kc, NN_row, NN_col, icol, irow, NN;
38      dim_t NN=elements->ReferenceElement->Type->numNodes;      if (elements != NULL)
39      index_t id[NN],*row_node,*col_node;      {
40      for (i=0;i<NN;i++) id[i]=i;      NN = elements->numNodes;
41      if (reduce_col_order) {      NN_col = (elements->numShapes);
42         col_node=elements->ReferenceElement->Type->linearNodes;      NN_row = (elements->numShapes);
43         NN_col=elements->LinearReferenceElement->Type->numNodes;  
44      } else {      for (color = elements->minColor; color <= elements->maxColor; color++)
45         col_node=id;      {
46         NN_col=elements->ReferenceElement->Type->numNodes;  #pragma omp for private(e,irow,kr,kc,icol) schedule(static)
47      }          for (e = 0; e < elements->numElements; e++)
48      if (reduce_row_order) {          {
49         row_node=elements->ReferenceElement->Type->linearNodes;          if (elements->Color[e] == color)
50         NN_row=elements->LinearReferenceElement->Type->numNodes;          {
51      } else {              for (kr = 0; kr < NN_row; kr++)
52         row_node=id;              {
53         NN_row=elements->ReferenceElement->Type->numNodes;              irow = row_map[elements->Nodes[INDEX2(kr, e, NN)]];
54      }              for (kc = 0; kc < NN_col; kc++)
55      for (color=elements->minColor;color<=elements->maxColor;color++) {              {
56          #pragma omp for private(e,irow,kr,kc,icol) schedule(static)                  icol = col_map[elements->Nodes[INDEX2(kc, e, NN)]];
57          for (e=0;e<elements->numElements;e++) {                  Dudley_IndexList_insertIndex(&(index_list[irow]), icol);
58              if (elements->Color[e]==color) {              }
59                  for (kr=0;kr<NN_row;kr++) {              }
60                    irow=row_Label[elements->Nodes[INDEX2(row_node[kr],e,NN)]];          }
61                    for (kc=0;kc<NN_col;kc++) {          }
62                         icol=col_Label[elements->Nodes[INDEX2(col_node[kc],e,NN)]];      }
63                         Finley_IndexList_insertIndex(&(index_list[irow]),icol);      }
64                    }      return;
65                  }  }
66              }  
67          }  void Dudley_IndexList_insertElementsWithRowRange(Dudley_IndexList * index_list, index_t firstRow, index_t lastRow,
68        }                           Dudley_ElementFile * elements, index_t * row_map, index_t * col_map)
69    }  {
70    return;  /* this does not resolve macro elements */
71  }      index_t color;
72        dim_t e, kr, kc, icol, irow, NN;
73  /* inserts row index row into the Finley_IndexList in if it does not exist */      if (elements != NULL)
74        {
75  void Finley_IndexList_insertIndex(Finley_IndexList* in, index_t index) {      NN = elements->numNodes;
76    dim_t i;      for (color = elements->minColor; color <= elements->maxColor; color++)
77    /* is index in in? */      {
78    for (i=0;i<in->n;i++) {  #pragma omp for private(e,irow,kr,kc,icol) schedule(static)
79      if (in->index[i]==index)  return;          for (e = 0; e < elements->numElements; e++)
80    }          {
81    /* index could not be found */          if (elements->Color[e] == color)
82    if (in->n==INDEXLIST_LENGTH) {          {
83       /* if in->index is full check the extension */              for (kr = 0; kr < NN; kr++)
84       if (in->extension==NULL) {              {
85          in->extension=TMPMEMALLOC(1,Finley_IndexList);              irow = row_map[elements->Nodes[INDEX2(kr, e, NN)]];
86          if (Finley_checkPtr(in->extension)) return;              if ((firstRow <= irow) && (irow < lastRow))
87          in->extension->n=0;              {
88          in->extension->extension=NULL;                  irow -= firstRow;
89       }                  for (kc = 0; kc < NN; kc++)
90       Finley_IndexList_insertIndex(in->extension,index);                  {
91    } else {                  icol = col_map[elements->Nodes[INDEX2(kc, e, NN)]];
92       /* insert index into in->index*/                  Dudley_IndexList_insertIndex(&(index_list[irow]), icol);
93       in->index[in->n]=index;                  }
94       in->n++;              }
95    }              }
96  }          }
97            }
98  /* counts the number of row indices in the Finley_IndexList in */      }
99        }
100  dim_t Finley_IndexList_count(Finley_IndexList* in) {  }
101    if (in==NULL) {  
102       return 0;  void Dudley_IndexList_insertElementsWithRowRangeNoMainDiagonal(Dudley_IndexList * index_list, index_t firstRow,
103    } else {                                     index_t lastRow, Dudley_ElementFile * elements,
104       return (in->n)+Finley_IndexList_count(in->extension);                                     index_t * row_map, index_t * col_map)
105    }  {
106  }      /* this does not resolve macro elements */
107        index_t color;
108  /* count the number of row indices in the Finley_IndexList in */      dim_t e, kr, kc, icol, irow, NN, irow_loc;
109        if (elements != NULL)
110  void Finley_IndexList_toArray(Finley_IndexList* in, index_t* array) {      {
111    dim_t i;      NN = elements->numNodes;
112    if (in!=NULL) {      for (color = elements->minColor; color <= elements->maxColor; color++)
113      for (i=0;i<in->n;i++) array[i]=in->index[i]+INDEX_OFFSET;      {
114      Finley_IndexList_toArray(in->extension,&(array[in->n]));  #pragma omp for private(e,irow,kr,kc,icol,irow_loc) schedule(static)
115    }          for (e = 0; e < elements->numElements; e++)
116  }          {
117            if (elements->Color[e] == color)
118  /* deallocates the Finley_IndexList in by recursive calls */          {
119                for (kr = 0; kr < NN; kr++)
120  void Finley_IndexList_free(Finley_IndexList* in) {              {
121    if (in!=NULL) {              irow = row_map[elements->Nodes[INDEX2(kr, e, NN)]];
122      Finley_IndexList_free(in->extension);              if ((firstRow <= irow) && (irow < lastRow))
123      TMPMEMFREE(in);              {
124    }                  irow_loc = irow - firstRow;
125  }                  for (kc = 0; kc < NN; kc++)
126                    {
127  /*                  icol = col_map[elements->Nodes[INDEX2(kc, e, NN)]];
128   * $Log$                  if (icol != irow)
129   * Revision 1.6  2005/09/15 03:44:22  jgs                      Dudley_IndexList_insertIndex(&(index_list[irow_loc]), icol);
130   * Merge of development branch dev-02 back to main trunk on 2005-09-15                  }
131   *              }
132   * Revision 1.5.2.1  2005/09/07 06:26:18  gross              }
133   * the solver from finley are put into the standalone package paso now          }
134   *          }
135   * Revision 1.5  2005/07/08 04:07:51  jgs      }
136   * Merge of development branch back to main trunk on 2005-07-08      }
137   *  }
138   * Revision 1.4  2004/12/15 07:08:32  jgs  
139   * *** empty log message ***  /* inserts row index row into the Dudley_IndexList in if it does not exist */
140   * Revision 1.1.1.1.2.3  2005/06/29 02:34:50  gross  
141   * some changes towards 64 integers in finley  void Dudley_IndexList_insertIndex(Dudley_IndexList * in, index_t index)
142   *  {
143   * Revision 1.1.1.1.2.2  2004/11/24 01:37:13  gross      dim_t i;
144   * some changes dealing with the integer overflow in memory allocation. Finley solves 4M unknowns now      /* is index in in? */
145   *      for (i = 0; i < in->n; i++)
146   *      {
147   *      if (in->index[i] == index)
148   */          return;
149        }
150        /* index could not be found */
151        if (in->n == INDEXLIST_LENGTH)
152        {
153        /* if in->index is full check the extension */
154        if (in->extension == NULL)
155        {
156            in->extension = TMPMEMALLOC(1, Dudley_IndexList);
157            if (Dudley_checkPtr(in->extension))
158            return;
159            in->extension->n = 0;
160            in->extension->extension = NULL;
161        }
162        Dudley_IndexList_insertIndex(in->extension, index);
163        }
164        else
165        {
166        /* insert index into in->index */
167        in->index[in->n] = index;
168        in->n++;
169        }
170    }
171    
172    /* counts the number of row indices in the Dudley_IndexList in */
173    
174    dim_t Dudley_IndexList_count(Dudley_IndexList * in, index_t range_min, index_t range_max)
175    {
176        dim_t i;
177        dim_t out = 0;
178        register index_t itmp;
179        if (in == NULL)
180        {
181        return 0;
182        }
183        else
184        {
185        for (i = 0; i < in->n; i++)
186        {
187            itmp = in->index[i];
188            if ((itmp >= range_min) && (range_max > itmp))
189            ++out;
190        }
191        return out + Dudley_IndexList_count(in->extension, range_min, range_max);
192        }
193    }
194    
195    /* count the number of row indices in the Dudley_IndexList in */
196    
197    void Dudley_IndexList_toArray(Dudley_IndexList * in, index_t * array, index_t range_min, index_t range_max,
198                      index_t index_offset)
199    {
200        dim_t i, ptr;
201        register index_t itmp;
202        if (in != NULL)
203        {
204        ptr = 0;
205        for (i = 0; i < in->n; i++)
206        {
207            itmp = in->index[i];
208            if ((itmp >= range_min) && (range_max > itmp))
209            {
210            array[ptr] = itmp + index_offset;
211            ptr++;
212            }
213    
214        }
215        Dudley_IndexList_toArray(in->extension, &(array[ptr]), range_min, range_max, index_offset);
216        }
217    }
218    
219    /* deallocates the Dudley_IndexList in by recursive calls */
220    
221    void Dudley_IndexList_free(Dudley_IndexList * in)
222    {
223        if (in != NULL)
224        {
225        Dudley_IndexList_free(in->extension);
226        TMPMEMFREE(in);
227        }
228    }
229    
230    /* creates a Paso_pattern from a range of indices */
231    Paso_Pattern *Dudley_IndexList_createPattern(dim_t n0, dim_t n, Dudley_IndexList * index_list, index_t range_min,
232                             index_t range_max, index_t index_offset)
233    {
234        dim_t *ptr = NULL;
235        register dim_t s, i, itmp;
236        index_t *index = NULL;
237        Paso_Pattern *out = NULL;
238    
239        ptr = MEMALLOC(n + 1 - n0, index_t);
240        if (!Dudley_checkPtr(ptr))
241        {
242        /* get the number of connections per row */
243    #pragma omp parallel for schedule(static) private(i)
244        for (i = n0; i < n; ++i)
245        {
246            ptr[i - n0] = Dudley_IndexList_count(&index_list[i], range_min, range_max);
247        }
248        /* accumulate ptr */
249        s = 0;
250        for (i = n0; i < n; ++i)
251        {
252            itmp = ptr[i - n0];
253            ptr[i - n0] = s;
254            s += itmp;
255        }
256        ptr[n - n0] = s;
257        /* fill index */
258        index = MEMALLOC(ptr[n - n0], index_t);
259        if (!Dudley_checkPtr(index))
260        {
261    #pragma omp parallel for schedule(static)
262            for (i = n0; i < n; ++i)
263            {
264            Dudley_IndexList_toArray(&index_list[i], &index[ptr[i - n0]], range_min, range_max, index_offset);
265            }
266            out = Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT, n - n0, range_max + index_offset, ptr, index);
267        }
268        }
269        if (!Dudley_noError())
270        {
271        MEMFREE(ptr);
272        MEMFREE(index);
273        Paso_Pattern_free(out);
274        }
275        return out;
276    }

Legend:
Removed from v.201  
changed lines
  Added in v.3981

  ViewVC Help
Powered by ViewVC 1.1.26