/[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 471 by jgs, Fri Jan 27 01:33:02 2006 UTC trunk/dudley/src/IndexList.c revision 3911 by jfenwick, Thu Jun 14 01:01:03 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  /**************************************************************/  /*******************************************************
3    *
4  /* Finley: Converting an element list into a matrix shape     */  * Copyright (c) 2003-2012 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    
20  #include "IndexList.h"  #include "IndexList.h"
21    
22    /* Translate from distributed/local array indices to global indices */
23    
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_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_t color;  {
33    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 */
34        index_t color;
35    if (elements!=NULL) {      dim_t e, kr, kc, NN_row, NN_col, icol, irow, NN;
36      dim_t NN=elements->ReferenceElement->Type->numNodes;      if (elements != NULL)
37      index_t id[NN],*row_node,*col_node;      {
38      for (i=0;i<NN;i++) id[i]=i;      NN = elements->numNodes;
39      if (reduce_col_order) {      NN_col = (elements->numShapes);
40         col_node=elements->ReferenceElement->Type->linearNodes;      NN_row = (elements->numShapes);
41         NN_col=elements->LinearReferenceElement->Type->numNodes;  
42      } else {      for (color = elements->minColor; color <= elements->maxColor; color++)
43         col_node=id;      {
44         NN_col=elements->ReferenceElement->Type->numNodes;  #pragma omp for private(e,irow,kr,kc,icol) schedule(static)
45      }          for (e = 0; e < elements->numElements; e++)
46      if (reduce_row_order) {          {
47         row_node=elements->ReferenceElement->Type->linearNodes;          if (elements->Color[e] == color)
48         NN_row=elements->LinearReferenceElement->Type->numNodes;          {
49      } else {              for (kr = 0; kr < NN_row; kr++)
50         row_node=id;              {
51         NN_row=elements->ReferenceElement->Type->numNodes;              irow = row_map[elements->Nodes[INDEX2(kr, e, NN)]];
52      }              for (kc = 0; kc < NN_col; kc++)
53      for (color=elements->minColor;color<=elements->maxColor;color++) {              {
54          #pragma omp for private(e,irow,kr,kc,icol) schedule(static)                  icol = col_map[elements->Nodes[INDEX2(kc, e, NN)]];
55          for (e=0;e<elements->numElements;e++) {                  Dudley_IndexList_insertIndex(&(index_list[irow]), icol);
56              if (elements->Color[e]==color) {              }
57                  for (kr=0;kr<NN_row;kr++) {              }
58                    irow=row_Label[elements->Nodes[INDEX2(row_node[kr],e,NN)]];          }
59                    for (kc=0;kc<NN_col;kc++) {          }
60                         icol=col_Label[elements->Nodes[INDEX2(col_node[kc],e,NN)]];      }
61                         Finley_IndexList_insertIndex(&(index_list[irow]),icol);      }
62                    }      return;
63                  }  }
64              }  
65          }  void Dudley_IndexList_insertElementsWithRowRange(Dudley_IndexList * index_list, index_t firstRow, index_t lastRow,
66        }                           Dudley_ElementFile * elements, index_t * row_map, index_t * col_map)
67    }  {
68    return;  /* this does not resolve macro elements */
69  }      index_t color;
70        dim_t e, kr, kc, icol, irow, NN;
71  /* inserts row index row into the Finley_IndexList in if it does not exist */      if (elements != NULL)
72        {
73  void Finley_IndexList_insertIndex(Finley_IndexList* in, index_t index) {      NN = elements->numNodes;
74    dim_t i;      for (color = elements->minColor; color <= elements->maxColor; color++)
75    /* is index in in? */      {
76    for (i=0;i<in->n;i++) {  #pragma omp for private(e,irow,kr,kc,icol) schedule(static)
77      if (in->index[i]==index)  return;          for (e = 0; e < elements->numElements; e++)
78    }          {
79    /* index could not be found */          if (elements->Color[e] == color)
80    if (in->n==INDEXLIST_LENGTH) {          {
81       /* if in->index is full check the extension */              for (kr = 0; kr < NN; kr++)
82       if (in->extension==NULL) {              {
83          in->extension=TMPMEMALLOC(1,Finley_IndexList);              irow = row_map[elements->Nodes[INDEX2(kr, e, NN)]];
84          if (Finley_checkPtr(in->extension)) return;              if ((firstRow <= irow) && (irow < lastRow))
85          in->extension->n=0;              {
86          in->extension->extension=NULL;                  irow -= firstRow;
87       }                  for (kc = 0; kc < NN; kc++)
88       Finley_IndexList_insertIndex(in->extension,index);                  {
89    } else {                  icol = col_map[elements->Nodes[INDEX2(kc, e, NN)]];
90       /* insert index into in->index*/                  Dudley_IndexList_insertIndex(&(index_list[irow]), icol);
91       in->index[in->n]=index;                  }
92       in->n++;              }
93    }              }
94  }          }
95            }
96  /* counts the number of row indices in the Finley_IndexList in */      }
97        }
98  dim_t Finley_IndexList_count(Finley_IndexList* in) {  }
99    if (in==NULL) {  
100       return 0;  void Dudley_IndexList_insertElementsWithRowRangeNoMainDiagonal(Dudley_IndexList * index_list, index_t firstRow,
101    } else {                                     index_t lastRow, Dudley_ElementFile * elements,
102       return (in->n)+Finley_IndexList_count(in->extension);                                     index_t * row_map, index_t * col_map)
103    }  {
104  }      /* this does not resolve macro elements */
105        index_t color;
106  /* count the number of row indices in the Finley_IndexList in */      dim_t e, kr, kc, icol, irow, NN, irow_loc;
107        if (elements != NULL)
108  void Finley_IndexList_toArray(Finley_IndexList* in, index_t* array) {      {
109    dim_t i;      NN = elements->numNodes;
110    if (in!=NULL) {      for (color = elements->minColor; color <= elements->maxColor; color++)
111      for (i=0;i<in->n;i++) array[i]=in->index[i];      {
112      Finley_IndexList_toArray(in->extension,&(array[in->n]));  #pragma omp for private(e,irow,kr,kc,icol,irow_loc) schedule(static)
113    }          for (e = 0; e < elements->numElements; e++)
114  }          {
115            if (elements->Color[e] == color)
116  /* deallocates the Finley_IndexList in by recursive calls */          {
117                for (kr = 0; kr < NN; kr++)
118  void Finley_IndexList_free(Finley_IndexList* in) {              {
119    if (in!=NULL) {              irow = row_map[elements->Nodes[INDEX2(kr, e, NN)]];
120      Finley_IndexList_free(in->extension);              if ((firstRow <= irow) && (irow < lastRow))
121      TMPMEMFREE(in);              {
122    }                  irow_loc = irow - firstRow;
123  }                  for (kc = 0; kc < NN; kc++)
124                    {
125  /*                  icol = col_map[elements->Nodes[INDEX2(kc, e, NN)]];
126   * $Log$                  if (icol != irow)
127   * Revision 1.6  2005/09/15 03:44:22  jgs                      Dudley_IndexList_insertIndex(&(index_list[irow_loc]), icol);
128   * Merge of development branch dev-02 back to main trunk on 2005-09-15                  }
129   *              }
130   * Revision 1.5.2.1  2005/09/07 06:26:18  gross              }
131   * the solver from finley are put into the standalone package paso now          }
132   *          }
133   * Revision 1.5  2005/07/08 04:07:51  jgs      }
134   * Merge of development branch back to main trunk on 2005-07-08      }
135   *  }
136   * Revision 1.4  2004/12/15 07:08:32  jgs  
137   * *** empty log message ***  /* inserts row index row into the Dudley_IndexList in if it does not exist */
138   * Revision 1.1.1.1.2.3  2005/06/29 02:34:50  gross  
139   * some changes towards 64 integers in finley  void Dudley_IndexList_insertIndex(Dudley_IndexList * in, index_t index)
140   *  {
141   * Revision 1.1.1.1.2.2  2004/11/24 01:37:13  gross      dim_t i;
142   * some changes dealing with the integer overflow in memory allocation. Finley solves 4M unknowns now      /* is index in in? */
143   *      for (i = 0; i < in->n; i++)
144   *      {
145   *      if (in->index[i] == index)
146   */          return;
147        }
148        /* index could not be found */
149        if (in->n == INDEXLIST_LENGTH)
150        {
151        /* if in->index is full check the extension */
152        if (in->extension == NULL)
153        {
154            in->extension = TMPMEMALLOC(1, Dudley_IndexList);
155            if (Dudley_checkPtr(in->extension))
156            return;
157            in->extension->n = 0;
158            in->extension->extension = NULL;
159        }
160        Dudley_IndexList_insertIndex(in->extension, index);
161        }
162        else
163        {
164        /* insert index into in->index */
165        in->index[in->n] = index;
166        in->n++;
167        }
168    }
169    
170    /* counts the number of row indices in the Dudley_IndexList in */
171    
172    dim_t Dudley_IndexList_count(Dudley_IndexList * in, index_t range_min, index_t range_max)
173    {
174        dim_t i;
175        dim_t out = 0;
176        register index_t itmp;
177        if (in == NULL)
178        {
179        return 0;
180        }
181        else
182        {
183        for (i = 0; i < in->n; i++)
184        {
185            itmp = in->index[i];
186            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.471  
changed lines
  Added in v.3911

  ViewVC Help
Powered by ViewVC 1.1.26