# Diff of /trunk/finley/src/IndexList.c

trunk/esys2/finley/src/finleyC/IndexList.c revision 123 by jgs, Fri Jul 8 04:08:13 2005 UTC trunk/finley/src/IndexList.c revision 2748 by gross, Tue Nov 17 07:32:59 2009 UTC
# Line 1  Line 1
/* \$Id\$ */
1
2  /**************************************************************/  /*******************************************************
3    *
4    * Copyright (c) 2003-2009 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
/* Finley: Converting an element list into a matrix shape     */
14
15  /**************************************************************/  /**************************************************************/
16
17  /* Copyrights by ACcESS Australia 2003,2004 */  /* Finley: Converting an element list into a matrix shape     */
/* Author: gross@access.edu.au */
18
19  /**************************************************************/  /**************************************************************/
20
#include "Finley.h"
#include "ElementFile.h"
#include "System.h"
21  #include "IndexList.h"  #include "IndexList.h"
22
23    /* Translate from distributed/local array indices to global indices */
24
25  /**************************************************************/  /**************************************************************/
26  /* inserts the contributions from the element matrices of elements  /* inserts the contributions from the element matrices of elements
27     into the row index col. If symmetric is set, only the upper     into the row index col. If symmetric is set, only the upper
28     triangle of the matrix is stored. */     triangle of the matrix is stored. */
29
30  void Finley_IndexList_insertElements(Finley_IndexList* index_list, Finley_ElementFile* elements,  void Finley_IndexList_insertElements(Finley_IndexList* index_list, Finley_ElementFile* elements,
31                                         bool_t reduce_row_order, index_t* row_Label,                                         bool_t reduce_row_order, index_t* row_map,
32                                         bool_t reduce_col_order, index_t* col_Label) {                                         bool_t reduce_col_order, index_t* col_map) {
33      /* 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;    index_t color;
35    dim_t e,kr,kc,NN_row,NN_col,i,icol,irow;    Finley_ReferenceElement*refElement;
36      dim_t e, kr, kc, NN_row, NN_col, icol, irow, NN, *row_node=NULL, *col_node=NULL, isub, numSub;
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;      refElement= Finley_ReferenceElementSet_borrowReferenceElement(elements->referenceElementSet, FALSE);
40      for (i=0;i<NN;i++) id[i]=i;
41      if (reduce_col_order) {      if (reduce_col_order) {
42         col_node=elements->ReferenceElement->Type->linearNodes;            numSub=1;
43         NN_col=elements->LinearReferenceElement->Type->numNodes;            col_node=refElement->Type->linearNodes;
44              NN_col=(refElement->LinearBasisFunctions->Type->numShapes) * (refElement->Type->numSides) ;
45      } else {      } else {
46         col_node=id;            numSub=refElement->Type->numSubElements;
47         NN_col=elements->ReferenceElement->Type->numNodes;            col_node=refElement->Type->subElementNodes;
48              NN_col=(refElement->BasisFunctions->Type->numShapes) * (refElement->Type->numSides) ;
49      }      }
50
51      if (reduce_row_order) {      if (reduce_row_order) {
52         row_node=elements->ReferenceElement->Type->linearNodes;            numSub=1;
53         NN_row=elements->LinearReferenceElement->Type->numNodes;            row_node=refElement->Type->linearNodes;
54              NN_row=(refElement->LinearBasisFunctions->Type->numShapes) * (refElement->Type->numSides) ;
55      } else {      } else {
56         row_node=id;            numSub=refElement->Type->numSubElements;
57         NN_row=elements->ReferenceElement->Type->numNodes;            row_node=refElement->Type->subElementNodes;
58              NN_row=(refElement->BasisFunctions->Type->numShapes) * (refElement->Type->numSides) ;
59        }
60
61        for (color=elements->minColor;color<=elements->maxColor;color++) {
62               #pragma omp for private(e,irow,kr,kc,icol,isub) schedule(static)
63               for (e=0;e<elements->numElements;e++) {
64                   if (elements->Color[e]==color) {
65                       for (isub=0;isub<numSub; isub++) {
66                           for (kr=0;kr<NN_row;kr++) {
67                               irow=row_map[elements->Nodes[INDEX2(row_node[INDEX2(kr,isub,NN_row)],e,NN)]];
68                               for (kc=0;kc<NN_col;kc++) {
69                                   icol=col_map[elements->Nodes[INDEX2(col_node[INDEX2(kc,isub,NN_col)],e,NN)]];
70                                   Finley_IndexList_insertIndex(&(index_list[irow]),icol);
71                               }
72                           }
73                       }
74                   }
75               }
76           }
77      }
78      return;
79    }
80
81
82    void Finley_IndexList_insertElementsWithRowRange(Finley_IndexList* index_list, index_t firstRow, index_t lastRow,
83                                                     Finley_ElementFile* elements, index_t* row_map, index_t* col_map)
84    {
85    /* this does not resolve macro elements */
86        index_t color;
87      dim_t e,kr,kc,icol,irow, NN;
88      if (elements!=NULL) {
89        NN=elements->numNodes;
90        for (color=elements->minColor;color<=elements->maxColor;color++) {
91               #pragma omp for private(e,irow,kr,kc,icol) schedule(static)
92               for (e=0;e<elements->numElements;e++) {
93                   if (elements->Color[e]==color) {
94                       for (kr=0;kr<NN;kr++) {
95                         irow=row_map[elements->Nodes[INDEX2(kr,e,NN)]];
96                         if ((firstRow<=irow) && (irow < lastRow)) {
97                              irow-=firstRow;
98                              for (kc=0;kc<NN;kc++) {
99                                  icol=col_map[elements->Nodes[INDEX2(kc,e,NN)]];
100                                  Finley_IndexList_insertIndex(&(index_list[irow]),icol);
101                              }
102                          }
103                      }
104                   }
105               }
106      }      }
107      }
108    }
109    void Finley_IndexList_insertElementsWithRowRangeNoMainDiagonal(Finley_IndexList* index_list, index_t firstRow, index_t lastRow,
110                                                                  Finley_ElementFile* elements, index_t* row_map, index_t* col_map)
111    {
112      /* this does not resolve macro elements */
113      index_t color;
114      dim_t e,kr,kc,icol,irow, NN,irow_loc;
115      if (elements!=NULL) {
116        NN=elements->numNodes;
117      for (color=elements->minColor;color<=elements->maxColor;color++) {      for (color=elements->minColor;color<=elements->maxColor;color++) {
118          #pragma omp for private(e,irow,kr,kc,icol) schedule(static)             #pragma omp for private(e,irow,kr,kc,icol,irow_loc) schedule(static)
119          for (e=0;e<elements->numElements;e++) {             for (e=0;e<elements->numElements;e++) {
120              if (elements->Color[e]==color) {                 if (elements->Color[e]==color) {
121                  for (kr=0;kr<NN_row;kr++) {                     for (kr=0;kr<NN;kr++) {
122                    irow=row_Label[elements->Nodes[INDEX2(row_node[kr],e,NN)]];                       irow=row_map[elements->Nodes[INDEX2(kr,e,NN)]];
123                    for (kc=0;kc<NN_col;kc++) {                       if ((firstRow<=irow) && (irow < lastRow)) {
124                         icol=col_Label[elements->Nodes[INDEX2(col_node[kc],e,NN)]];                            irow_loc=irow-firstRow;
125                         Finley_IndexList_insertIndex(&(index_list[irow]),icol);                            for (kc=0;kc<NN;kc++) {
126                                  icol=col_map[elements->Nodes[INDEX2(kc,e,NN)]];
127                                  if (icol != irow) Finley_IndexList_insertIndex(&(index_list[irow_loc]),icol);
128                              }
129                          }
130                    }                    }
131                  }                 }
132              }             }
133          }      }
}
134    }    }
return;
135  }  }
136
137  /* 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 90  void Finley_IndexList_insertIndex(Finley Line 161  void Finley_IndexList_insertIndex(Finley
161
162  /* counts the number of row indices in the Finley_IndexList in */  /* counts the number of row indices in the Finley_IndexList in */
163
164  dim_t Finley_IndexList_count(Finley_IndexList* in) {  dim_t Finley_IndexList_count(Finley_IndexList* in, index_t range_min,index_t range_max) {
165      dim_t i;
166      dim_t out=0;
167      register index_t itmp;
168    if (in==NULL) {    if (in==NULL) {
169       return 0;       return 0;
170    } else {    } else {
171       return (in->n)+Finley_IndexList_count(in->extension);      for (i=0;i<in->n;i++) {
172              itmp=in->index[i];
173              if ((itmp>=range_min) && (range_max>itmp)) ++out;
174        }
175         return out+Finley_IndexList_count(in->extension, range_min,range_max);
176    }    }
177  }  }
178
179  /* count the number of row indices in the Finley_IndexList in */  /* count the number of row indices in the Finley_IndexList in */
180
181  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) {
182    dim_t i;    dim_t i, ptr;
183      register index_t itmp;
184    if (in!=NULL) {    if (in!=NULL) {
185      for (i=0;i<in->n;i++) array[i]=in->index[i]+INDEX_OFFSET;      ptr=0;
186      Finley_IndexList_toArray(in->extension,&(array[in->n]));      for (i=0;i<in->n;i++) {
187              itmp=in->index[i];
188              if ((itmp>=range_min) && (range_max>itmp)) {
189                 array[ptr]=itmp+index_offset;
190                 ptr++;
191              }
192
193        }
194        Finley_IndexList_toArray(in->extension,&(array[ptr]), range_min, range_max, index_offset);
195    }    }
196  }  }
197
# Line 117  void Finley_IndexList_free(Finley_IndexL Line 204  void Finley_IndexList_free(Finley_IndexL
204    }    }
205  }  }
206
207  /*  /* creates a Paso_pattern from a range of indices */
208   * \$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)
209   * Revision 1.5  2005/07/08 04:07:51  jgs  {
210   * Merge of development branch back to main trunk on 2005-07-08     dim_t *ptr=NULL;
211   *     register dim_t s,i,itmp;
212   * Revision 1.4  2004/12/15 07:08:32  jgs     index_t *index=NULL;
213   * *** empty log message ***     Paso_Pattern* out=NULL;
214   * Revision 1.1.1.1.2.3  2005/06/29 02:34:50  gross
215   * some changes towards 64 integers in finley     ptr=MEMALLOC(n+1-n0,index_t);
216   *     if (! Finley_checkPtr(ptr) ) {
217   * Revision 1.1.1.1.2.2  2004/11/24 01:37:13  gross         /* get the number of connections per row */
218   * some changes dealing with the integer overflow in memory allocation. Finley solves 4M unknowns now         #pragma omp parallel for schedule(static) private(i)
219   *         for(i=n0;i<n;++i) {
220   *                ptr[i-n0]=Finley_IndexList_count(&index_list[i],range_min,range_max);
221   *         }
222   */         /* accumulate ptr */
223           s=0;
224           for(i=n0;i<n;++i) {
225                   itmp=ptr[i-n0];
226                   ptr[i-n0]=s;
227                   s+=itmp;
228           }
229           ptr[n-n0]=s;
230           /* fill index */
231           index=MEMALLOC(ptr[n-n0],index_t);
232           if (! Finley_checkPtr(index)) {
233                  #pragma omp parallel for schedule(static)
234                  for(i=n0;i<n;++i) {
235                      Finley_IndexList_toArray(&index_list[i],&index[ptr[i-n0]],range_min,range_max,index_offset);
236                  }
237                  out=Paso_Pattern_alloc(PATTERN_FORMAT_DEFAULT,n-n0,range_max+index_offset,ptr,index);
238           }
239      }
240      if (! Finley_noError()) {
241            MEMFREE(ptr);
242            MEMFREE(index);
243            Paso_Pattern_free(out);
244      }
245      return out;
246    }

Legend:
 Removed from v.123 changed lines Added in v.2748

 ViewVC Help Powered by ViewVC 1.1.26