/[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/esys2/finley/src/finleyC/IndexList.c revision 100 by jgs, Wed Dec 15 03:48:48 2004 UTC trunk/finley/src/IndexList.c revision 969 by ksteube, Tue Feb 13 23:02:23 2007 UTC
# Line 1  Line 1 
1  /* $Id$ */  /*
2     ************************************************************
3     *          Copyright 2006 by ACcESS MNRF                   *
4     *                                                          *
5     *              http://www.access.edu.au                    *
6     *       Primary Business: Queensland, Australia            *
7     *  Licensed under the Open Software License version 3.0    *
8     *     http://www.opensource.org/licenses/osl-3.0.php       *
9     *                                                          *
10     ************************************************************
11    */
12    
13  /**************************************************************/  /**************************************************************/
14    
# Line 6  Line 16 
16    
17  /**************************************************************/  /**************************************************************/
18    
19  /* Copyrights by ACcESS Australia 2003,2004 */  /*  Author: gross@access.edu.au */
20  /* Author: gross@access.edu.au */  /*  Version: $Id$ */
21    
22  /**************************************************************/  /**************************************************************/
23    
 #include "Finley.h"  
 #include "ElementFile.h"  
 #include "System.h"  
24  #include "IndexList.h"  #include "IndexList.h"
25    
26    /* Translate from distributed/local array indices to global indices */
27    
28    #ifdef PASO_MPI
29    int Finley_IndexList_localToGlobal(Finley_NodeDistribution *dofDistribution, int my_CPU, int localIndex) {
30      /*
31        get global id of icol
32        if icol is internal node (on this CPU): use icol+vtxdist[my_CPU]
33        else use indexExternal[icol-numLocal] to get global index of node
34        (actually DOF...the NodeDistribution structure should have been called DofDistribution)
35      */
36      if (localIndex < dofDistribution->numLocal) {
37        localIndex = localIndex + dofDistribution->vtxdist[my_CPU];
38      }
39      else {
40        localIndex = dofDistribution->indexExternal[localIndex-dofDistribution->numLocal];
41      }
42      return(localIndex);
43    }
44    #endif
45    
46  /**************************************************************/  /**************************************************************/
47  /* inserts the contributions from the element matrices of elements  /* inserts the contributions from the element matrices of elements
48     into the row index col. If symmetric is set, only the upper     into the row index col. If symmetric is set, only the upper
49     triangle of the matrix is stored. */     triangle of the matrix is stored. */
50    
51  void Finley_IndexList_insertElements(Finley_IndexList* index_list, Finley_ElementFile* elements,  void Finley_IndexList_insertElements(Finley_IndexList* index_list, Finley_Mesh* mesh, Finley_ElementFile* elements,
52                                         int reduce_row_order, int packed_row_block_size, maybelong* row_Label,                                         bool_t reduce_row_order, index_t* row_Label,
53                                         int reduce_col_order, int packed_col_block_size, maybelong* col_Label,                                         bool_t reduce_col_order, index_t* col_Label) {
54                                         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 */
55    maybelong e,kr,ir,kc,ic,NN_row,NN_col,i;    index_t color, num_CPUs = 1, my_CPU = 0;
56    Finley_IndexList * tmp_list;    dim_t e,kr,kc,NN_row,NN_col,i,icol,irow;
57    maybelong jr,jc,icol,irow,color;  #ifdef PASO_MPI
58        num_CPUs = mesh->MPIInfo->size;
59        my_CPU = mesh->MPIInfo->rank;
60    #endif
61      /* print_mesh_statistics( mesh, TRUE ); */
62    
63    if (elements!=NULL) {    if (elements!=NULL) {
64      maybelong NN=elements->ReferenceElement->Type->numNodes;      dim_t NN=elements->ReferenceElement->Type->numNodes;
65      maybelong id[NN],*row_node,*col_node;      index_t id[NN],*row_node,*col_node;
66      for (i=0;i<NN;i++) id[i]=i;      for (i=0;i<NN;i++) id[i]=i;
67      if (reduce_col_order) {      if (reduce_col_order) {
68         col_node=elements->ReferenceElement->Type->linearNodes;         col_node=elements->ReferenceElement->Type->linearNodes;
# Line 48  void Finley_IndexList_insertElements(Fin Line 78  void Finley_IndexList_insertElements(Fin
78         row_node=id;         row_node=id;
79         NN_row=elements->ReferenceElement->Type->numNodes;         NN_row=elements->ReferenceElement->Type->numNodes;
80      }      }
81        if (num_CPUs == 1) {
82      #pragma omp parallel private(color)      for (color=elements->minColor;color<=elements->maxColor;color++) {
83      {          #pragma omp for private(e,irow,kr,kc,icol) schedule(static)
84         switch(matType) {          for (e=0;e<elements->numElements;e++) {
85         case CSR:              if (elements->Color[e]==color) {
86           if (symmetric) {                  for (kr=0;kr<NN_row;kr++) {
87              for (color=0;color<elements->numColors;color++) {                    irow=row_Label[elements->Nodes[INDEX2(row_node[kr],e,NN)]];
88                 #pragma omp for private(e,kr,jr,kc,jc,ir,irow,ic,icol,tmp_list) schedule(static)                    for (kc=0;kc<NN_col;kc++) {
89                 for (e=0;e<elements->numElements;e++) {                         icol=col_Label[elements->Nodes[INDEX2(col_node[kc],e,NN)]];
90                    if (elements->Color[e]==color) {                         Finley_IndexList_insertIndex(&(index_list[irow]),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;  
                            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) {  
                      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;  
                            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);  
                               }  
                            }  
                         }  
                      }  
91                    }                    }
92                 }                  }
93              }              }
94           } else {          }
95             for (color=0;color<elements->numColors;color++) {        }
96                 #pragma omp for private (e,kr,jr,kc,jc,ir,irow,ic,icol,tmp_list) schedule(static)      }
97                 for (e=0;e<elements->numElements;e++) {      else {  /* More than one CPU (what's below should also work for one CPU, but let's build confidence in it first) */
98                    if (elements->Color[e]==color) {  #ifdef PASO_MPI
99                       for (kc=0;kc<NN_col;kc++) {      Finley_NodeDistribution *row_degreeOfFreedomDistribution;
100                          jc=col_Label[elements->Nodes[INDEX2(col_node[kc],e,NN)]];      Finley_NodeDistribution *col_degreeOfFreedomDistribution;
101                          for (ic=0;ic<packed_col_block_size;ic++) {      if (reduce_col_order) {
102                             icol=ic+packed_col_block_size*jc;        col_degreeOfFreedomDistribution = mesh->Nodes->reducedDegreeOfFreedomDistribution;
103                             tmp_list=&(index_list[icol]);      }
104                             for (kr=0;kr<NN_row;kr++) {      else {
105                                jr=row_Label[elements->Nodes[INDEX2(row_node[kr],e,NN)]];        col_degreeOfFreedomDistribution = mesh->Nodes->degreeOfFreedomDistribution;
                               for (ir=0;ir<packed_row_block_size;ir++) {  
                                  irow=ir+packed_row_block_size*jr;  
                                  Finley_IndexList_insertIndex(tmp_list,irow);  
                               }  
                            }  
                         }  
                      }  
                   }  
                }  
             }  
          } /* if else */  
        } /* switch matType */  
106      }      }
107        if (reduce_row_order) {
108          row_degreeOfFreedomDistribution = mesh->Nodes->reducedDegreeOfFreedomDistribution;
109        }
110        else {
111          row_degreeOfFreedomDistribution = mesh->Nodes->degreeOfFreedomDistribution;
112        }
113        /* Not using loop over colors as above */ {
114            #pragma omp for private(e,irow,kr,kc,icol) schedule(static)
115        for (e=0;e<elements->numElements;e++) {
116                    for (kr=0;kr<NN_row;kr++) {
117                      irow=row_Label[elements->Nodes[INDEX2(row_node[kr],e,NN)]];
118              if (irow < row_degreeOfFreedomDistribution->numLocal) {
119                        for (kc=0;kc<NN_col;kc++) {
120                   /* Get the local col ID */
121                           icol=col_Label[elements->Nodes[INDEX2(col_node[kc],e,NN)]];
122                   /* Convert to global col ID (row ID is saved as local value) */
123                   icol = Finley_IndexList_localToGlobal(col_degreeOfFreedomDistribution, my_CPU, icol);
124                           Finley_IndexList_insertIndex(&(index_list[irow]),icol);
125                   printf("ksteube Finley_IndexList_insertIndex cpu= %d irow= %d icol= %d\n", my_CPU, irow, icol);
126                        }
127              }
128                    }
129            }
130          }
131    #endif
132        }   /* More than one CPU */
133    }    }
134    return;    return;
135  }  }
136    
137  /* 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 */
138    
139  void Finley_IndexList_insertIndex(Finley_IndexList* in, maybelong index) {  void Finley_IndexList_insertIndex(Finley_IndexList* in, index_t index) {
140    int i;    dim_t i;
141    /* is index in in? */    /* is index in in? */
142    for (i=0;i<in->n;i++) {    for (i=0;i<in->n;i++) {
143      if (in->index[i]==index) return;      if (in->index[i]==index)  return;
144    }    }
145    /* index could not be found */    /* index could not be found */
146    if (in->n==INDEXLIST_LENGTH) {    if (in->n==INDEXLIST_LENGTH) {
147       /* if in->index is full check the extension */       /* if in->index is full check the extension */
148       if (in->extension==NULL) {       if (in->extension==NULL) {
149          in->extension=(Finley_IndexList*) TMPMEMALLOC(sizeof(Finley_IndexList));          in->extension=TMPMEMALLOC(1,Finley_IndexList);
150          if (Finley_checkPtr(in->extension)) return;          if (Finley_checkPtr(in->extension)) return;
151          in->extension->n=0;          in->extension->n=0;
152          in->extension->extension=NULL;          in->extension->extension=NULL;
# Line 177  void Finley_IndexList_insertIndex(Finley Line 161  void Finley_IndexList_insertIndex(Finley
161    
162  /* counts the number of row indices in the Finley_IndexList in */  /* counts the number of row indices in the Finley_IndexList in */
163    
164  int Finley_IndexList_count(Finley_IndexList* in) {  dim_t Finley_IndexList_count(Finley_IndexList* in) {
165    if (in==NULL) {    if (in==NULL) {
166       return 0;       return 0;
167    } else {    } else {
# Line 187  int Finley_IndexList_count(Finley_IndexL Line 171  int Finley_IndexList_count(Finley_IndexL
171    
172  /* count the number of row indices in the Finley_IndexList in */  /* count the number of row indices in the Finley_IndexList in */
173    
174  void Finley_IndexList_toArray(Finley_IndexList* in, maybelong* array) {  void Finley_IndexList_toArray(Finley_IndexList* in, index_t* array) {
175    int i;    dim_t i;
176    if (in!=NULL) {    if (in!=NULL) {
177      for (i=0;i<in->n;i++) array[i]=in->index[i]+INDEX_OFFSET;      for (i=0;i<in->n;i++) array[i]=in->index[i];
178      Finley_IndexList_toArray(in->extension,&(array[in->n]));      Finley_IndexList_toArray(in->extension,&(array[in->n]));
179    }    }
180  }  }
# Line 206  void Finley_IndexList_free(Finley_IndexL Line 190  void Finley_IndexList_free(Finley_IndexL
190    
191  /*  /*
192   * $Log$   * $Log$
193   * Revision 1.3  2004/12/15 03:48:45  jgs   * Revision 1.6  2005/09/15 03:44:22  jgs
194   * *** empty log message ***   * Merge of development branch dev-02 back to main trunk on 2005-09-15
195   *   *
196   * Revision 1.1.1.1  2004/10/26 06:53:57  jgs   * Revision 1.5.2.1  2005/09/07 06:26:18  gross
197   * initial import of project esys2   * the solver from finley are put into the standalone package paso now
198   *   *
199   * Revision 1.1.2.2  2004/10/26 06:36:39  jgs   * Revision 1.5  2005/07/08 04:07:51  jgs
200   * committing Lutz's changes to branch jgs   * Merge of development branch back to main trunk on 2005-07-08
201     *
202     * Revision 1.4  2004/12/15 07:08:32  jgs
203     * *** empty log message ***
204     * Revision 1.1.1.1.2.3  2005/06/29 02:34:50  gross
205     * some changes towards 64 integers in finley
206   *   *
207   * Revision 1.2  2004/10/13 01:53:42  gross   * Revision 1.1.1.1.2.2  2004/11/24 01:37:13  gross
208   * bug in CSC assembling fixed   * some changes dealing with the integer overflow in memory allocation. Finley solves 4M unknowns now
209   *   *
  * Revision 1.1  2004/07/02 04:21:13  gross  
  * Finley C code has been included  
210   *   *
211   *   *
212   */   */

Legend:
Removed from v.100  
changed lines
  Added in v.969

  ViewVC Help
Powered by ViewVC 1.1.26