# Diff of /trunk/dudley/src/IndexList.cpp

revision 971 by ksteube, Wed Feb 14 04:40:49 2007 UTC revision 1722 by gross, Fri Aug 22 04:20:30 2008 UTC
# Line 1  Line 1
/*
************************************************************
*          Copyright 2006 by ACcESS MNRF                   *
*                                                          *
*              http://www.access.edu.au                    *
*       Primary Business: Queensland, Australia            *
*                                                          *
************************************************************
*/
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  /*  Author: gross@access.edu.au */  /* Finley: Converting an element list into a matrix shape     */
/*  Version: \$Id\$ */
19
20  /**************************************************************/  /**************************************************************/
21
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                                         bool_t reduce_row_order, index_t* row_Label,                                         bool_t reduce_row_order, index_t* row_map,
33                                         bool_t reduce_col_order, index_t* col_Label) {                                         bool_t reduce_col_order, index_t* col_map) {
34    index_t color;    /* index_list is an array of linked lists. Each entry is a row (DOF) and contains the indices to the non-zero columns */
35    dim_t e,kr,kc,NN_row,NN_col,i,icol,irow;    index_t color, *id=NULL;
36      dim_t e,kr,kc,NN_row,NN_col,i,icol,irow, NN, *row_node=NULL,*col_node=NULL;
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;      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           }
49           if (reduce_row_order) {
50              row_node=elements->ReferenceElement->Type->linearNodes;
51              NN_row=elements->LinearReferenceElement->Type->numNodes;
52           } else {
53              row_node=id;
54              NN_row=elements->ReferenceElement->Type->numNodes;
55           }
56           for (color=elements->minColor;color<=elements->maxColor;color++) {
57               #pragma omp for private(e,irow,kr,kc,icol) schedule(static)
58               for (e=0;e<elements->numElements;e++) {
59                   if (elements->Color[e]==color) {
60                       for (kr=0;kr<NN_row;kr++) {
61                         irow=row_map[elements->Nodes[INDEX2(row_node[kr],e,NN)]];
62                         for (kc=0;kc<NN_col;kc++) {
63                              icol=col_map[elements->Nodes[INDEX2(col_node[kc],e,NN)]];
64                              Finley_IndexList_insertIndex(&(index_list[irow]),icol);
65                         }
66                       }
67                   }
68               }
69           }
70           TMPMEMFREE(id);
71      }      }
72      if (reduce_row_order) {    }
73         row_node=elements->ReferenceElement->Type->linearNodes;    return;
74         NN_row=elements->LinearReferenceElement->Type->numNodes;  }
75      } else {
76         row_node=id;  void Finley_IndexList_insertElementsWithRowRange(Finley_IndexList* index_list, index_t firstRow, index_t lastRow,
77         NN_row=elements->ReferenceElement->Type->numNodes;                                                   Finley_ElementFile* elements, index_t* row_map, index_t* col_map)
78    {
79      index_t color;
80      dim_t e,kr,kc,icol,irow, NN;
81      if (elements!=NULL) {
82        NN=elements->numNodes;
83        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               }
99      }      }
100      }
101    }
102    void Finley_IndexList_insertElementsWithRowRangeNoMainDiagonal(Finley_IndexList* index_list, index_t firstRow, index_t lastRow,
103                                                                  Finley_ElementFile* elements, index_t* row_map, index_t* col_map)
104    {
105      index_t color;
106      dim_t e,kr,kc,icol,irow, NN,irow_loc;
107      if (elements!=NULL) {
108        NN=elements->numNodes;
109      for (color=elements->minColor;color<=elements->maxColor;color++) {      for (color=elements->minColor;color<=elements->maxColor;color++) {
110          #pragma omp for private(e,irow,kr,kc,icol) schedule(static)             #pragma omp for private(e,irow,kr,kc,icol,irow_loc) schedule(static)
111          for (e=0;e<elements->numElements;e++) {             for (e=0;e<elements->numElements;e++) {
112              if (elements->Color[e]==color) {                 if (elements->Color[e]==color) {
113                  for (kr=0;kr<NN_row;kr++) {                     for (kr=0;kr<NN;kr++) {
114                    irow=row_Label[elements->Nodes[INDEX2(row_node[kr],e,NN)]];                       irow=row_map[elements->Nodes[INDEX2(kr,e,NN)]];
115                    for (kc=0;kc<NN_col;kc++) {                       if ((firstRow<=irow) && (irow < lastRow)) {
116                         icol=col_Label[elements->Nodes[INDEX2(col_node[kc],e,NN)]];                            irow_loc=irow-firstRow;
117                         Finley_IndexList_insertIndex(&(index_list[irow]),icol);                            for (kc=0;kc<NN;kc++) {
118                                  icol=col_map[elements->Nodes[INDEX2(kc,e,NN)]];
119                                  if (icol != irow) Finley_IndexList_insertIndex(&(index_list[irow_loc]),icol);
120                              }
121                          }
122                    }                    }
123                  }                 }
124              }             }
125          }      }
}
126    }    }
return;
127  }  }
128
129  /* 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 97  void Finley_IndexList_insertIndex(Finley Line 153  void Finley_IndexList_insertIndex(Finley
153
154  /* counts the number of row indices in the Finley_IndexList in */  /* counts the number of row indices in the Finley_IndexList in */
155
156  dim_t Finley_IndexList_count(Finley_IndexList* in) {  dim_t Finley_IndexList_count(Finley_IndexList* in, index_t range_min,index_t range_max) {
157      dim_t i;
158      dim_t out=0;
159      register index_t itmp;
160    if (in==NULL) {    if (in==NULL) {
161       return 0;       return 0;
162    } else {    } else {
163       return (in->n)+Finley_IndexList_count(in->extension);      for (i=0;i<in->n;i++) {
164              itmp=in->index[i];
165              if ((itmp>=range_min) && (range_max>itmp)) ++out;
166        }
167         return out+Finley_IndexList_count(in->extension, range_min,range_max);
168    }    }
169  }  }
170
171  /* count the number of row indices in the Finley_IndexList in */  /* count the number of row indices in the Finley_IndexList in */
172
173  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) {
174    dim_t i;    dim_t i, ptr;
175      register index_t itmp;
176    if (in!=NULL) {    if (in!=NULL) {
177      for (i=0;i<in->n;i++) array[i]=in->index[i];      ptr=0;
178      Finley_IndexList_toArray(in->extension,&(array[in->n]));      for (i=0;i<in->n;i++) {
179              itmp=in->index[i];
180              if ((itmp>=range_min) && (range_max>itmp)) {
181                 array[ptr]=itmp+index_offset;
182                 ptr++;
183              }
184
185        }
186        Finley_IndexList_toArray(in->extension,&(array[ptr]), range_min, range_max, index_offset);
187    }    }
188  }  }
189
# Line 124  void Finley_IndexList_free(Finley_IndexL Line 196  void Finley_IndexList_free(Finley_IndexL
196    }    }
197  }  }
198
199  /*  /* creates a Paso_pattern from a range of indices */
200   * \$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)
201   * Revision 1.6  2005/09/15 03:44:22  jgs  {
202   * Merge of development branch dev-02 back to main trunk on 2005-09-15     dim_t *ptr=NULL;
203   *     register dim_t s,i,itmp;
204   * Revision 1.5.2.1  2005/09/07 06:26:18  gross     index_t *index=NULL;
205   * the solver from finley are put into the standalone package paso now     Paso_Pattern* out=NULL;
206   *
207   * Revision 1.5  2005/07/08 04:07:51  jgs     ptr=MEMALLOC(n+1-n0,index_t);
208   * Merge of development branch back to main trunk on 2005-07-08     if (! Finley_checkPtr(ptr) ) {
209   *         /* get the number of connections per row */
210   * Revision 1.4  2004/12/15 07:08:32  jgs         #pragma omp parallel for schedule(static) private(i)
211   * *** empty log message ***         for(i=n0;i<n;++i) {
212   * Revision 1.1.1.1.2.3  2005/06/29 02:34:50  gross                ptr[i-n0]=Finley_IndexList_count(&index_list[i],range_min,range_max);
213   * some changes towards 64 integers in finley         }
214   *         /* accumulate ptr */
215   * Revision 1.1.1.1.2.2  2004/11/24 01:37:13  gross         s=0;
216   * some changes dealing with the integer overflow in memory allocation. Finley solves 4M unknowns now         for(i=n0;i<n;++i) {
217   *                 itmp=ptr[i-n0];
218   *                 ptr[i-n0]=s;
219   *                 s+=itmp;
220   */         }
221           ptr[n-n0]=s;
222           /* fill index */
223           index=MEMALLOC(ptr[n-n0],index_t);
224           if (! Finley_checkPtr(index)) {
225                  #pragma omp parallel for schedule(static)
226                  for(i=n0;i<n;++i) {
227                      Finley_IndexList_toArray(&index_list[i],&index[ptr[i-n0]],range_min,range_max,index_offset);
228                  }
229                  out=Paso_Pattern_alloc(PATTERN_FORMAT_DEFAULT,1,1,n-n0,ptr,index);
230           }
231      }
232      if (! Finley_noError()) {
233            MEMFREE(ptr);
234            MEMFREE(index);
235            Paso_Pattern_free(out);
236      }
237      return out;
238    }

Legend:
 Removed from v.971 changed lines Added in v.1722