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

trunk/esys2/finley/src/finleyC/IndexList.c revision 82 by jgs, Tue Oct 26 06:53:54 2004 UTC temp/finley/src/IndexList.c revision 1387 by trankine, Fri Jan 11 07:45:26 2008 UTC
# Line 1  Line 1
/* \$Id\$ */
1
2  /**************************************************************/  /* \$Id\$ */
3
4  /* Finley: Converting an element list into a matrix shape     */  /*******************************************************
5     *
6     *           Copyright 2003-2007 by ACceSS MNRF
7     *       Copyright 2007 by University of Queensland
8     *
9     *                http://esscc.uq.edu.au
10     *        Primary Business: Queensland, Australia
13     *
14     *******************************************************/
15
16  /**************************************************************/  /**************************************************************/
17
18  /* Copyrights by ACcESS Australia 2003,2004 */  /* Finley: Converting an element list into a matrix shape     */
/* Author: gross@access.edu.au */
19
20  /**************************************************************/  /**************************************************************/
21
#include "Finley.h"
#include "ElementFile.h"
#include "System.h"
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 Finley_IndexList_insertElements(Finley_IndexList* index_list, Finley_ElementFile* elements,
32                                         int reduce_row_order, int packed_row_block_size, maybelong* row_Label,                                         bool_t reduce_row_order, index_t* row_map,
33                                         int reduce_col_order, int packed_col_block_size, maybelong* col_Label,                                         bool_t reduce_col_order, index_t* col_map) {
34                                         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 */
35    maybelong e,kr,ir,kc,ic,NN_row,NN_col,i;    index_t color, *id=NULL;
36    Finley_IndexList * tmp_list;    dim_t e,kr,kc,NN_row,NN_col,i,icol,irow, NN, *row_node=NULL,*col_node=NULL;
maybelong jr,jc,icol,irow,color;

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;      id=TMPMEMALLOC(NN, index_t);
40      for (i=0;i<NN;i++) id[i]=i;      if (! Finley_checkPtr(id) ) {
41      if (reduce_col_order) {         for (i=0;i<NN;i++) id[i]=i;
42         col_node=elements->ReferenceElement->Type->linearNodes;         if (reduce_col_order) {
43         NN_col=elements->LinearReferenceElement->Type->numNodes;            col_node=elements->ReferenceElement->Type->linearNodes;
44      } else {            NN_col=elements->LinearReferenceElement->Type->numNodes;
45         col_node=id;         } else {
46         NN_col=elements->ReferenceElement->Type->numNodes;            col_node=id;
47      }            NN_col=elements->ReferenceElement->Type->numNodes;
48      if (reduce_row_order) {         }
49         row_node=elements->ReferenceElement->Type->linearNodes;         if (reduce_row_order) {
50         NN_row=elements->LinearReferenceElement->Type->numNodes;            row_node=elements->ReferenceElement->Type->linearNodes;
51      } else {            NN_row=elements->LinearReferenceElement->Type->numNodes;
52         row_node=id;         } else {
53         NN_row=elements->ReferenceElement->Type->numNodes;            row_node=id;
54      }            NN_row=elements->ReferenceElement->Type->numNodes;
55           }
56      #pragma omp parallel private(color)         for (color=elements->minColor;color<=elements->maxColor;color++) {
57      {             #pragma omp for private(e,irow,kr,kc,icol) schedule(static)
58         switch(matType) {             for (e=0;e<elements->numElements;e++) {
59         case CSR:                 if (elements->Color[e]==color) {
60           if (symmetric) {                     for (kr=0;kr<NN_row;kr++) {
61              for (color=0;color<elements->numColors;color++) {                       irow=row_map[elements->Nodes[INDEX2(row_node[kr],e,NN)]];
#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;
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);
}
}
}
}
}
}
}
} /* if else */
break;
case CSC:
if (symmetric) {
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) {
62                       for (kc=0;kc<NN_col;kc++) {                       for (kc=0;kc<NN_col;kc++) {
63                          jc=col_Label[elements->Nodes[INDEX2(col_node[kc],e,NN)]];                            icol=col_map[elements->Nodes[INDEX2(col_node[kc],e,NN)]];
64                          for (ic=0;ic<packed_col_block_size;ic++) {                            Finley_IndexList_insertIndex(&(index_list[irow]),icol);
icol=ic+packed_col_block_size*jc;
tmp_list=&(index_list[icol]);
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;
if (irow<=icol) Finley_IndexList_insertIndex(tmp_list,irow);
}
}
}
65                       }                       }
66                    }                     }
67                 }                 }
68              }             }
69           } else {         }
70             for (color=0;color<elements->numColors;color++) {         TMPMEMFREE(id);
71                 #pragma omp for private (e,kr,jr,kc,jc,ir,irow,ic,icol,tmp_list) schedule(static)      }
72                 for (e=0;e<elements->numElements;e++) {    }
73                    if (elements->Color[e]==color) {    return;
74                       for (kc=0;kc<NN_col;kc++) {  }
75                          jc=col_Label[elements->Nodes[INDEX2(col_node[kc],e,NN)]];
76                          for (ic=0;ic<packed_col_block_size;ic++) {  void Finley_IndexList_insertElementsWithRowRange(Finley_IndexList* index_list, index_t firstRow, index_t lastRow,
77                             icol=ic+packed_col_block_size*jc;                                                   Finley_ElementFile* elements, index_t* row_map, index_t* col_map)
78                             tmp_list=&(index_list[icol]);  {
79                             for (kr=0;kr<NN_row;kr++) {    index_t color;
80                                jr=row_Label[elements->Nodes[INDEX2(row_node[kr],e,NN)]];    dim_t e,kr,kc,i,icol,irow, NN;
81                                for (ir=0;ir<packed_row_block_size;ir++) {    if (elements!=NULL) {
82                                   irow=ir+packed_row_block_size*jr;      NN=elements->numNodes;
83                                   Finley_IndexList_insertIndex(tmp_list,irow);      for (color=elements->minColor;color<=elements->maxColor;color++) {
84                                }             #pragma omp for private(e,irow,kr,kc,icol) schedule(static)
85                             }             for (e=0;e<elements->numElements;e++) {
86                          }                 if (elements->Color[e]==color) {
87                       }                     for (kr=0;kr<NN;kr++) {
88                         irow=row_map[elements->Nodes[INDEX2(kr,e,NN)]];
89                         if ((firstRow<=irow) && (irow < lastRow)) {
90                              irow-=firstRow;
91                              for (kc=0;kc<NN;kc++) {
92                                  icol=col_map[elements->Nodes[INDEX2(kc,e,NN)]];
93                                  Finley_IndexList_insertIndex(&(index_list[irow]),icol);
94                              }
95                          }
96                    }                    }
97                 }                 }
98              }             }
} /* if else */
} /* switch matType */
99      }      }
100    }    }
return;
101  }  }
102
103  /* 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 */
104
105  void Finley_IndexList_insertIndex(Finley_IndexList* in, maybelong index) {  void Finley_IndexList_insertIndex(Finley_IndexList* in, index_t index) {
106    int i;    dim_t i;
107    /* is index in in? */    /* is index in in? */
108    for (i=0;i<in->n;i++) {    for (i=0;i<in->n;i++) {
109      if (in->index[i]==index) return;      if (in->index[i]==index)  return;
110    }    }
111    /* index could not be found */    /* index could not be found */
112    if (in->n==INDEXLIST_LENGTH) {    if (in->n==INDEXLIST_LENGTH) {
113       /* if in->index is full check the extension */       /* if in->index is full check the extension */
114       if (in->extension==NULL) {       if (in->extension==NULL) {
115          in->extension=(Finley_IndexList*) TMPMEMALLOC(sizeof(Finley_IndexList));          in->extension=TMPMEMALLOC(1,Finley_IndexList);
116          if (Finley_checkPtr(in->extension)) return;          if (Finley_checkPtr(in->extension)) return;
117          in->extension->n=0;          in->extension->n=0;
118          in->extension->extension=NULL;          in->extension->extension=NULL;
# Line 177  void Finley_IndexList_insertIndex(Finley Line 127  void Finley_IndexList_insertIndex(Finley
127
128  /* counts the number of row indices in the Finley_IndexList in */  /* counts the number of row indices in the Finley_IndexList in */
129
130  int Finley_IndexList_count(Finley_IndexList* in) {  dim_t Finley_IndexList_count(Finley_IndexList* in, index_t range_min,index_t range_max) {
131      dim_t i;
132      dim_t out=0;
133      register index_t itmp;
134    if (in==NULL) {    if (in==NULL) {
135       return 0;       return 0;
136    } else {    } else {
137       return (in->n)+Finley_IndexList_count(in->extension);      for (i=0;i<in->n;i++) {
138              itmp=in->index[i];
139              if ((itmp>=range_min) && (range_max>itmp)) ++out;
140        }
141         return out+Finley_IndexList_count(in->extension, range_min,range_max);
142    }    }
143  }  }
144
145  /* count the number of row indices in the Finley_IndexList in */  /* count the number of row indices in the Finley_IndexList in */
146
147  void Finley_IndexList_toArray(Finley_IndexList* in, maybelong* array) {  void Finley_IndexList_toArray(Finley_IndexList* in, index_t* array, index_t range_min,index_t range_max, index_t index_offset) {
148    int i;    dim_t i, ptr;
149      register index_t itmp;
150    if (in!=NULL) {    if (in!=NULL) {
151      for (i=0;i<in->n;i++) array[i]=in->index[i]+INDEX_OFFSET;      ptr=0;
152      Finley_IndexList_toArray(in->extension,&(array[in->n]));      for (i=0;i<in->n;i++) {
153              itmp=in->index[i];
154              if ((itmp>=range_min) && (range_max>itmp)) {
155                 array[ptr]=itmp+index_offset;
156                 ptr++;
157              }
158
159        }
160        Finley_IndexList_toArray(in->extension,&(array[ptr]), range_min, range_max, index_offset);
161    }    }
162  }  }
163
# Line 204  void Finley_IndexList_free(Finley_IndexL Line 170  void Finley_IndexList_free(Finley_IndexL
170    }    }
171  }  }
172
173  /*  /* creates a Paso_pattern from a range of indices */
174   * \$Log\$  Paso_Pattern* Finley_IndexList_createPattern(dim_t n,Finley_IndexList* index_list,index_t range_min,index_t range_max,index_t index_offset)
175   * Revision 1.1  2004/10/26 06:53:57  jgs  {
176   * Initial revision     dim_t *ptr=NULL;
177   *     register dim_t s,i,itmp;
178   * Revision 1.1.2.2  2004/10/26 06:36:39  jgs     index_t *index=NULL;
179   * committing Lutz's changes to branch jgs     Paso_Pattern* out=NULL;
180   *
181   * Revision 1.2  2004/10/13 01:53:42  gross     ptr=MEMALLOC(n+1,index_t);
182   * bug in CSC assembling fixed     if (! Finley_checkPtr(ptr) ) {
183   *         /* get the number of connections per row */
184   * Revision 1.1  2004/07/02 04:21:13  gross         #pragma omp parallel for schedule(static) private(i)
185   * Finley C code has been included         for(i=0;i<n;++i) {
186   *                ptr[i]=Finley_IndexList_count(&index_list[i],range_min,range_max);
187   *         }
188   */         /* accumulate ptr */
189           s=0;
190           for(i=0;i<n;++i) {
191                   itmp=ptr[i];
192                   ptr[i]=s;
193                   s+=itmp;
194           }
195           ptr[n]=s;
196           /* fill index */
197           index=MEMALLOC(ptr[n],index_t);
198           if (! Finley_checkPtr(index)) {
199                  #pragma omp parallel for schedule(static)
200                  for(i=0;i<n;++i) {
201                      Finley_IndexList_toArray(&index_list[i],&index[ptr[i]],range_min,range_max,index_offset);
202                  }
203                  out=Paso_Pattern_alloc(PATTERN_FORMAT_DEFAULT,1,1,n,ptr,index);
204           }
205      }
206      if (! Finley_noError()) {
207            MEMFREE(ptr);
208            MEMFREE(index);
209            Paso_Pattern_free(out);
210      }
211      return out;
212    }

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