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

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

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

trunk/esys2/finley/src/finleyC/Util.c revision 147 by jgs, Fri Aug 12 01:45:47 2005 UTC branches/domexper/dudley/src/Util.c revision 3080 by jfenwick, Tue Aug 3 04:28:03 2010 UTC
# Line 1  Line 1 
 /* $Id$ */  
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    
 /*   Some utility routines: */  
14    
15  /**************************************************************/  /**************************************************************/
16    
17  /*   Copyrights by ACcESS Australia, 2003 */  /*   Some utility routines: */
 /*   author: gross@access.edu.au */  
 /*   Version: $Id$ */  
18    
19  /**************************************************************/  /**************************************************************/
20    
 #include "Common.h"  
21  #include "Finley.h"  #include "Finley.h"
22  #include "Util.h"  #include "Util.h"
23    
24  #ifdef _OPENMP  #ifdef _OPENMP
25  #include <omp.h>  #include <omp.h>
26  #endif  #endif
# Line 63  void Finley_Util_Gather_int(dim_t len,in Line 69  void Finley_Util_Gather_int(dim_t len,in
69    
70  /*   adds a vector in into out using and index. */  /*   adds a vector in into out using and index. */
71    
72  /*        out(1:numData,index(1:len))+=in(1:numData,1:len) */  /*        out(1:numData,index[p])+=in(1:numData,p) where p = {k=1...len , index[k]<upperBound}*/
73    
74    
75  void Finley_Util_AddScatter(dim_t len,index_t* index,dim_t numData,double* in,double * out){  void Finley_Util_AddScatter(dim_t len,index_t* index,dim_t numData,double* in,double * out, index_t upperBound){
76     dim_t i,s;     dim_t i,s;
77     for (s=0;s<len;s++) {     for (s=0;s<len;s++) {
78         for(i=0;i<numData;i++) {         for(i=0;i<numData;i++) {
79            #pragma omp atomic            if( index[s]<upperBound ) {
80            out[INDEX2(i,index[s],numData)]+=in[INDEX2(i,s,numData)];              out[INDEX2(i,index[s],numData)]+=in[INDEX2(i,s,numData)];
81          }
82         }         }
83     }     }
84  }  }
# Line 81  void Finley_Util_AddScatter(dim_t len,in Line 89  void Finley_Util_AddScatter(dim_t len,in
89    
90  void Finley_Util_SmallMatMult(dim_t A1,dim_t A2, double* A, dim_t B2, double*B, double* C) {  void Finley_Util_SmallMatMult(dim_t A1,dim_t A2, double* A, dim_t B2, double*B, double* C) {
91      dim_t i,j,s;      dim_t i,j,s;
92      for (i=0;i<A1*A2;i++) A[i]=0;      register double rtmp;
93         for (i=0;i<A1;i++) {         for (i=0;i<A1;i++) {
94            for (j=0;j<A2;j++) {            for (j=0;j<A2;j++) {
95               for (s=0;s<B2;s++) {               rtmp=0;
96                  A[INDEX2(i,j,A1)]+=B[INDEX2(i,s,A1)]*C[INDEX2(s,j,B2)];               for (s=0;s<B2;s++) rtmp+=B[INDEX2(i,s,A1)]*C[INDEX2(s,j,B2)];
97               }               A[INDEX2(i,j,A1)]=rtmp;
98            }            }
99         }         }
100  }  }
# Line 97  void Finley_Util_SmallMatMult(dim_t A1,d Line 105  void Finley_Util_SmallMatMult(dim_t A1,d
105    
106  void Finley_Util_SmallMatSetMult(dim_t len,dim_t A1,dim_t A2, double* A, dim_t B2, double*B, double* C) {  void Finley_Util_SmallMatSetMult(dim_t len,dim_t A1,dim_t A2, double* A, dim_t B2, double*B, double* C) {
107      dim_t q,i,j,s;      dim_t q,i,j,s;
108      for (i=0;i<A1*A2*len;i++) A[i]=0;      register double rtmp;
109      for (q=0;q<len;q++) {      for (q=0;q<len;q++) {
110         for (i=0;i<A1;i++) {         for (i=0;i<A1;i++) {
111            for (j=0;j<A2;j++) {            for (j=0;j<A2;j++) {
112               for (s=0;s<B2;s++) {               rtmp=0;
113                  A[INDEX3(i,j,q,A1,A2)]+=B[INDEX3(i,s,q,A1,B2)]*C[INDEX3(s,j,q,B2,A2)];               for (s=0;s<B2;s++) rtmp+=B[INDEX3(i,s,q,A1,B2)]*C[INDEX3(s,j,q,B2,A2)];
114               }               A[INDEX3(i,j,q, A1,A2)]=rtmp;
115              }
116           }
117        }
118    }
119    /*    multiplies a set of matries with a single matrix: */
120    
121    /*        A(1:A1,1:A2,i)=B(1:A1,1:B2,i)*C(1:B2,1:A2) i=1,len */
122    
123    void Finley_Util_SmallMatSetMult1(dim_t len,dim_t A1,dim_t A2, double* A, dim_t B2, double*B, double* C) {
124        dim_t q,i,j,s;
125        register double rtmp;
126        for (q=0;q<len;q++) {
127           for (i=0;i<A1;i++) {
128              for (j=0;j<A2;j++) {
129                 rtmp=0;
130                 for (s=0;s<B2;s++) rtmp+=B[INDEX3(i,s,q, A1,B2)]*C[INDEX2(s,j,B2)];
131                 A[INDEX3(i,j,q,A1,A2)]=rtmp;
132            }            }
133         }         }
134      }      }
# Line 124  void Finley_Util_InvertSmallMat(dim_t le Line 149  void Finley_Util_InvertSmallMat(dim_t le
149                 D=1./D;                 D=1./D;
150                 invA[q]=D;                 invA[q]=D;
151              } else {              } else {
152                 Finley_ErrorCode=ZERO_DIVISION_ERROR;                 Finley_setError(ZERO_DIVISION_ERROR,"__FILE__: Non-regular matrix");
                sprintf(Finley_ErrorMsg,"Non-regular matrix");  
153                 return;                 return;
154              }              }
155           }           }
# Line 147  void Finley_Util_InvertSmallMat(dim_t le Line 171  void Finley_Util_InvertSmallMat(dim_t le
171                 invA[INDEX3(0,1,q,2,2)]=-A12*D;                 invA[INDEX3(0,1,q,2,2)]=-A12*D;
172                 invA[INDEX3(1,1,q,2,2)]= A11*D;                 invA[INDEX3(1,1,q,2,2)]= A11*D;
173              } else {              } else {
174                 Finley_ErrorCode=ZERO_DIVISION_ERROR;                 Finley_setError(ZERO_DIVISION_ERROR,"__FILE__: Non-regular matrix");
                sprintf(Finley_ErrorMsg,"Non-regular matrix");  
175                 return;                 return;
176              }              }
177           }           }
# Line 180  void Finley_Util_InvertSmallMat(dim_t le Line 203  void Finley_Util_InvertSmallMat(dim_t le
203                 invA[INDEX3(1,2,q,3,3)]=(A13*A21-A11*A23)*D;                 invA[INDEX3(1,2,q,3,3)]=(A13*A21-A11*A23)*D;
204                 invA[INDEX3(2,2,q,3,3)]=(A11*A22-A12*A21)*D;                 invA[INDEX3(2,2,q,3,3)]=(A11*A22-A12*A21)*D;
205              } else {              } else {
206                 Finley_ErrorCode=ZERO_DIVISION_ERROR;                 Finley_setError(ZERO_DIVISION_ERROR,"__FILE__: Non-regular matrix");
                sprintf(Finley_ErrorMsg,"Non-regular matrix");  
207                 return;                 return;
208              }              }
209           }           }
# Line 251  void  Finley_NormalVector(dim_t len, dim Line 273  void  Finley_NormalVector(dim_t len, dim
273              A21=A[INDEX3(1,0,q,2,dim1)];              A21=A[INDEX3(1,0,q,2,dim1)];
274              length = sqrt(A11*A11+A21*A21);              length = sqrt(A11*A11+A21*A21);
275              if (! length>0) {              if (! length>0) {
276                 Finley_ErrorCode=ZERO_DIVISION_ERROR;                 Finley_setError(ZERO_DIVISION_ERROR,"__FILE__: area equals zero.");
                sprintf(Finley_ErrorMsg,"area equals zero.");  
277                 return;                 return;
278              } else {              } else {
279                 invlength=1./length;                 invlength=1./length;
# Line 274  void  Finley_NormalVector(dim_t len, dim Line 295  void  Finley_NormalVector(dim_t len, dim
295              CO_A33=A11*A22-A21*A12;              CO_A33=A11*A22-A21*A12;
296              length=sqrt(CO_A13*CO_A13+CO_A23*CO_A23+CO_A33*CO_A33);              length=sqrt(CO_A13*CO_A13+CO_A23*CO_A23+CO_A33*CO_A33);
297              if (! length>0) {              if (! length>0) {
298                 Finley_ErrorCode=ZERO_DIVISION_ERROR;                 Finley_setError(ZERO_DIVISION_ERROR,"__FILE__: area equals zero.");
                sprintf(Finley_ErrorMsg,"area equals zero.");  
299                 return;                 return;
300              } else {              } else {
301                 invlength=1./length;                 invlength=1./length;
# Line 349  int Finley_Util_ValueAndIndex_compar(con Line 369  int Finley_Util_ValueAndIndex_compar(con
369     e2=(Finley_Util_ValueAndIndex*) arg2;     e2=(Finley_Util_ValueAndIndex*) arg2;
370     if (e1->value < e2->value) return -1;     if (e1->value < e2->value) return -1;
371     if (e1->value > e2->value) return  1;     if (e1->value > e2->value) return  1;
372       if (e1->index < e2->index) return -1;
373       if (e1->index > e2->index) return  1;
374     return 0;     return 0;
375  }  }
376    
377  void Finley_Util_sortValueAndIndex(dim_t n,Finley_Util_ValueAndIndex* array) {  void Finley_Util_sortValueAndIndex(dim_t n,Finley_Util_ValueAndIndex* array) {
378       /* OMP : needs parallelization !*/       /* OMP : needs parallelization !*/
379       qsort(array,n,sizeof(Finley_Util_ValueAndIndex),Finley_Util_ValueAndIndex_compar);       qsort(array,n,sizeof(Finley_Util_ValueAndIndex),Finley_Util_ValueAndIndex_compar);
# Line 402  index_t Finley_Util_getMaxInt(dim_t dim, Line 425  index_t Finley_Util_getMaxInt(dim_t dim,
425     }     }
426     return out;     return out;
427  }  }
428    /**************************************************************/
429    
430    /* calculates the minimum value from a dim X N integer array */
431    
432    index_t Finley_Util_getFlaggedMinInt(dim_t dim,dim_t N,index_t* values, index_t ignore) {
433       dim_t i,j;
434       index_t out,out_local;
435       out=INDEX_T_MAX;
436       if (values!=NULL && dim*N>0 ) {
437         out=values[0];
438         #pragma omp parallel private(out_local)
439         {
440             out_local=out;
441             #pragma omp for private(i,j) schedule(static)
442             for (j=0;j<N;j++) {
443               for (i=0;i<dim;i++) if (values[INDEX2(i,j,dim)]!=ignore) out_local=MIN(out_local,values[INDEX2(i,j,dim)]);
444             }
445             #pragma omp critical
446             out=MIN(out_local,out);
447         }
448       }
449       return out;
450    }
451                                                                                                                                                      
452    /* calculates the maximum value from a dim X N integer array */
453    
454    index_t Finley_Util_getFlaggedMaxInt(dim_t dim,dim_t N,index_t* values, index_t ignore) {
455       dim_t i,j;
456       index_t out,out_local;
457       out=-INDEX_T_MAX;
458       if (values!=NULL && dim*N>0 ) {
459         out=values[0];
460         #pragma omp parallel private(out_local)
461         {
462             out_local=out;
463             #pragma omp for private(i,j) schedule(static)
464             for (j=0;j<N;j++) {
465                 for (i=0;i<dim;i++) if (values[INDEX2(i,j,dim)]!=ignore) out_local=MAX(out_local,values[INDEX2(i,j,dim)]);
466             }
467             #pragma omp critical
468             out=MAX(out_local,out);
469          }
470       }
471       return out;
472    }
473    
474  /* set the index of the positive entries in mask. The length of index is returned. */  /* set the index of the positive entries in mask. The length of index is returned. */
475    
# Line 431  index_t Finley_Util_cumsum(dim_t N,index Line 499  index_t Finley_Util_cumsum(dim_t N,index
499     index_t out=0,tmp;     index_t out=0,tmp;
500     dim_t i;     dim_t i;
501     #ifdef _OPENMP     #ifdef _OPENMP
502        index_t partial_sums[omp_get_max_threads()],sum;        index_t *partial_sums=NULL, sum;
503          partial_sums=TMPMEMALLOC(omp_get_max_threads(),index_t);
504        #pragma omp parallel private(sum,i,tmp)        #pragma omp parallel private(sum,i,tmp)
505        {        {
506          sum=0;          sum=0;
# Line 457  index_t Finley_Util_cumsum(dim_t N,index Line 526  index_t Finley_Util_cumsum(dim_t N,index
526            array[i]=tmp;            array[i]=tmp;
527          }          }
528        }        }
529          TMPMEMFREE(partial_sums);
530     #else     #else
531        for (i=0;i<N;++i) {        for (i=0;i<N;++i) {
532           tmp=out;           tmp=out;
# Line 466  index_t Finley_Util_cumsum(dim_t N,index Line 536  index_t Finley_Util_cumsum(dim_t N,index
536     #endif     #endif
537     return out;     return out;
538  }  }
539    void Finley_Util_setValuesInUse(const index_t *values, const dim_t numValues, dim_t *numValuesInUse, index_t **valuesInUse, Paso_MPIInfo* mpiinfo)
540    {
541       dim_t i;
542       index_t lastFoundValue=INDEX_T_MIN, minFoundValue, local_minFoundValue, *newValuesInUse=NULL;
543       register index_t itmp;
544       bool_t allFound=FALSE;
545       dim_t nv=0;
546    
547       while (! allFound) {
548           /*
549            *  find smallest value bigger than lastFoundValue
550            */
551            minFoundValue=INDEX_T_MAX;
552            #pragma omp parallel private(local_minFoundValue)
553            {
554                local_minFoundValue=minFoundValue;
555                #pragma omp for private(i,itmp) schedule(static)
556                for (i=0;i< numValues;i++) {
557                   itmp=values[i];
558                   if ((itmp>lastFoundValue) && (itmp<local_minFoundValue)) local_minFoundValue=itmp;
559                }
560                #pragma omp critical
561                {
562                   if (local_minFoundValue<minFoundValue) minFoundValue=local_minFoundValue;
563                }
564    
565             }
566             #ifdef PASO_MPI
567             local_minFoundValue=minFoundValue;
568             MPI_Allreduce(&local_minFoundValue,&minFoundValue, 1, MPI_INT, MPI_MIN, mpiinfo->comm );
569             #endif
570             /* if we found a new tag we need to add this too the valuesInUseList */
571    
572             if (minFoundValue < INDEX_T_MAX) {
573                 newValuesInUse=MEMALLOC(nv+1,index_t);
574                 if (*valuesInUse!=NULL) {
575                     memcpy(newValuesInUse,*valuesInUse,sizeof(index_t)*nv);
576                     MEMFREE(*valuesInUse);
577                 }
578                 newValuesInUse[nv]=minFoundValue;
579                 *valuesInUse=newValuesInUse;
580                 newValuesInUse=NULL;
581                 nv++;
582                 lastFoundValue=minFoundValue;
583             } else {
584                 allFound=TRUE;
585             }
586       }
587       *numValuesInUse=nv;
588    }
589    
590    
591  void Finley_copyDouble(dim_t n,double* source, double* target) {  #ifdef PASO_MPI
592    dim_t i;  void Finley_printDoubleArray( FILE *fid, dim_t n, double *array, char *name  )
593    for (i=0;i<n;i++) target[i]=source[i];  {
594  }    index_t i;
595      
596  /*    if( name )
597   * $Log$      fprintf( fid, "%s [ ", name );
598   * Revision 1.8  2005/08/12 01:45:43  jgs    else
599   * erge of development branch dev-02 back to main trunk on 2005-08-12      fprintf( fid, "[ " );  
600   *    for( i=0; i<(n<60 ? n : 60); i++ )
601   * Revision 1.7.2.1  2005/08/04 22:41:11  gross      fprintf( fid, "%g ", array[i] );
602   * some extra routines for finley that might speed-up RHS assembling in some cases (not actived right now)    if( n>=30 )
603   *      fprintf( fid, "... " );
604   * Revision 1.7  2005/07/08 04:07:59  jgs    fprintf( fid, "]\n" );
605   * Merge of development branch back to main trunk on 2005-07-08  }
606   *  void Finley_printIntArray( FILE *fid, dim_t n, int *array, char *name  )
607   * Revision 1.1.1.1.2.4  2005/06/29 02:34:57  gross  {
608   * some changes towards 64 integers in finley    index_t i;
609   *    
610   * Revision 1.1.1.1.2.3  2005/03/02 23:35:06  gross    if( name )
611   * reimplementation of the ILU in Finley. block size>1 still needs some testing      fprintf( fid, "%s [ ", name );
612   *    else
613   * Revision 1.1.1.1.2.2  2005/02/18 02:27:31  gross      fprintf( fid, "[ " );  
614   * two function that will be used for a reimplementation of the ILU preconditioner    for( i=0; i<(n<60 ? n : 60); i++ )
615   *      fprintf( fid, "%d ", array[i] );
616   * Revision 1.1.1.1.2.1  2004/11/12 06:58:19  gross    if( n>=30 )
617   * a lot of changes to get the linearPDE class running: most important change is that there is no matrix format exposed to the user anymore. the format is chosen by the Domain according to the solver and symmetry      fprintf( fid, "... " );
618   *    fprintf( fid, "]\n" );
619   * Revision 1.1.1.1  2004/10/26 06:53:57  jgs  }
620   * initial import of project esys2  void Finley_printMaskArray( FILE *fid, dim_t n, int *array, char *name  )
621   *  {
622   * Revision 1.3  2004/08/26 12:03:52  gross    index_t i;
623   * Some other bug in Finley_Assemble_gradient fixed.    
624   *    if( name )
625   * Revision 1.2  2004/07/02 04:21:13  gross      fprintf( fid, "%s [ ", name );
626   * Finley C code has been included    else
627   *      fprintf( fid, "[ " );  
628   * Revision 1.1.1.1  2004/06/24 04:00:40  johng    for( i=0; i<(n<60 ? n : 60); i++ )
629   * Initial version of eys using boost-python.      if( array[i]!=-1 )
630   *        fprintf( fid, "%3d ", array[i] );
631   *      else
632   */        fprintf( fid, "  * " );
633      if( n>=30 )
634        fprintf( fid, "... " );
635      fprintf( fid, "]\n" );
636    }
637    #endif

Legend:
Removed from v.147  
changed lines
  Added in v.3080

  ViewVC Help
Powered by ViewVC 1.1.26