/[escript]/trunk/finley/src/IndexList.c
ViewVC logotype

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

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 2551 by gross, Thu Jul 23 09:19:15 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_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 */
34    dim_t e,kr,kc,NN_row,NN_col,i,icol,irow;    index_t color, *id=NULL;
35      dim_t e,kr,kc,NN_row,NN_col,i,icol,irow, NN, *row_node=NULL,*col_node=NULL;
36    if (elements!=NULL) {    if (elements!=NULL) {
37      dim_t NN=elements->ReferenceElement->Type->numNodes;      NN=elements->numNodes;
38      index_t id[NN],*row_node,*col_node;      id=TMPMEMALLOC(NN, index_t);
39      for (i=0;i<NN;i++) id[i]=i;      if (! Finley_checkPtr(id) ) {
40      if (reduce_col_order) {         for (i=0;i<NN;i++) id[i]=i;
41         col_node=elements->ReferenceElement->Type->linearNodes;         if (reduce_col_order) {
42         NN_col=elements->LinearReferenceElement->Type->numNodes;            col_node=elements->ReferenceElement->Type->linearNodes;
43      } else {            NN_col=elements->LinearReferenceElement->Type->numNodes;
44         col_node=id;         } else {
45         NN_col=elements->ReferenceElement->Type->numNodes;            col_node=id;
46              NN_col=elements->ReferenceElement->Type->numNodes;
47           }
48           if (reduce_row_order) {
49              row_node=elements->ReferenceElement->Type->linearNodes;
50              NN_row=elements->LinearReferenceElement->Type->numNodes;
51           } else {
52              row_node=id;
53              NN_row=elements->ReferenceElement->Type->numNodes;
54           }
55           for (color=elements->minColor;color<=elements->maxColor;color++) {
56               #pragma omp for private(e,irow,kr,kc,icol) schedule(static)
57               for (e=0;e<elements->numElements;e++) {
58                   if (elements->Color[e]==color) {
59                       for (kr=0;kr<NN_row;kr++) {
60                         irow=row_map[elements->Nodes[INDEX2(row_node[kr],e,NN)]];
61                         for (kc=0;kc<NN_col;kc++) {
62                              icol=col_map[elements->Nodes[INDEX2(col_node[kc],e,NN)]];
63                              Finley_IndexList_insertIndex(&(index_list[irow]),icol);
64                         }
65                       }
66                   }
67               }
68           }
69           TMPMEMFREE(id);
70      }      }
71      if (reduce_row_order) {    }
72         row_node=elements->ReferenceElement->Type->linearNodes;    return;
73         NN_row=elements->LinearReferenceElement->Type->numNodes;  }
74      } else {  
75         row_node=id;  void Finley_IndexList_insertElementsWithRowRange(Finley_IndexList* index_list, index_t firstRow, index_t lastRow,
76         NN_row=elements->ReferenceElement->Type->numNodes;                                                   Finley_ElementFile* elements, index_t* row_map, index_t* col_map)
77    {
78      index_t color;
79      dim_t e,kr,kc,icol,irow, NN;
80      if (elements!=NULL) {
81        NN=elements->numNodes;
82        for (color=elements->minColor;color<=elements->maxColor;color++) {
83               #pragma omp for private(e,irow,kr,kc,icol) schedule(static)
84               for (e=0;e<elements->numElements;e++) {
85                   if (elements->Color[e]==color) {
86                       for (kr=0;kr<NN;kr++) {
87                         irow=row_map[elements->Nodes[INDEX2(kr,e,NN)]];
88                         if ((firstRow<=irow) && (irow < lastRow)) {
89                              irow-=firstRow;
90                              for (kc=0;kc<NN;kc++) {
91                                  icol=col_map[elements->Nodes[INDEX2(kc,e,NN)]];
92                                  Finley_IndexList_insertIndex(&(index_list[irow]),icol);
93                              }
94                          }
95                      }
96                   }
97               }
98      }      }
99      }
100    }
101    void Finley_IndexList_insertElementsWithRowRangeNoMainDiagonal(Finley_IndexList* index_list, index_t firstRow, index_t lastRow,
102                                                                  Finley_ElementFile* elements, index_t* row_map, index_t* col_map)
103    {
104      index_t color;
105      dim_t e,kr,kc,icol,irow, NN,irow_loc;
106      if (elements!=NULL) {
107        NN=elements->numNodes;
108      for (color=elements->minColor;color<=elements->maxColor;color++) {      for (color=elements->minColor;color<=elements->maxColor;color++) {
109          #pragma omp for private(e,irow,kr,kc,icol) schedule(static)             #pragma omp for private(e,irow,kr,kc,icol,irow_loc) schedule(static)
110          for (e=0;e<elements->numElements;e++) {             for (e=0;e<elements->numElements;e++) {
111              if (elements->Color[e]==color) {                 if (elements->Color[e]==color) {
112                  for (kr=0;kr<NN_row;kr++) {                     for (kr=0;kr<NN;kr++) {
113                    irow=row_Label[elements->Nodes[INDEX2(row_node[kr],e,NN)]];                       irow=row_map[elements->Nodes[INDEX2(kr,e,NN)]];
114                    for (kc=0;kc<NN_col;kc++) {                       if ((firstRow<=irow) && (irow < lastRow)) {
115                         icol=col_Label[elements->Nodes[INDEX2(col_node[kc],e,NN)]];                            irow_loc=irow-firstRow;
116                         Finley_IndexList_insertIndex(&(index_list[irow]),icol);                            for (kc=0;kc<NN;kc++) {
117                                  icol=col_map[elements->Nodes[INDEX2(kc,e,NN)]];
118                                  if (icol != irow) Finley_IndexList_insertIndex(&(index_list[irow_loc]),icol);
119                              }
120                          }
121                    }                    }
122                  }                 }
123              }             }
124          }      }
       }  
125    }    }
   return;  
