/[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 3114 by jfenwick, Fri Aug 27 05:26:25 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"  
 #include "Finley.h"  
21  #include "Util.h"  #include "Util.h"
22    
23  #ifdef _OPENMP  #ifdef _OPENMP
24  #include <omp.h>  #include <omp.h>
25  #endif  #endif
# Line 23  Line 28 
28    
29  /*   returns true if any of the values in the short array values is not equalt to Zero */  /*   returns true if any of the values in the short array values is not equalt to Zero */
30    
31  bool_t Finley_Util_anyNonZeroDouble(dim_t N, double* values) {  bool_t Dudley_Util_anyNonZeroDouble(dim_t N, double* values) {
32     dim_t q;     dim_t q;
33     for (q=0;q<N;++q) if (ABS(values[q])>0) return TRUE;     for (q=0;q<N;++q) if (ABS(values[q])>0) return TRUE;
34     return FALSE;     return FALSE;
# Line 34  bool_t Finley_Util_anyNonZeroDouble(dim_ Line 39  bool_t Finley_Util_anyNonZeroDouble(dim_
39    
40  /*        out(1:numData,1:len)=in(1:numData,index(1:len)) */  /*        out(1:numData,1:len)=in(1:numData,index(1:len)) */
41    
42  void Finley_Util_Gather_double(dim_t len,index_t* index,dim_t numData,double* in, double * out){  void Dudley_Util_Gather_double(dim_t len,index_t* index,dim_t numData,double* in, double * out){
43      dim_t s,i;      dim_t s,i;
44      for (s=0;s<len;s++) {      for (s=0;s<len;s++) {
45         for (i=0;i<numData;i++) {         for (i=0;i<numData;i++) {
# Line 50  void Finley_Util_Gather_double(dim_t len Line 55  void Finley_Util_Gather_double(dim_t len
55    
56  /*        out(1:numData,1:len)=in(1:numData,index(1:len)) */  /*        out(1:numData,1:len)=in(1:numData,index(1:len)) */
57    
58  void Finley_Util_Gather_int(dim_t len,index_t* index,dim_t numData, index_t* in, index_t * out){  void Dudley_Util_Gather_int(dim_t len,index_t* index,dim_t numData, index_t* in, index_t * out){
59      dim_t s,i;      dim_t s,i;
60      for (s=0;s<len;s++) {      for (s=0;s<len;s++) {
61         for (i=0;i<numData;i++) {         for (i=0;i<numData;i++) {
# Line 63  void Finley_Util_Gather_int(dim_t len,in Line 68  void Finley_Util_Gather_int(dim_t len,in
68    
69  /*   adds a vector in into out using and index. */  /*   adds a vector in into out using and index. */
70    
71  /*        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}*/
72    
73    
74  void Finley_Util_AddScatter(dim_t len,index_t* index,dim_t numData,double* in,double * out){  void Dudley_Util_AddScatter(dim_t len,index_t* index,dim_t numData,double* in,double * out, index_t upperBound){
75     dim_t i,s;     dim_t i,s;
76     for (s=0;s<len;s++) {     for (s=0;s<len;s++) {
77         for(i=0;i<numData;i++) {         for(i=0;i<numData;i++) {
78            #pragma omp atomic            if( index[s]<upperBound ) {
79            out[INDEX2(i,index[s],numData)]+=in[INDEX2(i,s,numData)];              out[INDEX2(i,index[s],numData)]+=in[INDEX2(i,s,numData)];
80          }
81         }         }
82     }     }
83  }  }
# Line 79  void Finley_Util_AddScatter(dim_t len,in Line 86  void Finley_Util_AddScatter(dim_t len,in
86    
87  /*          A(1:A1,1:A2)=B(1:A1,1:B2)*C(1:B2,1:A2) */  /*          A(1:A1,1:A2)=B(1:A1,1:B2)*C(1:B2,1:A2) */
88    
89  void Finley_Util_SmallMatMult(dim_t A1,dim_t A2, double* A, dim_t B2, double*B, double* C) {  void Dudley_Util_SmallMatMult(dim_t A1,dim_t A2, double* A, dim_t B2, double*B, double* C) {
90      dim_t i,j,s;      dim_t i,j,s;
91      for (i=0;i<A1*A2;i++) A[i]=0;      register double rtmp;
92         for (i=0;i<A1;i++) {         for (i=0;i<A1;i++) {
93            for (j=0;j<A2;j++) {            for (j=0;j<A2;j++) {
94               for (s=0;s<B2;s++) {               rtmp=0;
95                  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)];
96               }               A[INDEX2(i,j,A1)]=rtmp;
97            }            }
98         }         }
99  }  }
# Line 95  void Finley_Util_SmallMatMult(dim_t A1,d Line 102  void Finley_Util_SmallMatMult(dim_t A1,d
102    
103  /*        A(1:A1,1:A2,i)=B(1:A1,1:B2,i)*C(1:B2,1:A2,i) i=1,len */  /*        A(1:A1,1:A2,i)=B(1:A1,1:B2,i)*C(1:B2,1:A2,i) i=1,len */
104    
105  void Finley_Util_SmallMatSetMult(dim_t len,dim_t A1,dim_t A2, double* A, dim_t B2, double*B, double* C) {  void Dudley_Util_SmallMatSetMult(dim_t len,dim_t A1,dim_t A2, double* A, dim_t B2, double*B, double* C) {
106      dim_t q,i,j,s;      dim_t q,i,j,s;
107      for (i=0;i<A1*A2*len;i++) A[i]=0;      register double rtmp;
108      for (q=0;q<len;q++) {      for (q=0;q<len;q++) {
109         for (i=0;i<A1;i++) {         for (i=0;i<A1;i++) {
110            for (j=0;j<A2;j++) {            for (j=0;j<A2;j++) {
111               for (s=0;s<B2;s++) {               rtmp=0;
112                  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)];
113               }               A[INDEX3(i,j,q, A1,A2)]=rtmp;
114              }
115           }
116        }
117    }
118    /*    multiplies a set of matries with a single matrix: */
119    
120    /*        A(1:A1,1:A2,i)=B(1:A1,1:B2,i)*C(1:B2,1:A2) i=1,len */
121    
122    void Dudley_Util_SmallMatSetMult1(dim_t len,dim_t A1,dim_t A2, double* A, dim_t B2, double*B, double* C) {
123        dim_t q,i,j,s;
124        register double rtmp;
125        for (q=0;q<len;q++) {
126           for (i=0;i<A1;i++) {
127              for (j=0;j<A2;j++) {
128                 rtmp=0;
129                 for (s=0;s<B2;s++) rtmp+=B[INDEX3(i,s,q, A1,B2)]*C[INDEX2(s,j,B2)];
130                 A[INDEX3(i,j,q,A1,A2)]=rtmp;
131            }            }
132         }         }
133      }      }
# Line 111  void Finley_Util_SmallMatSetMult(dim_t l Line 135  void Finley_Util_SmallMatSetMult(dim_t l
135  /*    inverts the set of dim x dim matrices A(:,:,1:len) with dim=1,2,3 */  /*    inverts the set of dim x dim matrices A(:,:,1:len) with dim=1,2,3 */
136  /*    the determinante is returned. */  /*    the determinante is returned. */
137    
138  void Finley_Util_InvertSmallMat(dim_t len,dim_t dim,double* A,double *invA, double* det){  void Dudley_Util_InvertSmallMat(dim_t len,dim_t dim,double* A,double *invA, double* det){
139     dim_t q;     dim_t q;
140     register double D,A11,A12,A13,A21,A22,A23,A31,A32,A33;     register double D,A11,A12,A13,A21,A22,A23,A31,A32,A33;
141    
# Line 124  void Finley_Util_InvertSmallMat(dim_t le Line 148  void Finley_Util_InvertSmallMat(dim_t le
148                 D=1./D;                 D=1./D;
149                 invA[q]=D;                 invA[q]=D;
150              } else {              } else {
151                 Finley_ErrorCode=ZERO_DIVISION_ERROR;                 Dudley_setError(ZERO_DIVISION_ERROR,"__FILE__: Non-regular matrix");
                sprintf(Finley_ErrorMsg,"Non-regular matrix");  
152                 return;                 return;
153              }              }
154           }           }
# Line 147  void Finley_Util_InvertSmallMat(dim_t le Line 170  void Finley_Util_InvertSmallMat(dim_t le
170                 invA[INDEX3(0,1,q,2,2)]=-A12*D;                 invA[INDEX3(0,1,q,2,2)]=-A12*D;
171                 invA[INDEX3(1,1,q,2,2)]= A11*D;                 invA[INDEX3(1,1,q,2,2)]= A11*D;
172              } else {              } else {
173                 Finley_ErrorCode=ZERO_DIVISION_ERROR;                 Dudley_setError(ZERO_DIVISION_ERROR,"__FILE__: Non-regular matrix");
                sprintf(Finley_ErrorMsg,"Non-regular matrix");  
174                 return;                 return;
175              }              }
176           }           }
# Line 180  void Finley_Util_InvertSmallMat(dim_t le Line 202  void Finley_Util_InvertSmallMat(dim_t le
202                 invA[INDEX3(1,2,q,3,3)]=(A13*A21-A11*A23)*D;                 invA[INDEX3(1,2,q,3,3)]=(A13*A21-A11*A23)*D;
203                 invA[INDEX3(2,2,q,3,3)]=(A11*A22-A12*A21)*D;                 invA[INDEX3(2,2,q,3,3)]=(A11*A22-A12*A21)*D;
204              } else {              } else {
205                 Finley_ErrorCode=ZERO_DIVISION_ERROR;                 Dudley_setError(ZERO_DIVISION_ERROR,"__FILE__: Non-regular matrix");
                sprintf(Finley_ErrorMsg,"Non-regular matrix");  
206                 return;                 return;
207              }              }
208           }           }
# Line 193  void Finley_Util_InvertSmallMat(dim_t le Line 214  void Finley_Util_InvertSmallMat(dim_t le
214    
215  /*    sets the derterminate of a set of dim x dim matrices A(:,:,1:len) with dim=1,2,3 */  /*    sets the derterminate of a set of dim x dim matrices A(:,:,1:len) with dim=1,2,3 */
216    
217  void Finley_Util_DetOfSmallMat(dim_t len,dim_t dim,double* A, double* det){  void Dudley_Util_DetOfSmallMat(dim_t len,dim_t dim,double* A, double* det){
218     dim_t q;     dim_t q;
219     register double A11,A12,A13,A21,A22,A23,A31,A32,A33;     register double A11,A12,A13,A21,A22,A23,A31,A32,A33;
220    
# Line 237  void Finley_Util_DetOfSmallMat(dim_t len Line 258  void Finley_Util_DetOfSmallMat(dim_t len
258  /*    returns the normalized vector Normal[dim,len] orthogonal to A(:,0,q) and A(:,1,q) in the case of dim=3  */  /*    returns the normalized vector Normal[dim,len] orthogonal to A(:,0,q) and A(:,1,q) in the case of dim=3  */
259  /*    or the vector A(:,0,q) in the case of dim=2                                             */  /*    or the vector A(:,0,q) in the case of dim=2                                             */
260    
261  void  Finley_NormalVector(dim_t len, dim_t dim, dim_t dim1, double* A,double* Normal) {  void  Dudley_NormalVector(dim_t len, dim_t dim, dim_t dim1, double* A,double* Normal) {
262     dim_t q;     dim_t q;
263     register double A11,A12,CO_A13,A21,A22,CO_A23,A31,A32,CO_A33,length,invlength;     register double A11,A12,CO_A13,A21,A22,CO_A23,A31,A32,CO_A33,length,invlength;
264    
# Line 251  void  Finley_NormalVector(dim_t len, dim Line 272  void  Finley_NormalVector(dim_t len, dim
272              A21=A[INDEX3(1,0,q,2,dim1)];              A21=A[INDEX3(1,0,q,2,dim1)];
273              length = sqrt(A11*A11+A21*A21);              length = sqrt(A11*A11+A21*A21);
274              if (! length>0) {              if (! length>0) {
275                 Finley_ErrorCode=ZERO_DIVISION_ERROR;                 Dudley_setError(ZERO_DIVISION_ERROR,"__FILE__: area equals zero.");
                sprintf(Finley_ErrorMsg,"area equals zero.");  
276                 return;                 return;
277              } else {              } else {
278                 invlength=1./length;                 invlength=1./length;
# Line 274  void  Finley_NormalVector(dim_t len, dim Line 294  void  Finley_NormalVector(dim_t len, dim
294              CO_A33=A11*A22-A21*A12;              CO_A33=A11*A22-A21*A12;
295              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);
296              if (! length>0) {              if (! length>0) {
297                 Finley_ErrorCode=ZERO_DIVISION_ERROR;                 Dudley_setError(ZERO_DIVISION_ERROR,"__FILE__: area equals zero.");
                sprintf(Finley_ErrorMsg,"area equals zero.");  
298                 return;                 return;
299              } else {              } else {
300                 invlength=1./length;                 invlength=1./length;
# Line 294  void  Finley_NormalVector(dim_t len, dim Line 313  void  Finley_NormalVector(dim_t len, dim
313  /*    return the length of the vector which is orthogonal to the vectors A(:,0,q) and A(:,1,q) in the case of dim=3 */  /*    return the length of the vector which is orthogonal to the vectors A(:,0,q) and A(:,1,q) in the case of dim=3 */
314  /*    or the vector A(:,0,q) in the case of dim=2                                                                   */  /*    or the vector A(:,0,q) in the case of dim=2                                                                   */
315    
316  void  Finley_LengthOfNormalVector(dim_t len, dim_t dim, dim_t dim1, double* A,double* length) {  void  Dudley_LengthOfNormalVector(dim_t len, dim_t dim, dim_t dim1, double* A,double* length) {
317     dim_t q;     dim_t q;
318     double A11,A12,CO_A13,A21,A22,CO_A23,A31,A32,CO_A33;     double A11,A12,CO_A13,A21,A22,CO_A23,A31,A32,CO_A33;
319    
# Line 332  void  Finley_LengthOfNormalVector(dim_t Line 351  void  Finley_LengthOfNormalVector(dim_t
351  /* there is no range checking! */  /* there is no range checking! */
352  /* at output Map[invMap[i]]=i for i=0:lenInvMap */  /* at output Map[invMap[i]]=i for i=0:lenInvMap */
353    
354  void Finley_Util_InvertMap(dim_t lenInvMap, index_t* invMap,dim_t lenMap, index_t* Map) {  void Dudley_Util_InvertMap(dim_t lenInvMap, index_t* invMap,dim_t lenMap, index_t* Map) {
355     dim_t i;     dim_t i;
356     for (i=0;i<lenInvMap;i++) invMap[i]=0;     for (i=0;i<lenInvMap;i++) invMap[i]=0;
357     for (i=0;i<lenMap;i++) {     for (i=0;i<lenMap;i++) {
# Line 340  void Finley_Util_InvertMap(dim_t lenInvM Line 359  void Finley_Util_InvertMap(dim_t lenInvM
359     }     }
360  }  }
361    
362  /* orders a Finley_Util_ValueAndIndex array by value */  /* orders a Dudley_Util_ValueAndIndex array by value */
363  /* it is assumed that n is large */  /* it is assumed that n is large */
364    
365  int Finley_Util_ValueAndIndex_compar(const void *arg1 , const void *arg2 ) {  int Dudley_Util_ValueAndIndex_compar(const void *arg1 , const void *arg2 ) {
366     Finley_Util_ValueAndIndex *e1,*e2;     Dudley_Util_ValueAndIndex *e1,*e2;
367     e1=(Finley_Util_ValueAndIndex*) arg1;     e1=(Dudley_Util_ValueAndIndex*) arg1;
368     e2=(Finley_Util_ValueAndIndex*) arg2;     e2=(Dudley_Util_ValueAndIndex*) arg2;
369     if (e1->value < e2->value) return -1;     if (e1->value < e2->value) return -1;
370     if (e1->value > e2->value) return  1;     if (e1->value > e2->value) return  1;
371       if (e1->index < e2->index) return -1;
372       if (e1->index > e2->index) return  1;
373     return 0;     return 0;
374  }  }
375  void Finley_Util_sortValueAndIndex(dim_t n,Finley_Util_ValueAndIndex* array) {  
376    void Dudley_Util_sortValueAndIndex(dim_t n,Dudley_Util_ValueAndIndex* array) {
377       /* OMP : needs parallelization !*/       /* OMP : needs parallelization !*/
378       qsort(array,n,sizeof(Finley_Util_ValueAndIndex),Finley_Util_ValueAndIndex_compar);       qsort(array,n,sizeof(Dudley_Util_ValueAndIndex),Dudley_Util_ValueAndIndex_compar);
379  }  }
380    
381    
# Line 361  void Finley_Util_sortValueAndIndex(dim_t Line 383  void Finley_Util_sortValueAndIndex(dim_t
383    
384  /* calculates the minimum value from a dim X N integer array */  /* calculates the minimum value from a dim X N integer array */
385    
386  index_t Finley_Util_getMinInt(dim_t dim,dim_t N,index_t* values) {  index_t Dudley_Util_getMinInt(dim_t dim,dim_t N,index_t* values) {
387     dim_t i,j;     dim_t i,j;
388     index_t out,out_local;     index_t out,out_local;
389     out=INDEX_T_MAX;     out=INDEX_T_MAX;
# Line 383  index_t Finley_Util_getMinInt(dim_t dim, Line 405  index_t Finley_Util_getMinInt(dim_t dim,
405                                                                                                                                                                                                                                                                                                        
406  /* calculates the maximum value from a dim X N integer array */  /* calculates the maximum value from a dim X N integer array */
407    
408  index_t Finley_Util_getMaxInt(dim_t dim,dim_t N,index_t* values) {  index_t Dudley_Util_getMaxInt(dim_t dim,dim_t N,index_t* values) {
409       dim_t i,j;
410       index_t out,out_local;
411       out=-INDEX_T_MAX;
412       if (values!=NULL && dim*N>0 ) {
413         out=values[0];
414         #pragma omp parallel private(out_local)
415         {
416             out_local=out;
417             #pragma omp for private(i,j) schedule(static)
418             for (j=0;j<N;j++) {
419                 for (i=0;i<dim;i++)
420    {
421    //printf("%d,%d,%d[%d] %d\n",i,j,dim,INDEX2(i,j,dim),  values[INDEX2(i,j,dim)]);
422    out_local=MAX(out_local,values[INDEX2(i,j,dim)]);
423    
424    }
425             }
426             #pragma omp critical
427             out=MAX(out_local,out);
428          }
429       }
430       return out;
431    }
432    /**************************************************************/
433    
434    /* calculates the minimum value from a dim X N integer array */
435    
436    index_t Dudley_Util_getFlaggedMinInt(dim_t dim,dim_t N,index_t* values, index_t ignore) {
437       dim_t i,j;
438       index_t out,out_local;
439       out=INDEX_T_MAX;
440       if (values!=NULL && dim*N>0 ) {
441         out=values[0];
442         #pragma omp parallel private(out_local)
443         {
444             out_local=out;
445             #pragma omp for private(i,j) schedule(static)
446             for (j=0;j<N;j++) {
447               for (i=0;i<dim;i++) if (values[INDEX2(i,j,dim)]!=ignore) out_local=MIN(out_local,values[INDEX2(i,j,dim)]);
448             }
449             #pragma omp critical
450             out=MIN(out_local,out);
451         }
452       }
453       return out;
454    }
455                                                                                                                                                      
456    /* calculates the maximum value from a dim X N integer array */
457    
458    index_t Dudley_Util_getFlaggedMaxInt(dim_t dim,dim_t N,index_t* values, index_t ignore) {
459     dim_t i,j;     dim_t i,j;
460     index_t out,out_local;     index_t out,out_local;
461     out=-INDEX_T_MAX;     out=-INDEX_T_MAX;
# Line 394  index_t Finley_Util_getMaxInt(dim_t dim, Line 466  index_t Finley_Util_getMaxInt(dim_t dim,
466           out_local=out;           out_local=out;
467           #pragma omp for private(i,j) schedule(static)           #pragma omp for private(i,j) schedule(static)
468           for (j=0;j<N;j++) {           for (j=0;j<N;j++) {
469               for (i=0;i<dim;i++) out_local=MAX(out_local,values[INDEX2(i,j,dim)]);               for (i=0;i<dim;i++) if (values[INDEX2(i,j,dim)]!=ignore) out_local=MAX(out_local,values[INDEX2(i,j,dim)]);
470           }           }
471           #pragma omp critical           #pragma omp critical
472           out=MAX(out_local,out);           out=MAX(out_local,out);
# Line 405  index_t Finley_Util_getMaxInt(dim_t dim, Line 477  index_t Finley_Util_getMaxInt(dim_t dim,
477    
478  /* 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. */
479    
480  dim_t Finley_Util_packMask(dim_t N,index_t* mask,index_t* index) {  dim_t Dudley_Util_packMask(dim_t N,index_t* mask,index_t* index) {
481        dim_t out,k;        dim_t out,k;
482        out=0;        out=0;
483        /*OMP */        /*OMP */
# Line 419  dim_t Finley_Util_packMask(dim_t N,index Line 491  dim_t Finley_Util_packMask(dim_t N,index
491  }  }
492    
493  /* returns true if array contains value */  /* returns true if array contains value */
494  bool_t Finley_Util_isAny(dim_t N,index_t* array,index_t value) {  bool_t Dudley_Util_isAny(dim_t N,index_t* array,index_t value) {
495     bool_t out=FALSE;     bool_t out=FALSE;
496     dim_t i;     dim_t i;
497     #pragma omp parallel for private(i) schedule(static) reduction(||:out)     #pragma omp parallel for private(i) schedule(static) reduction(||:out)
# Line 427  bool_t Finley_Util_isAny(dim_t N,index_t Line 499  bool_t Finley_Util_isAny(dim_t N,index_t
499     return out;     return out;
500  }  }
501  /* calculates the cummultative sum in array and returns the total sum */  /* calculates the cummultative sum in array and returns the total sum */
502  index_t Finley_Util_cumsum(dim_t N,index_t* array) {  index_t Dudley_Util_cumsum(dim_t N,index_t* array) {
503     index_t out=0,tmp;     index_t out=0,tmp;
504     dim_t i;     dim_t i;
505     #ifdef _OPENMP     #ifdef _OPENMP
506        index_t partial_sums[omp_get_max_threads()],sum;        index_t *partial_sums=NULL, sum;
507          partial_sums=TMPMEMALLOC(omp_get_max_threads(),index_t);
508        #pragma omp parallel private(sum,i,tmp)        #pragma omp parallel private(sum,i,tmp)
509        {        {
510          sum=0;          sum=0;
# Line 457  index_t Finley_Util_cumsum(dim_t N,index Line 530  index_t Finley_Util_cumsum(dim_t N,index
530            array[i]=tmp;            array[i]=tmp;
531          }          }
532        }        }
533          TMPMEMFREE(partial_sums);
534     #else     #else
535        for (i=0;i<N;++i) {        for (i=0;i<N;++i) {
536           tmp=out;           tmp=out;
# Line 466  index_t Finley_Util_cumsum(dim_t N,index Line 540  index_t Finley_Util_cumsum(dim_t N,index
540     #endif     #endif
541     return out;     return out;
542  }  }
543    void Dudley_Util_setValuesInUse(const index_t *values, const dim_t numValues, dim_t *numValuesInUse, index_t **valuesInUse, Paso_MPIInfo* mpiinfo)
544    {
545       dim_t i;
546       index_t lastFoundValue=INDEX_T_MIN, minFoundValue, local_minFoundValue, *newValuesInUse=NULL;
547       register index_t itmp;
548       bool_t allFound=FALSE;
549       dim_t nv=0;
550    
551       while (! allFound) {
552           /*
553            *  find smallest value bigger than lastFoundValue
554            */
555            minFoundValue=INDEX_T_MAX;
556            #pragma omp parallel private(local_minFoundValue)
557            {
558                local_minFoundValue=minFoundValue;
559                #pragma omp for private(i,itmp) schedule(static)
560                for (i=0;i< numValues;i++) {
561                   itmp=values[i];
562                   if ((itmp>lastFoundValue) && (itmp<local_minFoundValue)) local_minFoundValue=itmp;
563                }
564                #pragma omp critical
565                {
566                   if (local_minFoundValue<minFoundValue) minFoundValue=local_minFoundValue;
567                }
568    
569             }
570             #ifdef PASO_MPI
571             local_minFoundValue=minFoundValue;
572             MPI_Allreduce(&local_minFoundValue,&minFoundValue, 1, MPI_INT, MPI_MIN, mpiinfo->comm );
573             #endif
574             /* if we found a new tag we need to add this too the valuesInUseList */
575    
576             if (minFoundValue < INDEX_T_MAX) {
577                 newValuesInUse=MEMALLOC(nv+1,index_t);
578                 if (*valuesInUse!=NULL) {
579                     memcpy(newValuesInUse,*valuesInUse,sizeof(index_t)*nv);
580                     MEMFREE(*valuesInUse);
581                 }
582                 newValuesInUse[nv]=minFoundValue;
583                 *valuesInUse=newValuesInUse;
584                 newValuesInUse=NULL;
585                 nv++;
586                 lastFoundValue=minFoundValue;
587             } else {
588                 allFound=TRUE;
589             }
590       }
591       *numValuesInUse=nv;
592    }
593    
594    
595  void Finley_copyDouble(dim_t n,double* source, double* target) {  #ifdef PASO_MPI
596    dim_t i;  void Dudley_printDoubleArray( FILE *fid, dim_t n, double *array, char *name  )
597    for (i=0;i<n;i++) target[i]=source[i];  {
598  }    index_t i;
599      
600  /*    if( name )
601   * $Log$      fprintf( fid, "%s [ ", name );
602   * Revision 1.8  2005/08/12 01:45:43  jgs    else
603   * erge of development branch dev-02 back to main trunk on 2005-08-12      fprintf( fid, "[ " );  
604   *    for( i=0; i<(n<60 ? n : 60); i++ )
605   * Revision 1.7.2.1  2005/08/04 22:41:11  gross      fprintf( fid, "%g ", array[i] );
606   * some extra routines for finley that might speed-up RHS assembling in some cases (not actived right now)    if( n>=30 )
607   *      fprintf( fid, "... " );
608   * Revision 1.7  2005/07/08 04:07:59  jgs    fprintf( fid, "]\n" );
609   * Merge of development branch back to main trunk on 2005-07-08  }
610   *  void Dudley_printIntArray( FILE *fid, dim_t n, int *array, char *name  )
611   * Revision 1.1.1.1.2.4  2005/06/29 02:34:57  gross  {
612   * some changes towards 64 integers in finley    index_t i;
613   *    
614   * Revision 1.1.1.1.2.3  2005/03/02 23:35:06  gross    if( name )
615   * reimplementation of the ILU in Finley. block size>1 still needs some testing      fprintf( fid, "%s [ ", name );
616   *    else
617   * Revision 1.1.1.1.2.2  2005/02/18 02:27:31  gross      fprintf( fid, "[ " );  
618   * two function that will be used for a reimplementation of the ILU preconditioner    for( i=0; i<(n<60 ? n : 60); i++ )
619   *      fprintf( fid, "%d ", array[i] );
620   * Revision 1.1.1.1.2.1  2004/11/12 06:58:19  gross    if( n>=30 )
621   * 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, "... " );
622   *    fprintf( fid, "]\n" );
623   * Revision 1.1.1.1  2004/10/26 06:53:57  jgs  }
624   * initial import of project esys2  void Dudley_printMaskArray( FILE *fid, dim_t n, int *array, char *name  )
625   *  {
626   * Revision 1.3  2004/08/26 12:03:52  gross    index_t i;
627   * Some other bug in Finley_Assemble_gradient fixed.    
628   *    if( name )
629   * Revision 1.2  2004/07/02 04:21:13  gross      fprintf( fid, "%s [ ", name );
630   * Finley C code has been included    else
631   *      fprintf( fid, "[ " );  
632   * Revision 1.1.1.1  2004/06/24 04:00:40  johng    for( i=0; i<(n<60 ? n : 60); i++ )
633   * Initial version of eys using boost-python.      if( array[i]!=-1 )
634   *        fprintf( fid, "%3d ", array[i] );
635   *      else
636   */        fprintf( fid, "  * " );
637      if( n>=30 )
638        fprintf( fid, "... " );
639      fprintf( fid, "]\n" );
640    }
641    #endif

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

  ViewVC Help
Powered by ViewVC 1.1.26