/[escript]/trunk/dudley/src/IndexList.cpp
ViewVC logotype

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

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

trunk/finley/src/finley/IndexList.c revision 201 by jgs, Wed Nov 23 04:10:21 2005 UTC branches/domexper/dudley/src/IndexList.c revision 3141 by jfenwick, Thu Sep 2 23:52:41 2010 UTC
# Line 1  Line 1 
 /*  
  ******************************************************************************  
  *                                                                            *  
  *       COPYRIGHT  ACcESS 2003,2004,2005 -  All Rights Reserved              *  
  *                                                                            *  
  * This software is the property of ACcESS. No part of this code              *  
  * may be copied in any form or by any means without the expressed written    *  
  * consent of ACcESS.  Copying, use or modification of this software          *  
  * by any unauthorised person is illegal unless that person has a software    *  
  * license agreement with ACcESS.                                             *  
  *                                                                            *  
  ******************************************************************************  
 */  
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
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  /*  Author: gross@access.edu.au */  /* Dudley: Converting an element list into a matrix shape     */
 /*  Version: $Id$ */  
18    
19  /**************************************************************/  /**************************************************************/
20    
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                                         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    {
34      /* index_list is an array of linked lists. Each entry is a row (DOF) and contains the indices to the non-zero columns */
35    index_t color;    index_t color;
36    dim_t e,kr,kc,NN_row,NN_col,i,icol,irow;    Dudley_ReferenceElement*refElement;
37      dim_t e, kr, kc, NN_row, NN_col, icol, irow, NN, *row_node=NULL, *col_node=NULL;
38    if (elements!=NULL) {    if (elements!=NULL)
39      dim_t NN=elements->ReferenceElement->Type->numNodes;    {
40      index_t id[NN],*row_node,*col_node;      NN=elements->numNodes;
41      for (i=0;i<NN;i++) id[i]=i;      refElement= Dudley_ReferenceElementSet_borrowReferenceElement(elements->referenceElementSet, FALSE);
42      if (reduce_col_order) {      if (reduce_col_order)
43         col_node=elements->ReferenceElement->Type->linearNodes;      {
44         NN_col=elements->LinearReferenceElement->Type->numNodes;      col_node=refElement->Type->linearNodes;
45      } else {      NN_col=(refElement->LinearBasisFunctions->Type->numShapes);
46         col_node=id;      }
47         NN_col=elements->ReferenceElement->Type->numNodes;      else
48        {
49        col_node=refElement->Type->subElementNodes;
50        NN_col=(refElement->BasisFunctions->Type->numShapes);
51      }      }
52      if (reduce_row_order) {      if (reduce_row_order)
53         row_node=elements->ReferenceElement->Type->linearNodes;      {
54         NN_row=elements->LinearReferenceElement->Type->numNodes;      row_node=refElement->Type->linearNodes;
55        NN_row=(refElement->LinearBasisFunctions->Type->numShapes);
56      } else {      } else {
57         row_node=id;      row_node=refElement->Type->subElementNodes;
58         NN_row=elements->ReferenceElement->Type->numNodes;      NN_row=(refElement->BasisFunctions->Type->numShapes) ;
59      }      }
60    
61        for (color=elements->minColor;color<=elements->maxColor;color++)
62        {
63        #pragma omp for private(e,irow,kr,kc,icol,isub) schedule(static)
64        for (e=0;e<elements->numElements;e++)
65        {
66            if (elements->Color[e]==color)
67            {
68                for (kr=0;kr<NN_row;kr++)
69                {
70                    irow=row_map[elements->Nodes[INDEX2(row_node[INDEX2(kr,0,NN_row)],e,NN)]];
71                    for (kc=0;kc<NN_col;kc++)
72                    {
73                    icol=col_map[elements->Nodes[INDEX2(col_node[INDEX2(kc,0,NN_col)],e,NN)]];
74                    Dudley_IndexList_insertIndex(&(index_list[irow]),icol);
75                    }
76                }
77            }
78            }
79        }
80      }
81      return;
82    }
83    
84    
85    void Dudley_IndexList_insertElementsWithRowRange(Dudley_IndexList* index_list, index_t firstRow, index_t lastRow,
86                                                     Dudley_ElementFile* elements, index_t* row_map, index_t* col_map)
87    {
88    /* this does not resolve macro elements */
89        index_t color;
90      dim_t e,kr,kc,icol,irow, NN;
91      if (elements!=NULL) {
92        NN=elements->numNodes;
93      for (color=elements->minColor;color<=elements->maxColor;color++) {      for (color=elements->minColor;color<=elements->maxColor;color++) {
94          #pragma omp for private(e,irow,kr,kc,icol) schedule(static)             #pragma omp for private(e,irow,kr,kc,icol) schedule(static)
95          for (e=0;e<elements->numElements;e++) {             for (e=0;e<elements->numElements;e++) {
96              if (elements->Color[e]==color) {                 if (elements->Color[e]==color) {
97                  for (kr=0;kr<NN_row;kr++) {                     for (kr=0;kr<NN;kr++) {
98                    irow=row_Label[elements->Nodes[INDEX2(row_node[kr],e,NN)]];                       irow=row_map[elements->Nodes[INDEX2(kr,e,NN)]];
99                    for (kc=0;kc<NN_col;kc++) {                       if ((firstRow<=irow) && (irow < lastRow)) {
100                         icol=col_Label[elements->Nodes[INDEX2(col_node[kc],e,NN)]];                            irow-=firstRow;
101                         Finley_IndexList_insertIndex(&(index_list[irow]),icol);                            for (kc=0;kc<NN;kc++) {
102                                  icol=col_map[elements->Nodes[INDEX2(kc,e,NN)]];
103                                  Dudley_IndexList_insertIndex(&(index_list[irow]),icol);
104                              }
105                          }
106                    }                    }
107                  }                 }
108              }             }
109          }      }
110        }    }
111    }
112    void Dudley_IndexList_insertElementsWithRowRangeNoMainDiagonal(Dudley_IndexList* index_list, index_t firstRow, index_t lastRow,
113                                                                  Dudley_ElementFile* elements, index_t* row_map, index_t* col_map)
114    {
115      /* this does not resolve macro elements */
116      index_t color;
117      dim_t e,kr,kc,icol,irow, NN,irow_loc;
118      if (elements!=NULL) {
119        NN=elements->numNodes;
120        for (color=elements->minColor;color<=elements->maxColor;color++) {
121               #pragma omp for private(e,irow,kr,kc,icol,irow_loc) schedule(static)
122               for (e=0;e<elements->numElements;e++) {
123                   if (elements->Color[e]==color) {
124                       for (kr=0;kr<NN;kr++) {
125                         irow=row_map[elements->Nodes[INDEX2(kr,e,NN)]];
126                         if ((firstRow<=irow) && (irow < lastRow)) {
127                              irow_loc=irow-firstRow;
128                              for (kc=0;kc<NN;kc++) {
129                                  icol=col_map[elements->Nodes[INDEX2(kc,e,NN)]];
130                                  if (icol != irow) Dudley_IndexList_insertIndex(&(index_list[irow_loc]),icol);
131                              }
132                          }
133                      }
134                   }
135               }
136        }
137    }    }
   return;  
138  }  }
139    
140  /* 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 */
141    
142  void Finley_IndexList_insertIndex(Finley_IndexList* in, index_t index) {  void Dudley_IndexList_insertIndex(Dudley_IndexList* in, index_t index) {
143    dim_t i;    dim_t i;
144    /* is index in in? */    /* is index in in? */
145    for (i=0;i<in->n;i++) {    for (i=0;i<in->n;i++) {
# Line 84  void Finley_IndexList_insertIndex(Finley Line 149  void Finley_IndexList_insertIndex(Finley
149    if (in->n==INDEXLIST_LENGTH) {    if (in->n==INDEXLIST_LENGTH) {
150       /* if in->index is full check the extension */       /* if in->index is full check the extension */
151       if (in->extension==NULL) {       if (in->extension==NULL) {
152          in->extension=TMPMEMALLOC(1,Finley_IndexList);          in->extension=TMPMEMALLOC(1,Dudley_IndexList);
153          if (Finley_checkPtr(in->extension)) return;          if (Dudley_checkPtr(in->extension)) return;
154          in->extension->n=0;          in->extension->n=0;
155          in->extension->extension=NULL;          in->extension->extension=NULL;
156       }       }
157       Finley_IndexList_insertIndex(in->extension,index);       Dudley_IndexList_insertIndex(in->extension,index);
158    } else {    } else {
159       /* insert index into in->index*/       /* insert index into in->index*/
160       in->index[in->n]=index;       in->index[in->n]=index;
# Line 97  void Finley_IndexList_insertIndex(Finley Line 162  void Finley_IndexList_insertIndex(Finley
162    }    }
163  }  }
164    
165  /* counts the number of row indices in the Finley_IndexList in */  /* counts the number of row indices in the Dudley_IndexList in */
166    
167  dim_t Finley_IndexList_count(Finley_IndexList* in) {  dim_t Dudley_IndexList_count(Dudley_IndexList* in, index_t range_min,index_t range_max) {
168      dim_t i;
169      dim_t out=0;
170      register index_t itmp;
171    if (in==NULL) {    if (in==NULL) {
172       return 0;       return 0;
173    } else {    } else {
174       return (in->n)+Finley_IndexList_count(in->extension);      for (i=0;i<in->n;i++) {
175              itmp=in->index[i];
176              if ((itmp>=range_min) && (range_max>itmp)) ++out;
177        }
178         return out+Dudley_IndexList_count(in->extension, range_min,range_max);
179    }    }
180  }  }
181    
182  /* count the number of row indices in the Finley_IndexList in */  /* count the number of row indices in the Dudley_IndexList in */
183    
184  void Finley_IndexList_toArray(Finley_IndexList* in, index_t* array) {  void Dudley_IndexList_toArray(Dudley_IndexList* in, index_t* array, index_t range_min,index_t range_max, index_t index_offset) {
185    dim_t i;    dim_t i, ptr;
186      register index_t itmp;
187    if (in!=NULL) {    if (in!=NULL) {
188      for (i=0;i<in->n;i++) array[i]=in->index[i]+INDEX_OFFSET;      ptr=0;
189      Finley_IndexList_toArray(in->extension,&(array[in->n]));      for (i=0;i<in->n;i++) {
190              itmp=in->index[i];
191              if ((itmp>=range_min) && (range_max>itmp)) {
192                 array[ptr]=itmp+index_offset;
193                 ptr++;
194              }
195    
196        }
197        Dudley_IndexList_toArray(in->extension,&(array[ptr]), range_min, range_max, index_offset);
198    }    }
199  }  }
200    
201  /* deallocates the Finley_IndexList in by recursive calls */  /* deallocates the Dudley_IndexList in by recursive calls */
202    
203  void Finley_IndexList_free(Finley_IndexList* in) {  void Dudley_IndexList_free(Dudley_IndexList* in) {
204    if (in!=NULL) {    if (in!=NULL) {
205      Finley_IndexList_free(in->extension);      Dudley_IndexList_free(in->extension);
206      TMPMEMFREE(in);      TMPMEMFREE(in);
207    }    }
208  }  }
209    
210  /*  /* creates a Paso_pattern from a range of indices */
211   * $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)
212   * Revision 1.6  2005/09/15 03:44:22  jgs  {
213   * Merge of development branch dev-02 back to main trunk on 2005-09-15     dim_t *ptr=NULL;
214   *     register dim_t s,i,itmp;
215   * Revision 1.5.2.1  2005/09/07 06:26:18  gross     index_t *index=NULL;
216   * the solver from finley are put into the standalone package paso now     Paso_Pattern* out=NULL;
217   *  
218   * Revision 1.5  2005/07/08 04:07:51  jgs     ptr=MEMALLOC(n+1-n0,index_t);
219   * Merge of development branch back to main trunk on 2005-07-08     if (! Dudley_checkPtr(ptr) ) {
220   *         /* get the number of connections per row */
221   * Revision 1.4  2004/12/15 07:08:32  jgs         #pragma omp parallel for schedule(static) private(i)
222   * *** empty log message ***         for(i=n0;i<n;++i) {
223   * Revision 1.1.1.1.2.3  2005/06/29 02:34:50  gross                ptr[i-n0]=Dudley_IndexList_count(&index_list[i],range_min,range_max);
224   * some changes towards 64 integers in finley         }
225   *         /* accumulate ptr */
226   * Revision 1.1.1.1.2.2  2004/11/24 01:37:13  gross         s=0;
227   * some changes dealing with the integer overflow in memory allocation. Finley solves 4M unknowns now         for(i=n0;i<n;++i) {
228   *                 itmp=ptr[i-n0];
229   *                 ptr[i-n0]=s;
230   *                 s+=itmp;
231   */         }
232           ptr[n-n0]=s;
233           /* fill index */
234           index=MEMALLOC(ptr[n-n0],index_t);
235           if (! Dudley_checkPtr(index)) {
236                  #pragma omp parallel for schedule(static)
237                  for(i=n0;i<n;++i) {
238                      Dudley_IndexList_toArray(&index_list[i],&index[ptr[i-n0]],range_min,range_max,index_offset);
239                  }
240                  out=Paso_Pattern_alloc(PATTERN_FORMAT_DEFAULT,n-n0,range_max+index_offset,ptr,index);
241           }
242      }
243      if (! Dudley_noError()) {
244            MEMFREE(ptr);
245            MEMFREE(index);
246            Paso_Pattern_free(out);
247      }
248      return out;
249    }

Legend:
Removed from v.201  
changed lines
  Added in v.3141

  ViewVC Help
Powered by ViewVC 1.1.26