/[escript]/branches/domexper/dudley/src/IndexList.c
ViewVC logotype

Diff of /branches/domexper/dudley/src/IndexList.c

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

trunk/esys2/finley/src/finleyC/IndexList.c revision 100 by jgs, Wed Dec 15 03:48:48 2004 UTC trunk/finley/src/IndexList.c revision 1722 by gross, Fri Aug 22 04:20:30 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
11     *  Licensed under the Open Software License version 3.0
12     *     http://www.opensource.org/licenses/osl-3.0.php
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           }
70           TMPMEMFREE(id);
71        }
72      }
73      return;
74    }
75    
76    void Finley_IndexList_insertElementsWithRowRange(Finley_IndexList* index_list, index_t firstRow, index_t lastRow,
77                                                     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           } else {      }
100             for (color=0;color<elements->numColors;color++) {    }
101                 #pragma omp for private (e,kr,jr,kc,jc,ir,irow,ic,icol,tmp_list) schedule(static)  }
102                 for (e=0;e<elements->numElements;e++) {  void Finley_IndexList_insertElementsWithRowRangeNoMainDiagonal(Finley_IndexList* index_list, index_t firstRow, index_t lastRow,
103                    if (elements->Color[e]==color) {                                                                Finley_ElementFile* elements, index_t* row_map, index_t* col_map)
104                       for (kc=0;kc<NN_col;kc++) {  {
105                          jc=col_Label[elements->Nodes[INDEX2(col_node[kc],e,NN)]];    index_t color;
106                          for (ic=0;ic<packed_col_block_size;ic++) {    dim_t e,kr,kc,icol,irow, NN,irow_loc;
107                             icol=ic+packed_col_block_size*jc;    if (elements!=NULL) {
108                             tmp_list=&(index_list[icol]);      NN=elements->numNodes;
109                             for (kr=0;kr<NN_row;kr++) {      for (color=elements->minColor;color<=elements->maxColor;color++) {
110                                jr=row_Label[elements->Nodes[INDEX2(row_node[kr],e,NN)]];             #pragma omp for private(e,irow,kr,kc,icol,irow_loc) schedule(static)
111                                for (ir=0;ir<packed_row_block_size;ir++) {             for (e=0;e<elements->numElements;e++) {
112                                   irow=ir+packed_row_block_size*jr;                 if (elements->Color[e]==color) {
113                                   Finley_IndexList_insertIndex(tmp_list,irow);                     for (kr=0;kr<NN;kr++) {
114                                }                       irow=row_map[elements->Nodes[INDEX2(kr,e,NN)]];
115                             }                       if ((firstRow<=irow) && (irow < lastRow)) {
116                          }                            irow_loc=irow-firstRow;
117                       }                            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              }             }
          } /* if else */  
        } /* switch matType */  
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 */
130    
131  void Finley_IndexList_insertIndex(Finley_IndexList* in, maybelong index) {  void Finley_IndexList_insertIndex(Finley_IndexList* in, index_t index) {
132    int i;    dim_t i;
133    /* is index in in? */    /* is index in in? */
134    for (i=0;i<in->n;i++) {    for (i=0;i<in->n;i++) {
135      if (in->index[i]==index) return;      if (in->index[i]==index)  return;
136    }    }
137    /* index could not be found */    /* index could not be found */
138    if (in->n==INDEXLIST_LENGTH) {    if (in->n==INDEXLIST_LENGTH) {
139       /* if in->index is full check the extension */       /* if in->index is full check the extension */
140       if (in->extension==NULL) {       if (in->extension==NULL) {
141          in->extension=(Finley_IndexList*) TMPMEMALLOC(sizeof(Finley_IndexList));          in->extension=TMPMEMALLOC(1,Finley_IndexList);
142          if (Finley_checkPtr(in->extension)) return;          if (Finley_checkPtr(in->extension)) return;
143          in->extension->n=0;          in->extension->n=0;
144          in->extension->extension=NULL;          in->extension->extension=NULL;
# Line 177  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  int 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, maybelong* array) {  void Finley_IndexList_toArray(Finley_IndexList* in, index_t* array, index_t range_min,index_t range_max, index_t index_offset) {
174    int 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]+INDEX_OFFSET;      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 204  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.3  2004/12/15 03:48:45  jgs  {
202   * *** empty log message ***     dim_t *ptr=NULL;
203   *     register dim_t s,i,itmp;
204   * Revision 1.1.1.1  2004/10/26 06:53:57  jgs     index_t *index=NULL;
205   * initial import of project esys2     Paso_Pattern* out=NULL;
206   *  
207   * Revision 1.1.2.2  2004/10/26 06:36:39  jgs     ptr=MEMALLOC(n+1-n0,index_t);
208   * committing Lutz's changes to branch jgs     if (! Finley_checkPtr(ptr) ) {
209   *         /* get the number of connections per row */
210   * Revision 1.2  2004/10/13 01:53:42  gross         #pragma omp parallel for schedule(static) private(i)
211   * bug in CSC assembling fixed         for(i=n0;i<n;++i) {
212   *                ptr[i-n0]=Finley_IndexList_count(&index_list[i],range_min,range_max);
213   * Revision 1.1  2004/07/02 04:21:13  gross         }
214   * Finley C code has been included         /* accumulate ptr */
215   *         s=0;
216   *         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.100  
changed lines
  Added in v.1722

  ViewVC Help
Powered by ViewVC 1.1.26