/[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/IndexList.c revision 971 by ksteube, Wed Feb 14 04:40:49 2007 UTC branches/domexper/dudley/src/IndexList.c revision 3221 by jfenwick, Wed Sep 29 01:00:21 2010 UTC
# Line 1  Line 1 
 /*  
  ************************************************************  
  *          Copyright 2006 by ACcESS MNRF                   *  
  *                                                          *  
  *              http://www.access.edu.au                    *  
  *       Primary Business: Queensland, Australia            *  
  *  Licensed under the Open Software License version 3.0    *  
  *     http://www.opensource.org/licenses/osl-3.0.php       *  
  *                                                          *  
  ************************************************************  
 */  
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;    dim_t e, kr, kc, NN_row, NN_col, icol, irow, NN;
37      if (elements!=NULL)
38      {
39        NN=elements->numNodes;
40        NN_col=(elements->numShapes);
41        NN_row=(elements->numShapes) ;
42    
43        for (color=elements->minColor;color<=elements->maxColor;color++)
44        {
45        #pragma omp for private(e,irow,kr,kc,icol) schedule(static)
46        for (e=0;e<elements->numElements;e++)
47        {
48            if (elements->Color[e]==color)
49            {
50                for (kr=0;kr<NN_row;kr++)
51                {
52                    irow=row_map[elements->Nodes[INDEX2(kr,e,NN)]];
53                    for (kc=0;kc<NN_col;kc++)
54                    {
55                    icol=col_map[elements->Nodes[INDEX2(kc,e,NN)]];
56                    Dudley_IndexList_insertIndex(&(index_list[irow]),icol);
57                    }
58                }
59            }
60        }
61         }
62      }
63      return;
64    }
65    
66    
67    void Dudley_IndexList_insertElementsWithRowRange(Dudley_IndexList* index_list, index_t firstRow, index_t lastRow,
68                                                     Dudley_ElementFile* elements, index_t* row_map, index_t* col_map)
69    {
70    /* this does not resolve macro elements */
71        index_t color;
72      dim_t e,kr,kc,icol,irow, NN;
73    if (elements!=NULL) {    if (elements!=NULL) {
74      dim_t NN=elements->ReferenceElement->Type->numNodes;      NN=elements->numNodes;
75      index_t id[NN],*row_node,*col_node;      for (color=elements->minColor;color<=elements->maxColor;color++) {
76      for (i=0;i<NN;i++) id[i]=i;             #pragma omp for private(e,irow,kr,kc,icol) schedule(static)
77      if (reduce_col_order) {             for (e=0;e<elements->numElements;e++) {
78         col_node=elements->ReferenceElement->Type->linearNodes;                 if (elements->Color[e]==color) {
79         NN_col=elements->LinearReferenceElement->Type->numNodes;                     for (kr=0;kr<NN;kr++) {
80      } else {                       irow=row_map[elements->Nodes[INDEX2(kr,e,NN)]];
81         col_node=id;                       if ((firstRow<=irow) && (irow < lastRow)) {
82         NN_col=elements->ReferenceElement->Type->numNodes;                            irow-=firstRow;
83      }                            for (kc=0;kc<NN;kc++) {
84      if (reduce_row_order) {                                icol=col_map[elements->Nodes[INDEX2(kc,e,NN)]];
85         row_node=elements->ReferenceElement->Type->linearNodes;                                Dudley_IndexList_insertIndex(&(index_list[irow]),icol);
86         NN_row=elements->LinearReferenceElement->Type->numNodes;                            }
87      } else {                        }
88         row_node=id;                    }
89         NN_row=elements->ReferenceElement->Type->numNodes;                 }
90               }
91      }      }
92      }
93    }
94    void Dudley_IndexList_insertElementsWithRowRangeNoMainDiagonal(Dudley_IndexList* index_list, index_t firstRow, index_t lastRow,
95                                                                  Dudley_ElementFile* elements, index_t* row_map, index_t* col_map)
96    {
97      /* this does not resolve macro elements */
98      index_t color;
99      dim_t e,kr,kc,icol,irow, NN,irow_loc;
100      if (elements!=NULL) {
101        NN=elements->numNodes;
102      for (color=elements->minColor;color<=elements->maxColor;color++) {      for (color=elements->minColor;color<=elements->maxColor;color++) {
103          #pragma omp for private(e,irow,kr,kc,icol) schedule(static)             #pragma omp for private(e,irow,kr,kc,icol,irow_loc) schedule(static)
104          for (e=0;e<elements->numElements;e++) {             for (e=0;e<elements->numElements;e++) {
105              if (elements->Color[e]==color) {                 if (elements->Color[e]==color) {
106                  for (kr=0;kr<NN_row;kr++) {                     for (kr=0;kr<NN;kr++) {
107                    irow=row_Label[elements->Nodes[INDEX2(row_node[kr],e,NN)]];                       irow=row_map[elements->Nodes[INDEX2(kr,e,NN)]];
108                    for (kc=0;kc<NN_col;kc++) {                       if ((firstRow<=irow) && (irow < lastRow)) {
109                         icol=col_Label[elements->Nodes[INDEX2(col_node[kc],e,NN)]];                            irow_loc=irow-firstRow;
110                         Finley_IndexList_insertIndex(&(index_list[irow]),icol);                            for (kc=0;kc<NN;kc++) {
111                                  icol=col_map[elements->Nodes[INDEX2(kc,e,NN)]];
112                                  if (icol != irow) Dudley_IndexList_insertIndex(&(index_list[irow_loc]),icol);
113                              }
114                          }
115                    }                    }
116                  }                 }
117              }             }
118          }      }
       }  
119    }    }
   return;  
120  }  }
121    
122  /* 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 */
123    
124  void Finley_IndexList_insertIndex(Finley_IndexList* in, index_t index) {  void Dudley_IndexList_insertIndex(Dudley_IndexList* in, index_t index) {
125    dim_t i;    dim_t i;
126    /* is index in in? */    /* is index in in? */
127    for (i=0;i<in->n;i++) {    for (i=0;i<in->n;i++) {
# Line 82  void Finley_IndexList_insertIndex(Finley Line 131  void Finley_IndexList_insertIndex(Finley
131    if (in->n==INDEXLIST_LENGTH) {    if (in->n==INDEXLIST_LENGTH) {
132       /* if in->index is full check the extension */       /* if in->index is full check the extension */
133       if (in->extension==NULL) {       if (in->extension==NULL) {
134          in->extension=TMPMEMALLOC(1,Finley_IndexList);          in->extension=TMPMEMALLOC(1,Dudley_IndexList);
135          if (Finley_checkPtr(in->extension)) return;          if (Dudley_checkPtr(in->extension)) return;
136          in->extension->n=0;          in->extension->n=0;
137          in->extension->extension=NULL;          in->extension->extension=NULL;
138       }       }
139       Finley_IndexList_insertIndex(in->extension,index);       Dudley_IndexList_insertIndex(in->extension,index);
140    } else {    } else {
141       /* insert index into in->index*/       /* insert index into in->index*/
142       in->index[in->n]=index;       in->index[in->n]=index;
# Line 95  void Finley_IndexList_insertIndex(Finley Line 144  void Finley_IndexList_insertIndex(Finley
144    }    }
145  }  }
146    
147  /* counts the number of row indices in the Finley_IndexList in */  /* counts the number of row indices in the Dudley_IndexList in */
148    
149  dim_t Finley_IndexList_count(Finley_IndexList* in) {  dim_t Dudley_IndexList_count(Dudley_IndexList* in, index_t range_min,index_t range_max) {
150      dim_t i;
151      dim_t out=0;
152      register index_t itmp;
153    if (in==NULL) {    if (in==NULL) {
154       return 0;       return 0;
155    } else {    } else {
156       return (in->n)+Finley_IndexList_count(in->extension);      for (i=0;i<in->n;i++) {
157              itmp=in->index[i];
158              if ((itmp>=range_min) && (range_max>itmp)) ++out;
159        }
160         return out+Dudley_IndexList_count(in->extension, range_min,range_max);
161    }    }
162  }  }
163    
164  /* count the number of row indices in the Finley_IndexList in */  /* count the number of row indices in the Dudley_IndexList in */
165    
166  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) {
167    dim_t i;    dim_t i, ptr;
168      register index_t itmp;
169    if (in!=NULL) {    if (in!=NULL) {
170      for (i=0;i<in->n;i++) array[i]=in->index[i];      ptr=0;
171      Finley_IndexList_toArray(in->extension,&(array[in->n]));      for (i=0;i<in->n;i++) {
172              itmp=in->index[i];
173              if ((itmp>=range_min) && (range_max>itmp)) {
174                 array[ptr]=itmp+index_offset;
175                 ptr++;
176              }
177    
178        }
179        Dudley_IndexList_toArray(in->extension,&(array[ptr]), range_min, range_max, index_offset);
180    }    }
181  }  }
182    
183  /* deallocates the Finley_IndexList in by recursive calls */  /* deallocates the Dudley_IndexList in by recursive calls */
184    
185  void Finley_IndexList_free(Finley_IndexList* in) {  void Dudley_IndexList_free(Dudley_IndexList* in) {
186    if (in!=NULL) {    if (in!=NULL) {
187      Finley_IndexList_free(in->extension);      Dudley_IndexList_free(in->extension);
188      TMPMEMFREE(in);      TMPMEMFREE(in);
189    }    }
190  }  }
191    
192  /*  /* creates a Paso_pattern from a range of indices */
193   * $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)
194   * Revision 1.6  2005/09/15 03:44:22  jgs  {
195   * Merge of development branch dev-02 back to main trunk on 2005-09-15     dim_t *ptr=NULL;
196   *     register dim_t s,i,itmp;
197   * Revision 1.5.2.1  2005/09/07 06:26:18  gross     index_t *index=NULL;
198   * the solver from finley are put into the standalone package paso now     Paso_Pattern* out=NULL;
199   *  
200   * Revision 1.5  2005/07/08 04:07:51  jgs     ptr=MEMALLOC(n+1-n0,index_t);
201   * Merge of development branch back to main trunk on 2005-07-08     if (! Dudley_checkPtr(ptr) ) {
202   *         /* get the number of connections per row */
203   * Revision 1.4  2004/12/15 07:08:32  jgs         #pragma omp parallel for schedule(static) private(i)
204   * *** empty log message ***         for(i=n0;i<n;++i) {
205   * 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);
206   * some changes towards 64 integers in finley         }
207   *         /* accumulate ptr */
208   * Revision 1.1.1.1.2.2  2004/11/24 01:37:13  gross         s=0;
209   * some changes dealing with the integer overflow in memory allocation. Finley solves 4M unknowns now         for(i=n0;i<n;++i) {
210   *                 itmp=ptr[i-n0];
211   *                 ptr[i-n0]=s;
212   *                 s+=itmp;
213   */         }
214           ptr[n-n0]=s;
215           /* fill index */
216           index=MEMALLOC(ptr[n-n0],index_t);
217           if (! Dudley_checkPtr(index)) {
218                  #pragma omp parallel for schedule(static)
219                  for(i=n0;i<n;++i) {
220                      Dudley_IndexList_toArray(&index_list[i],&index[ptr[i-n0]],range_min,range_max,index_offset);
221                  }
222                  out=Paso_Pattern_alloc(PATTERN_FORMAT_DEFAULT,n-n0,range_max+index_offset,ptr,index);
223           }
224      }
225      if (! Dudley_noError()) {
226            MEMFREE(ptr);
227            MEMFREE(index);
228            Paso_Pattern_free(out);
229      }
230      return out;
231    }

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

  ViewVC Help
Powered by ViewVC 1.1.26