/[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/esys2/finley/src/finleyC/IndexList.c revision 82 by jgs, Tue Oct 26 06:53:54 2004 UTC trunk/dudley/src/IndexList.c revision 3911 by jfenwick, Thu Jun 14 01:01:03 2012 UTC
# Line 1  Line 1 
 /* $Id$ */  
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  /* Copyrights by ACcESS Australia 2003,2004 */  /* Dudley: Converting an element list into a matrix shape     */
 /* Author: gross@access.edu.au */  
17    
18  /**************************************************************/  /**************************************************************/
19    
 #include "Finley.h"  
 #include "ElementFile.h"  
 #include "System.h"  
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                                         int reduce_row_order, int packed_row_block_size, maybelong* row_Label,                       bool_t reduce_row_order, index_t * row_map,
31                                         int reduce_col_order, int packed_col_block_size, maybelong* col_Label,                       bool_t reduce_col_order, index_t * col_map)
32                                         int symmetric, Finley_SystemMatrixType matType) {  {
33    maybelong e,kr,ir,kc,ic,NN_row,NN_col,i;      /* index_list is an array of linked lists. Each entry is a row (DOF) and contains the indices to the non-zero columns */
34    Finley_IndexList * tmp_list;      index_t color;
35    maybelong jr,jc,icol,irow,color;      dim_t e, kr, kc, NN_row, NN_col, icol, irow, NN;
36        if (elements != NULL)
37        {
38    if (elements!=NULL) {      NN = elements->numNodes;
39      maybelong NN=elements->ReferenceElement->Type->numNodes;      NN_col = (elements->numShapes);
40      maybelong id[NN],*row_node,*col_node;      NN_row = (elements->numShapes);
41      for (i=0;i<NN;i++) id[i]=i;  
42      if (reduce_col_order) {      for (color = elements->minColor; color <= elements->maxColor; color++)
43         col_node=elements->ReferenceElement->Type->linearNodes;      {
44         NN_col=elements->LinearReferenceElement->Type->numNodes;  #pragma omp for private(e,irow,kr,kc,icol) schedule(static)
45      } else {          for (e = 0; e < elements->numElements; e++)
46         col_node=id;          {
47         NN_col=elements->ReferenceElement->Type->numNodes;          if (elements->Color[e] == color)
48      }          {
49      if (reduce_row_order) {              for (kr = 0; kr < NN_row; kr++)
50         row_node=elements->ReferenceElement->Type->linearNodes;              {
51         NN_row=elements->LinearReferenceElement->Type->numNodes;              irow = row_map[elements->Nodes[INDEX2(kr, e, NN)]];
52      } else {              for (kc = 0; kc < NN_col; kc++)
53         row_node=id;              {
54         NN_row=elements->ReferenceElement->Type->numNodes;                  icol = col_map[elements->Nodes[INDEX2(kc, e, NN)]];
55      }                  Dudley_IndexList_insertIndex(&(index_list[irow]), icol);
56                }
57      #pragma omp parallel private(color)              }
58      {          }
59         switch(matType) {          }
60         case CSR:      }
61           if (symmetric) {      }
62              for (color=0;color<elements->numColors;color++) {      return;
63                 #pragma omp for private(e,kr,jr,kc,jc,ir,irow,ic,icol,tmp_list) schedule(static)  }
64                 for (e=0;e<elements->numElements;e++) {  
65                    if (elements->Color[e]==color) {  void Dudley_IndexList_insertElementsWithRowRange(Dudley_IndexList * index_list, index_t firstRow, index_t lastRow,
66                       for (kr=0;kr<NN_row;kr++) {                           Dudley_ElementFile * elements, index_t * row_map, index_t * col_map)
67                          jr=row_Label[elements->Nodes[INDEX2(row_node[kr],e,NN)]];  {
68                          for (ir=0;ir<packed_row_block_size;ir++) {  /* this does not resolve macro elements */
69                             irow=ir+packed_row_block_size*jr;      index_t color;
70                             tmp_list=&(index_list[irow]);      dim_t e, kr, kc, icol, irow, NN;
71                             for (kc=0;kc<NN_col;kc++) {      if (elements != NULL)
72                                jc=col_Label[elements->Nodes[INDEX2(col_node[kc],e,NN)]];      {
73                                for (ic=0;ic<packed_col_block_size;ic++) {      NN = elements->numNodes;
74                                   icol=ic+packed_col_block_size*jc;      for (color = elements->minColor; color <= elements->maxColor; color++)
75                                   if (irow<=icol) Finley_IndexList_insertIndex(tmp_list,icol);      {
76                                }  #pragma omp for private(e,irow,kr,kc,icol) schedule(static)
77                             }          for (e = 0; e < elements->numElements; e++)
78                          }          {
79                       }          if (elements->Color[e] == color)
80                    }          {
81                 }              for (kr = 0; kr < NN; kr++)
82              }              {
83           } else {              irow = row_map[elements->Nodes[INDEX2(kr, e, NN)]];
84              for (color=0;color<elements->numColors;color++) {              if ((firstRow <= irow) && (irow < lastRow))
85                 #pragma omp for private(e,kr,jr,kc,jc,ir,irow,ic,icol,tmp_list) schedule(static)              {
86                 for (e=0;e<elements->numElements;e++) {                  irow -= firstRow;
87                    if (elements->Color[e]==color) {                  for (kc = 0; kc < NN; kc++)
88                       for (kr=0;kr<NN_row;kr++) {                  {
89                          jr=row_Label[elements->Nodes[INDEX2(row_node[kr],e,NN)]];                  icol = col_map[elements->Nodes[INDEX2(kc, e, NN)]];
90                          for (ir=0;ir<packed_row_block_size;ir++) {                  Dudley_IndexList_insertIndex(&(index_list[irow]), icol);
91                             irow=ir+packed_row_block_size*jr;                  }
92                             tmp_list=&(index_list[irow]);              }
93                             for (kc=0;kc<NN_col;kc++) {              }
94                                jc=col_Label[elements->Nodes[INDEX2(col_node[kc],e,NN)]];          }
95                                for (ic=0;ic<packed_col_block_size;ic++) {                                  icol=ic+packed_col_block_size*jc;          }
96                                   Finley_IndexList_insertIndex(tmp_list,icol);      }
97                                }      }
98                             }  }
99                          }  
100                       }  void Dudley_IndexList_insertElementsWithRowRangeNoMainDiagonal(Dudley_IndexList * index_list, index_t firstRow,
101                    }                                     index_t lastRow, Dudley_ElementFile * elements,
102                 }                                     index_t * row_map, index_t * col_map)
103              }  {
104           } /* if else */      /* this does not resolve macro elements */
105           break;      index_t color;
106         case CSC:      dim_t e, kr, kc, icol, irow, NN, irow_loc;
107           if (symmetric) {      if (elements != NULL)
108              for (color=0;color<elements->numColors;color++) {      {
109                 #pragma omp for private (e,kr,jr,kc,jc,ir,irow,ic,icol,tmp_list) schedule(static)      NN = elements->numNodes;
110                 for (e=0;e<elements->numElements;e++) {      for (color = elements->minColor; color <= elements->maxColor; color++)
111                    if (elements->Color[e]==color) {      {
112                       for (kc=0;kc<NN_col;kc++) {  #pragma omp for private(e,irow,kr,kc,icol,irow_loc) schedule(static)
113                          jc=col_Label[elements->Nodes[INDEX2(col_node[kc],e,NN)]];          for (e = 0; e < elements->numElements; e++)
114                          for (ic=0;ic<packed_col_block_size;ic++) {          {
115                             icol=ic+packed_col_block_size*jc;          if (elements->Color[e] == color)
116                             tmp_list=&(index_list[icol]);          {
117                             for (kr=0;kr<NN_row;kr++) {              for (kr = 0; kr < NN; kr++)
118                                jr=row_Label[elements->Nodes[INDEX2(row_node[kr],e,NN)]];              {
119                                for (ir=0;ir<packed_row_block_size;ir++) {              irow = row_map[elements->Nodes[INDEX2(kr, e, NN)]];
120                                   irow=ir+packed_row_block_size*jr;              if ((firstRow <= irow) && (irow < lastRow))
121                                   if (irow<=icol) Finley_IndexList_insertIndex(tmp_list,irow);              {
122                                }                  irow_loc = irow - firstRow;
123                             }                  for (kc = 0; kc < NN; kc++)
124                          }                  {
125                       }                  icol = col_map[elements->Nodes[INDEX2(kc, e, NN)]];
126                    }                  if (icol != irow)
127                 }                      Dudley_IndexList_insertIndex(&(index_list[irow_loc]), icol);
128              }                  }
129           } else {              }
130             for (color=0;color<elements->numColors;color++) {              }
131                 #pragma omp for private (e,kr,jr,kc,jc,ir,irow,ic,icol,tmp_list) schedule(static)          }
132                 for (e=0;e<elements->numElements;e++) {          }
133                    if (elements->Color[e]==color) {      }
134                       for (kc=0;kc<NN_col;kc++) {      }
135                          jc=col_Label[elements->Nodes[INDEX2(col_node[kc],e,NN)]];  }
136                          for (ic=0;ic<packed_col_block_size;ic++) {  
137                             icol=ic+packed_col_block_size*jc;  /* inserts row index row into the Dudley_IndexList in if it does not exist */
138                             tmp_list=&(index_list[icol]);  
139                             for (kr=0;kr<NN_row;kr++) {  void Dudley_IndexList_insertIndex(Dudley_IndexList * in, index_t index)
140                                jr=row_Label[elements->Nodes[INDEX2(row_node[kr],e,NN)]];  {
141                                for (ir=0;ir<packed_row_block_size;ir++) {      dim_t i;
142                                   irow=ir+packed_row_block_size*jr;      /* is index in in? */
143                                   Finley_IndexList_insertIndex(tmp_list,irow);      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 else */      /* if in->index is full check the extension */
152         } /* switch matType */      if (in->extension == NULL)
153      }      {
154    }          in->extension = TMPMEMALLOC(1, Dudley_IndexList);
155    return;          if (Dudley_checkPtr(in->extension))
156  }          return;
157            in->extension->n = 0;
158  /* inserts row index row into the Finley_IndexList in if it does not exist */          in->extension->extension = NULL;
159        }
160  void Finley_IndexList_insertIndex(Finley_IndexList* in, maybelong index) {      Dudley_IndexList_insertIndex(in->extension, index);
161    int i;      }
162    /* is index in in? */      else
163    for (i=0;i<in->n;i++) {      {
164      if (in->index[i]==index) return;      /* insert index into in->index */
165    }      in->index[in->n] = index;
166    /* index could not be found */      in->n++;
167    if (in->n==INDEXLIST_LENGTH) {      }
168       /* if in->index is full check the extension */  }
169       if (in->extension==NULL) {  
170          in->extension=(Finley_IndexList*) TMPMEMALLOC(sizeof(Finley_IndexList));  /* counts the number of row indices in the Dudley_IndexList in */
171          if (Finley_checkPtr(in->extension)) return;  
172          in->extension->n=0;  dim_t Dudley_IndexList_count(Dudley_IndexList * in, index_t range_min, index_t range_max)
173          in->extension->extension=NULL;  {
174       }      dim_t i;
175       Finley_IndexList_insertIndex(in->extension,index);      dim_t out = 0;
176    } else {      register index_t itmp;
177       /* insert index into in->index*/      if (in == NULL)
178       in->index[in->n]=index;      {
179       in->n++;      return 0;
180    }      }
181  }      else
182        {
183  /* counts the number of row indices in the Finley_IndexList in */      for (i = 0; i < in->n; i++)
184        {
185  int Finley_IndexList_count(Finley_IndexList* in) {          itmp = in->index[i];
186    if (in==NULL) {          if ((itmp >= range_min) && (range_max > itmp))
187       return 0;          ++out;
188    } else {      }
189       return (in->n)+Finley_IndexList_count(in->extension);      return out + Dudley_IndexList_count(in->extension, range_min, range_max);
190    }      }
191  }  }
192    
193  /* count the number of row indices in the Finley_IndexList in */  /* count the number of row indices in the Dudley_IndexList in */
194    
195  void Finley_IndexList_toArray(Finley_IndexList* in, maybelong* array) {  void Dudley_IndexList_toArray(Dudley_IndexList * in, index_t * array, index_t range_min, index_t range_max,
196    int i;                    index_t index_offset)
197    if (in!=NULL) {  {
198      for (i=0;i<in->n;i++) array[i]=in->index[i]+INDEX_OFFSET;      dim_t i, ptr;
199      Finley_IndexList_toArray(in->extension,&(array[in->n]));      register index_t itmp;
200    }      if (in != NULL)
201  }      {
202        ptr = 0;
203  /* deallocates the Finley_IndexList in by recursive calls */      for (i = 0; i < in->n; i++)
204        {
205  void Finley_IndexList_free(Finley_IndexList* in) {          itmp = in->index[i];
206    if (in!=NULL) {          if ((itmp >= range_min) && (range_max > itmp))
207      Finley_IndexList_free(in->extension);          {
208      TMPMEMFREE(in);          array[ptr] = itmp + index_offset;
209    }          ptr++;
210  }          }
211    
212  /*      }
213   * $Log$      Dudley_IndexList_toArray(in->extension, &(array[ptr]), range_min, range_max, index_offset);
214   * Revision 1.1  2004/10/26 06:53:57  jgs      }
215   * Initial revision  }
216   *  
217   * Revision 1.1.2.2  2004/10/26 06:36:39  jgs  /* deallocates the Dudley_IndexList in by recursive calls */
218   * committing Lutz's changes to branch jgs  
219   *  void Dudley_IndexList_free(Dudley_IndexList * in)
220   * Revision 1.2  2004/10/13 01:53:42  gross  {
221   * bug in CSC assembling fixed      if (in != NULL)
222   *      {
223   * Revision 1.1  2004/07/02 04:21:13  gross      Dudley_IndexList_free(in->extension);
224   * Finley C code has been included      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.82  
changed lines
  Added in v.3911

  ViewVC Help
Powered by ViewVC 1.1.26