126  }  }
127    
128  /* 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 152  void Finley_IndexList_insertIndex(Finley
152    
153  /* counts the number of row indices in the Finley_IndexList in */  /* counts the number of row indices in the Finley_IndexList in */
154    
155  dim_t Finley_IndexList_count(Finley_IndexList* in) {  dim_t Finley_IndexList_count(Finley_IndexList* in, index_t range_min,index_t range_max) {
156      dim_t i;
157      dim_t out=0;
158      register index_t itmp;
159    if (in==NULL) {    if (in==NULL) {
160       return 0;       return 0;
161    } else {    } else {
162       return (in->n)+Finley_IndexList_count(in->extension);      for (i=0;i<in->n;i++) {
163              itmp=in->index[i];
164              if ((itmp>=range_min) && (range_max>itmp)) ++out;
165        }
166         return out+Finley_IndexList_count(in->extension, range_min,range_max);
167    }    }
168  }  }
169    
170  /* count the number of row indices in the Finley_IndexList in */  /* count the number of row indices in the Finley_IndexList in */
171    
172  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) {
173    dim_t i;    dim_t i, ptr;
174      register index_t itmp;
175    if (in!=NULL) {    if (in!=NULL) {
176      for (i=0;i<in->n;i++) array[i]=in->index[i]+INDEX_OFFSET;      ptr=0;
177      Finley_IndexList_toArray(in->extension,&(array[in->n]));      for (i=0;i<in->n;i++) {
178              itmp=in->index[i];
179              if ((itmp>=range_min) && (range_max>itmp)) {
180                 array[ptr]=itmp+index_offset;
181                 ptr++;
182              }
183    
184        }
185        Finley_IndexList_toArray(in->extension,&(array[ptr]), range_min, range_max, index_offset);
186    }    }
187  }  }
188    
# Line 117  void Finley_IndexList_free(Finley_IndexL Line 195  void Finley_IndexList_free(Finley_IndexL
195    }    }
196  }  }
197    
198  /*  /* creates a Paso_pattern from a range of indices */
199   * $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)
200   * Revision 1.5  2005/07/08 04:07:51  jgs  {
201   * Merge of development branch back to main trunk on 2005-07-08     dim_t *ptr=NULL;
202   *     register dim_t s,i,itmp;
203   * Revision 1.4  2004/12/15 07:08:32  jgs     index_t *index=NULL;
204   * *** empty log message ***     Paso_Pattern* out=NULL;
205   * Revision 1.1.1.1.2.3  2005/06/29 02:34:50  gross  
206   * some changes towards 64 integers in finley     ptr=MEMALLOC(n+1-n0,index_t);
207   *     if (! Finley_checkPtr(ptr) ) {
208   * Revision 1.1.1.1.2.2  2004/11/24 01:37:13  gross         /* get the number of connections per row */
209   * some changes dealing with the integer overflow in memory allocation. Finley solves 4M unknowns now         #pragma omp parallel for schedule(static) private(i)
210   *         for(i=n0;i<n;++i) {
211   *                ptr[i-n0]=Finley_IndexList_count(&index_list[i],range_min,range_max);
212   *         }
213   */         /* accumulate ptr */
214           s=0;
215           for(i=n0;i<n;++i) {
216                   itmp=ptr[i-n0];
217                   ptr[i-n0]=s;
218                   s+=itmp;
219           }
220           ptr[n-n0]=s;
221           /* fill index */
222           index=MEMALLOC(ptr[n-n0],index_t);
223           if (! Finley_checkPtr(index)) {
224                  #pragma omp parallel for schedule(static)
225                  for(i=n0;i<n;++i) {
226                      Finley_IndexList_toArray(&index_list[i],&index[ptr[i-n0]],range_min,range_max,index_offset);
227                  }
228                  out=Paso_Pattern_alloc(PATTERN_FORMAT_DEFAULT,n-n0,range_max+index_offset,ptr,index);
229           }
230      }
231      if (! Finley_noError()) {
232            MEMFREE(ptr);
233            MEMFREE(index);
234            Paso_Pattern_free(out);
235      }
236      return out;
237    }

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

  ViewVC Help
Powered by ViewVC 1.1.26