# Diff of /branches/domexper/dudley/src/IndexList.c

trunk/esys2/finley/src/finleyC/IndexList.c revision 82 by jgs, Tue Oct 26 06:53:54 2004 UTC branches/domexper/dudley/src/IndexList.c revision 3127 by jfenwick, Wed Sep 1 01:30:28 2010 UTC
# Line 1  Line 1
/* \$Id\$ */
1
2  /**************************************************************/  /*******************************************************
3    *
4    * 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
11    *
12    *******************************************************/
13
/* Finley: Converting an element list into a matrix shape     */
14
15  /**************************************************************/  /**************************************************************/
16
17  /* Copyrights by ACcESS Australia 2003,2004 */  /* Dudley: 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 Dudley_IndexList_insertElements(Dudley_IndexList* index_list, Dudley_ElementFile* elements,
31                                         int reduce_row_order, int packed_row_block_size, maybelong* row_Label,                                         bool_t reduce_row_order, index_t* row_map,
32                                         int reduce_col_order, int packed_col_block_size, maybelong* col_Label,                                         bool_t reduce_col_order, index_t* col_map) {
33                                         int symmetric, Finley_SystemMatrixType matType) {    /* index_list is an array of linked lists. Each entry is a row (DOF) and contains the indices to the non-zero columns */
34    maybelong e,kr,ir,kc,ic,NN_row,NN_col,i;    index_t color;
35    Finley_IndexList * tmp_list;    Dudley_ReferenceElement*refElement;
36    maybelong jr,jc,icol,irow,color;    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      maybelong NN=elements->ReferenceElement->Type->numNodes;      NN=elements->numNodes;
39      maybelong id[NN],*row_node,*col_node;      refElement= Dudley_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);
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);
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);
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) ;
59        }
60      #pragma omp parallel private(color)
61      {      for (color=elements->minColor;color<=elements->maxColor;color++) {
62         switch(matType) {             #pragma omp for private(e,irow,kr,kc,icol,isub) schedule(static)
63         case CSR:             for (e=0;e<elements->numElements;e++) {
64           if (symmetric) {                 if (elements->Color[e]==color) {
65              for (color=0;color<elements->numColors;color++) {                     for (isub=0;isub<numSub; isub++) {
66                 #pragma omp for private(e,kr,jr,kc,jc,ir,irow,ic,icol,tmp_list) schedule(static)                         for (kr=0;kr<NN_row;kr++) {
67                 for (e=0;e<elements->numElements;e++) {                             irow=row_map[elements->Nodes[INDEX2(row_node[INDEX2(kr,isub,NN_row)],e,NN)]];
68                    if (elements->Color[e]==color) {                             for (kc=0;kc<NN_col;kc++) {
69                       for (kr=0;kr<NN_row;kr++) {                                 icol=col_map[elements->Nodes[INDEX2(col_node[INDEX2(kc,isub,NN_col)],e,NN)]];
70                          jr=row_Label[elements->Nodes[INDEX2(row_node[kr],e,NN)]];                                 Dudley_IndexList_insertIndex(&(index_list[irow]),icol);
71                          for (ir=0;ir<packed_row_block_size;ir++) {                             }
72                             irow=ir+packed_row_block_size*jr;                         }
73                             tmp_list=&(index_list[irow]);                     }
for (kc=0;kc<NN_col;kc++) {
jc=col_Label[elements->Nodes[INDEX2(col_node[kc],e,NN)]];
for (ic=0;ic<packed_col_block_size;ic++) {
icol=ic+packed_col_block_size*jc;
if (irow<=icol) Finley_IndexList_insertIndex(tmp_list,icol);
}
}
}
}
}
}
}
} else {
for (color=0;color<elements->numColors;color++) {
#pragma omp for private(e,kr,jr,kc,jc,ir,irow,ic,icol,tmp_list) schedule(static)
for (e=0;e<elements->numElements;e++) {
if (elements->Color[e]==color) {
for (kr=0;kr<NN_row;kr++) {
jr=row_Label[elements->Nodes[INDEX2(row_node[kr],e,NN)]];
for (ir=0;ir<packed_row_block_size;ir++) {
irow=ir+packed_row_block_size*jr;
tmp_list=&(index_list[irow]);
for (kc=0;kc<NN_col;kc++) {
jc=col_Label[elements->Nodes[INDEX2(col_node[kc],e,NN)]];
for (ic=0;ic<packed_col_block_size;ic++) {                                  icol=ic+packed_col_block_size*jc;
Finley_IndexList_insertIndex(tmp_list,icol);
}
}
}
}
}
74                 }                 }
75              }             }
76           } /* if else */         }
77           break;    }
78         case CSC:    return;
79           if (symmetric) {  }
80              for (color=0;color<elements->numColors;color++) {
81                 #pragma omp for private (e,kr,jr,kc,jc,ir,irow,ic,icol,tmp_list) schedule(static)
82                 for (e=0;e<elements->numElements;e++) {  void Dudley_IndexList_insertElementsWithRowRange(Dudley_IndexList* index_list, index_t firstRow, index_t lastRow,
83                    if (elements->Color[e]==color) {                                                   Dudley_ElementFile* elements, index_t* row_map, index_t* col_map)
84                       for (kc=0;kc<NN_col;kc++) {  {
85                          jc=col_Label[elements->Nodes[INDEX2(col_node[kc],e,NN)]];  /* this does not resolve macro elements */
86                          for (ic=0;ic<packed_col_block_size;ic++) {      index_t color;
87                             icol=ic+packed_col_block_size*jc;    dim_t e,kr,kc,icol,irow, NN;
88                             tmp_list=&(index_list[icol]);    if (elements!=NULL) {
89                             for (kr=0;kr<NN_row;kr++) {      NN=elements->numNodes;
90                                jr=row_Label[elements->Nodes[INDEX2(row_node[kr],e,NN)]];      for (color=elements->minColor;color<=elements->maxColor;color++) {
91                                for (ir=0;ir<packed_row_block_size;ir++) {             #pragma omp for private(e,irow,kr,kc,icol) schedule(static)
92                                   irow=ir+packed_row_block_size*jr;             for (e=0;e<elements->numElements;e++) {
93                                   if (irow<=icol) Finley_IndexList_insertIndex(tmp_list,irow);                 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                                  Dudley_IndexList_insertIndex(&(index_list[irow]),icol);
101                              }
102                          }
103                    }                    }
104                 }                 }
105              }             }
106           } else {      }
107             for (color=0;color<elements->numColors;color++) {    }
108                 #pragma omp for private (e,kr,jr,kc,jc,ir,irow,ic,icol,tmp_list) schedule(static)  }
109                 for (e=0;e<elements->numElements;e++) {  void Dudley_IndexList_insertElementsWithRowRangeNoMainDiagonal(Dudley_IndexList* index_list, index_t firstRow, index_t lastRow,
110                    if (elements->Color[e]==color) {                                                                Dudley_ElementFile* elements, index_t* row_map, index_t* col_map)
111                       for (kc=0;kc<NN_col;kc++) {  {
112                          jc=col_Label[elements->Nodes[INDEX2(col_node[kc],e,NN)]];    /* this does not resolve macro elements */
113                          for (ic=0;ic<packed_col_block_size;ic++) {    index_t color;
114                             icol=ic+packed_col_block_size*jc;    dim_t e,kr,kc,icol,irow, NN,irow_loc;
115                             tmp_list=&(index_list[icol]);    if (elements!=NULL) {
116                             for (kr=0;kr<NN_row;kr++) {      NN=elements->numNodes;
117                                jr=row_Label[elements->Nodes[INDEX2(row_node[kr],e,NN)]];      for (color=elements->minColor;color<=elements->maxColor;color++) {
118                                for (ir=0;ir<packed_row_block_size;ir++) {             #pragma omp for private(e,irow,kr,kc,icol,irow_loc) schedule(static)
119                                   irow=ir+packed_row_block_size*jr;             for (e=0;e<elements->numElements;e++) {
120                                   Finley_IndexList_insertIndex(tmp_list,irow);                 if (elements->Color[e]==color) {
121                                }                     for (kr=0;kr<NN;kr++) {
122                             }                       irow=row_map[elements->Nodes[INDEX2(kr,e,NN)]];
123                          }                       if ((firstRow<=irow) && (irow < lastRow)) {
124                       }                            irow_loc=irow-firstRow;
125                              for (kc=0;kc<NN;kc++) {
126                                  icol=col_map[elements->Nodes[INDEX2(kc,e,NN)]];
127                                  if (icol != irow) Dudley_IndexList_insertIndex(&(index_list[irow_loc]),icol);
128                              }
129                          }
130                    }                    }
131                 }                 }
132              }             }
} /* if else */
} /* switch matType */
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 Dudley_IndexList in if it does not exist */
138
139  void Finley_IndexList_insertIndex(Finley_IndexList* in, maybelong index) {  void Dudley_IndexList_insertIndex(Dudley_IndexList* in, index_t index) {
140    int i;    dim_t i;
141    /* is index in in? */    /* is index in in? */
142    for (i=0;i<in->n;i++) {    for (i=0;i<in->n;i++) {
143      if (in->index[i]==index) return;      if (in->index[i]==index)  return;
144    }    }
145    /* index could not be found */    /* index could not be found */
146    if (in->n==INDEXLIST_LENGTH) {    if (in->n==INDEXLIST_LENGTH) {
147       /* if in->index is full check the extension */       /* if in->index is full check the extension */
148       if (in->extension==NULL) {       if (in->extension==NULL) {
149          in->extension=(Finley_IndexList*) TMPMEMALLOC(sizeof(Finley_IndexList));          in->extension=TMPMEMALLOC(1,Dudley_IndexList);
150          if (Finley_checkPtr(in->extension)) return;          if (Dudley_checkPtr(in->extension)) return;
151          in->extension->n=0;          in->extension->n=0;
152          in->extension->extension=NULL;          in->extension->extension=NULL;
153       }       }
154       Finley_IndexList_insertIndex(in->extension,index);       Dudley_IndexList_insertIndex(in->extension,index);
155    } else {    } else {
156       /* insert index into in->index*/       /* insert index into in->index*/
157       in->index[in->n]=index;       in->index[in->n]=index;
# Line 175  void Finley_IndexList_insertIndex(Finley Line 159  void Finley_IndexList_insertIndex(Finley
159    }    }
160  }  }
161
162  /* counts the number of row indices in the Finley_IndexList in */  /* counts the number of row indices in the Dudley_IndexList in */
163
164  int Finley_IndexList_count(Finley_IndexList* in) {  dim_t Dudley_IndexList_count(Dudley_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+Dudley_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 Dudley_IndexList in */
180
181  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, index_t index_offset) {
182    int 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        Dudley_IndexList_toArray(in->extension,&(array[ptr]), range_min, range_max, index_offset);
195    }    }
196  }  }
197
198  /* deallocates the Finley_IndexList in by recursive calls */  /* deallocates the Dudley_IndexList in by recursive calls */
199
200  void Finley_IndexList_free(Finley_IndexList* in) {  void Dudley_IndexList_free(Dudley_IndexList* in) {
201    if (in!=NULL) {    if (in!=NULL) {
202      Finley_IndexList_free(in->extension);      Dudley_IndexList_free(in->extension);
203      TMPMEMFREE(in);      TMPMEMFREE(in);
204    }    }
205  }  }
206
207  /*  /* creates a Paso_pattern from a range of indices */
208   * \$Log\$  Paso_Pattern* Dudley_IndexList_createPattern(dim_t n0, dim_t n,Dudley_IndexList* index_list,index_t range_min,index_t range_max,index_t index_offset)
209   * Revision 1.1  2004/10/26 06:53:57  jgs  {
210   * Initial revision     dim_t *ptr=NULL;
211   *     register dim_t s,i,itmp;
212   * Revision 1.1.2.2  2004/10/26 06:36:39  jgs     index_t *index=NULL;
213   * committing Lutz's changes to branch jgs     Paso_Pattern* out=NULL;
214   *
215   * Revision 1.2  2004/10/13 01:53:42  gross     ptr=MEMALLOC(n+1-n0,index_t);
216   * bug in CSC assembling fixed     if (! Dudley_checkPtr(ptr) ) {
217   *         /* get the number of connections per row */
218   * Revision 1.1  2004/07/02 04:21:13  gross         #pragma omp parallel for schedule(static) private(i)
219   * Finley C code has been included         for(i=n0;i<n;++i) {
220   *                ptr[i-n0]=Dudley_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 (! Dudley_checkPtr(index)) {
233                  #pragma omp parallel for schedule(static)
234                  for(i=n0;i<n;++i) {
235                      Dudley_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 (! Dudley_noError()) {
241            MEMFREE(ptr);
242            MEMFREE(index);
243            Paso_Pattern_free(out);
244      }
245      return out;
246    }

Legend:
 Removed from v.82 changed lines Added in v.3127