/[escript]/branches/trilinos_from_5897/dudley/src/Util.cpp
ViewVC logotype

Diff of /branches/trilinos_from_5897/dudley/src/Util.cpp

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

revision 6078 by caltinay, Wed Mar 2 04:13:26 2016 UTC revision 6079 by caltinay, Mon Mar 21 12:22:38 2016 UTC
# Line 18  Line 18 
18    
19  #include <escript/index.h>  #include <escript/index.h>
20    
21  #include <limits.h>  #include <algorithm> // std::sort
 #include <cstring>  /* for memcpy*/  
22    
23  namespace dudley {  namespace dudley {
24    namespace util {
25    
26  /* returns true if any of the values in the short array values is not equal to zero */  /// comparison function for sortValueAndIndex
27  bool Dudley_Util_anyNonZeroDouble(dim_t N, double *values)  bool ValueAndIndexCompare(const std::pair<int,int> &i, const std::pair<int, int> &j)
28  {  {
29      dim_t q;      // to ensure we have a strict ordering as required by std
30      for (q = 0; q < N; ++q)      if (i.first == j.first)
31          if (std::abs(values[q]) > 0)          return i.second < j.second;
32              return true;      return i.first < j.first;
     return false;  
33  }  }
34    
35    void sortValueAndIndex(ValueAndIndexList& array)
 /*   gathers double values out from in by index: */  
 /*        out(1:numData,1:len)=in(1:numData,index(1:len)) */  
 void Dudley_Util_Gather_double(dim_t len, index_t * index, dim_t numData, double *in, double *out)  
36  {  {
37      dim_t s, i;      std::sort(array.begin(), array.end(), ValueAndIndexCompare);
     for (s = 0; s < len; s++) {  
         for (i = 0; i < numData; i++) {  
             out[INDEX2(i, s, numData)] = in[INDEX2(i, index[s], numData)];  
         }  
     }  
38  }  }
39    
40  /*   gathers maybelong values out from in by index: */  void gather(int len, const index_t* index, int numData, const double* in,
41  /*        out(1:numData,1:len)=in(1:numData,index(1:len)) */              double* out)
42  void Dudley_Util_Gather_int(dim_t len, index_t * index, dim_t numData, index_t * in, index_t * out)  {
43  {      for (int s = 0; s < len; s++) {
44      dim_t s, i;          for (int i = 0; i < numData; i++) {
     for (s = 0; s < len; s++) {  
         for (i = 0; i < numData; i++) {  
45              out[INDEX2(i, s, numData)] = in[INDEX2(i, index[s], numData)];              out[INDEX2(i, s, numData)] = in[INDEX2(i, index[s], numData)];
46          }          }
47      }      }
48  }  }
49    
50  /*   adds a vector in into out using and index. */  void addScatter(int len, const index_t* index, int numData,
51  /*        out(1:numData,index[p])+=in(1:numData,p) where p = {k=1...len , index[k]<upperBound}*/                  const double* in, double* out, index_t upperBound)
52  void Dudley_Util_AddScatter(dim_t len, const index_t* index, dim_t numData, const double *in, double *out, index_t upperBound)  {
53  {      for (int s = 0; s < len; s++) {
54      dim_t i, s;          for (int i = 0; i < numData; i++) {
     for (s = 0; s < len; s++) {  
         for (i = 0; i < numData; i++) {  
55              if (index[s] < upperBound) {              if (index[s] < upperBound) {
56                  out[INDEX2(i, index[s], numData)] += in[INDEX2(i, s, numData)];                  out[INDEX2(i, index[s], numData)] += in[INDEX2(i, s, numData)];
57              }              }
# Line 72  void Dudley_Util_AddScatter(dim_t len, c Line 59  void Dudley_Util_AddScatter(dim_t len, c
59      }      }
60  }  }
61    
62  /*    multiplies two matrices */  void smallMatMult(int A1, int A2, double* A, int B2, const double* B,
63  /*          A(1:A1,1:A2)=B(1:A1,1:B2)*C(1:B2,1:A2) */                    const double* C)
 void Dudley_Util_SmallMatMult(dim_t A1, dim_t A2, double *A, dim_t B2, const double *B, const double *C)  
 {  
     dim_t i, j, s;  
     double rtmp;  
     for (i = 0; i < A1; i++) {  
         for (j = 0; j < A2; j++) {  
             rtmp = 0;  
             for (s = 0; s < B2; s++) {  
                 rtmp += B[INDEX2(i, s, A1)] * C[INDEX2(s, j, B2)];  
             }  
             A[INDEX2(i, j, A1)] = rtmp;  
         }  
     }  
 }  
   
 /*    multiplies two sets of matrices: */  
 /*        A(1:A1,1:A2,i)=B(1:A1,1:B2,i)*C(1:B2,1:A2,i) i=1,len */  
 void Dudley_Util_SmallMatSetMult(dim_t len, dim_t A1, dim_t A2, double *A, dim_t B2, const double *B, const double *C)  
 {  
     dim_t q, i, j, s;  
     double rtmp;  
     for (q = 0; q < len; q++) {  
         for (i = 0; i < A1; i++) {  
             for (j = 0; j < A2; j++) {  
                 rtmp = 0;  
                 for (s = 0; s < B2; s++)  
                     rtmp += B[INDEX3(i, s, q, A1, B2)] * C[INDEX3(s, j, q, B2, A2)];  
                 A[INDEX3(i, j, q, A1, A2)] = rtmp;  
             }  
         }  
     }  
 }  
   
 /*    multiplies a set of matrices with a single matrix: */  
 /*        A(1:A1,1:A2,i)=B(1:A1,1:B2,i)*C(1:B2,1:A2) i=1,len */  
 void Dudley_Util_SmallMatSetMult1(dim_t len, dim_t A1, dim_t A2, double *A, dim_t B2, const double *B, const double *C)  
 {  
     dim_t q, i, j, s;  
     double rtmp;  
     for (q = 0; q < len; q++) {  
         for (i = 0; i < A1; i++) {  
             for (j = 0; j < A2; j++) {  
                 rtmp = 0;  
                 for (s = 0; s < B2; s++)  
                     rtmp += B[INDEX3(i, s, q, A1, B2)] * C[INDEX2(s, j, B2)];  
                 A[INDEX3(i, j, q, A1, A2)] = rtmp;  
             }  
         }  
     }  
 }  
   
 /*    inverts the set of dim x dim matrices A(:,:,1:len) with dim=1,2,3 */  
 /*    the determinant is returned. */  
 void Dudley_Util_InvertSmallMat(dim_t len, dim_t dim, double *A, double *invA, double *det)  
64  {  {
65      dim_t q;      for (int i = 0; i < A1; i++) {
66      double D, A11, A12, A13, A21, A22, A23, A31, A32, A33;          for (int j = 0; j < A2; j++) {
67                double sum = 0.;
68      switch (dim) {              for (int s = 0; s < B2; s++)
69      case 1:                  sum += B[INDEX2(i,s,A1)] * C[INDEX2(s,j,B2)];
70          for (q = 0; q < len; q++) {              A[INDEX2(i,j,A1)] = sum;
             D = A[q];  
             if (std::abs(D) > 0) {  
                 det[q] = D;  
                 D = 1. / D;  
                 invA[q] = D;  
             } else {  
                 throw DudleyException("InvertSmallMat: Non-regular matrix");  
             }  
71          }          }
         break;  
   
     case 2:  
         for (q = 0; q < len; q++) {  
             A11 = A[INDEX3(0, 0, q, 2, 2)];  
             A12 = A[INDEX3(0, 1, q, 2, 2)];  
             A21 = A[INDEX3(1, 0, q, 2, 2)];  
             A22 = A[INDEX3(1, 1, q, 2, 2)];  
   
             D = A11 * A22 - A12 * A21;  
             if (std::abs(D) > 0) {  
                 det[q] = D;  
                 D = 1. / D;  
                 invA[INDEX3(0, 0, q, 2, 2)] = A22 * D;  
                 invA[INDEX3(1, 0, q, 2, 2)] = -A21 * D;  
                 invA[INDEX3(0, 1, q, 2, 2)] = -A12 * D;  
                 invA[INDEX3(1, 1, q, 2, 2)] = A11 * D;  
             } else {  
                 throw DudleyException("InvertSmallMat: Non-regular matrix");  
             }  
         }  
         break;  
   
     case 3:  
         for (q = 0; q < len; q++) {  
             A11 = A[INDEX3(0, 0, q, 3, 3)];  
             A21 = A[INDEX3(1, 0, q, 3, 3)];  
             A31 = A[INDEX3(2, 0, q, 3, 3)];  
             A12 = A[INDEX3(0, 1, q, 3, 3)];  
             A22 = A[INDEX3(1, 1, q, 3, 3)];  
             A32 = A[INDEX3(2, 1, q, 3, 3)];  
             A13 = A[INDEX3(0, 2, q, 3, 3)];  
             A23 = A[INDEX3(1, 2, q, 3, 3)];  
             A33 = A[INDEX3(2, 2, q, 3, 3)];  
   
             D = A11 * (A22 * A33 - A23 * A32) + A12 * (A31 * A23 - A21 * A33) + A13 * (A21 * A32 - A31 * A22);  
             if (std::abs(D) > 0) {  
                 det[q] = D;  
                 D = 1. / D;  
                 invA[INDEX3(0, 0, q, 3, 3)] = (A22 * A33 - A23 * A32) * D;  
                 invA[INDEX3(1, 0, q, 3, 3)] = (A31 * A23 - A21 * A33) * D;  
                 invA[INDEX3(2, 0, q, 3, 3)] = (A21 * A32 - A31 * A22) * D;  
                 invA[INDEX3(0, 1, q, 3, 3)] = (A13 * A32 - A12 * A33) * D;  
                 invA[INDEX3(1, 1, q, 3, 3)] = (A11 * A33 - A31 * A13) * D;  
                 invA[INDEX3(2, 1, q, 3, 3)] = (A12 * A31 - A11 * A32) * D;  
                 invA[INDEX3(0, 2, q, 3, 3)] = (A12 * A23 - A13 * A22) * D;  
                 invA[INDEX3(1, 2, q, 3, 3)] = (A13 * A21 - A11 * A23) * D;  
                 invA[INDEX3(2, 2, q, 3, 3)] = (A11 * A22 - A12 * A21) * D;  
             } else {  
                 throw DudleyException("InvertSmallMat: Non-regular matrix");  
             }  
         }  
         break;  
72      }      }
73  }  }
74    
75  /*    sets the derterminat of a set of dim x dim matrices A(:,:,1:len) with dim=1,2,3 */  void smallMatSetMult1(int len, int A1, int A2, double* A, int B2,
76  void Dudley_Util_DetOfSmallMat(dim_t len, dim_t dim, double *A, double *det)                        const double* B, const double* C)
77  {  {
78      dim_t q;      for (int q = 0; q < len; q++) {
79      double A11, A12, A13, A21, A22, A23, A31, A32, A33;          for (int i = 0; i < A1; i++) {
80                for (int j = 0; j < A2; j++) {
81      switch (dim) {                  double sum = 0.;
82      case 1:                  for (int s = 0; s < B2; s++)
83          for (q = 0; q < len; q++) {                      sum += B[INDEX3(i,s,q,A1,B2)] * C[INDEX2(s,j,B2)];
84              det[q] = A[q];                  A[INDEX3(i,j,q,A1,A2)] = sum;
         }  
         break;  
   
     case 2:  
         for (q = 0; q < len; q++) {  
             A11 = A[INDEX3(0, 0, q, 2, 2)];  
             A12 = A[INDEX3(0, 1, q, 2, 2)];  
             A21 = A[INDEX3(1, 0, q, 2, 2)];  
             A22 = A[INDEX3(1, 1, q, 2, 2)];  
   
             det[q] = A11 * A22 - A12 * A21;  
         }  
         break;  
   
     case 3:  
         for (q = 0; q < len; q++) {  
             A11 = A[INDEX3(0, 0, q, 3, 3)];  
             A21 = A[INDEX3(1, 0, q, 3, 3)];  
             A31 = A[INDEX3(2, 0, q, 3, 3)];  
             A12 = A[INDEX3(0, 1, q, 3, 3)];  
             A22 = A[INDEX3(1, 1, q, 3, 3)];  
             A32 = A[INDEX3(2, 1, q, 3, 3)];  
             A13 = A[INDEX3(0, 2, q, 3, 3)];  
             A23 = A[INDEX3(1, 2, q, 3, 3)];  
             A33 = A[INDEX3(2, 2, q, 3, 3)];  
   
             det[q] = A11 * (A22 * A33 - A23 * A32) + A12 * (A31 * A23 - A21 * A33) + A13 * (A21 * A32 - A31 * A22);  
         }  
         break;  
   
     }  
 }  
   
 /*    returns the normalized vector Normal[dim,len] orthogonal to A(:,0,q) and A(:,1,q) in the case of dim=3  */  
 /*    or the vector A(:,0,q) in the case of dim=2 */  
 void Dudley_NormalVector(dim_t len, dim_t dim, dim_t dim1, double *A, double *Normal)  
 {  
     dim_t q;  
     double A11, A12, CO_A13, A21, A22, CO_A23, A31, A32, CO_A33, length, invlength;  
   
     switch (dim) {  
     case 1:  
         for (q = 0; q < len; q++)  
             Normal[q] = 1;  
         break;  
     case 2:  
         for (q = 0; q < len; q++) {  
             A11 = A[INDEX3(0, 0, q, 2, dim1)];  
             A21 = A[INDEX3(1, 0, q, 2, dim1)];  
             length = sqrt(A11 * A11 + A21 * A21);  
             if (length <= 0) {  
                 throw DudleyException("NormalVector: area equals zero.");  
             } else {  
                 invlength = 1. / length;  
                 Normal[INDEX2(0, q, 2)] = A21 * invlength;  
                 Normal[INDEX2(1, q, 2)] = -A11 * invlength;  
85              }              }
86          }          }
         break;  
     case 3:  
         for (q = 0; q < len; q++)  
         {  
             A11 = A[INDEX3(0, 0, q, 3, dim1)];  
             A21 = A[INDEX3(1, 0, q, 3, dim1)];  
             A31 = A[INDEX3(2, 0, q, 3, dim1)];  
             A12 = A[INDEX3(0, 1, q, 3, dim1)];  
             A22 = A[INDEX3(1, 1, q, 3, dim1)];  
             A32 = A[INDEX3(2, 1, q, 3, dim1)];  
             CO_A13 = A21 * A32 - A31 * A22;  
             CO_A23 = A31 * A12 - A11 * A32;  
             CO_A33 = A11 * A22 - A21 * A12;  
             length = sqrt(CO_A13 * CO_A13 + CO_A23 * CO_A23 + CO_A33 * CO_A33);  
             if (length <= 0) {  
                 throw DudleyException("NormalVector: area equals zero.");  
             } else {  
                 invlength = 1. / length;  
                 Normal[INDEX2(0, q, 3)] = CO_A13 * invlength;  
                 Normal[INDEX2(1, q, 3)] = CO_A23 * invlength;  
                 Normal[INDEX2(2, q, 3)] = CO_A33 * invlength;  
             }  
   
         }  
         break;  
87      }      }
88  }  }
89    
90  /*    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 */  void normalVector(int len, int dim, int dim1, const double* A, double* Normal)
 /*    or the vector A(:,0,q) in the case of dim=2                                                                   */  
   
 void Dudley_LengthOfNormalVector(dim_t len, dim_t dim, dim_t dim1, double *A, double *length)  
91  {  {
92      dim_t q;      int q;
     double A11, A12, CO_A13, A21, A22, CO_A23, A31, A32, CO_A33;  
93    
94      switch (dim) {      switch (dim) {
95      case 1:          case 1:
96          for (q = 0; q < len; q++)              for (q = 0; q < len; q++)
97              length[q] = 1;                  Normal[q] = 1.;
98          break;              break;
99      case 2:          case 2:
100          for (q = 0; q < len; q++) {              for (q = 0; q < len; q++) {
101              A11 = A[INDEX3(0, 0, q, 2, dim1)];                  const double A11 = A[INDEX3(0,0,q,2,dim1)];
102              A21 = A[INDEX3(1, 0, q, 2, dim1)];                  const double A21 = A[INDEX3(1,0,q,2,dim1)];
103              length[q] = sqrt(A11 * A11 + A21 * A21);                  const double length = sqrt(A11*A11+A21*A21);
104          }                  if (length <= 0) {
105          break;                      throw DudleyException("normalVector: area equals zero.");
106      case 3:                  } else {
107          for (q = 0; q < len; q++) {                      const double invlength = 1./length;
108              A11 = A[INDEX3(0, 0, q, 3, dim1)];                      Normal[INDEX2(0,q,2)] =  A21*invlength;
109              A21 = A[INDEX3(1, 0, q, 3, dim1)];                      Normal[INDEX2(1,q,2)] = -A11*invlength;
110              A31 = A[INDEX3(2, 0, q, 3, dim1)];                  }
             A12 = A[INDEX3(0, 1, q, 3, dim1)];  
             A22 = A[INDEX3(1, 1, q, 3, dim1)];  
             A32 = A[INDEX3(2, 1, q, 3, dim1)];  
             CO_A13 = A21 * A32 - A31 * A22;  
             CO_A23 = A31 * A12 - A11 * A32;  
             CO_A33 = A11 * A22 - A21 * A12;  
             length[q] = sqrt(CO_A13 * CO_A13 + CO_A23 * CO_A23 + CO_A33 * CO_A33);  
         }  
         break;  
     }  
 }  
   
 /* inverts the map map of length len */  
 /* there is no range checking! */  
 /* at output Map[invMap[i]]=i for i=0:lenInvMap */  
 void Dudley_Util_InvertMap(dim_t lenInvMap, index_t * invMap, dim_t lenMap, index_t * Map)  
 {  
     dim_t i;  
     for (i = 0; i < lenInvMap; i++)  
         invMap[i] = 0;  
     for (i = 0; i < lenMap; i++)  
         if (Map[i] >= 0)  
             invMap[Map[i]] = i;  
 }  
   
 /* orders a Dudley_Util_ValueAndIndex array by value */  
 /* it is assumed that n is large */  
 int Dudley_Util_ValueAndIndex_compar(const void *arg1, const void *arg2)  
 {  
     Dudley_Util_ValueAndIndex *e1, *e2;  
     e1 = (Dudley_Util_ValueAndIndex *) arg1;  
     e2 = (Dudley_Util_ValueAndIndex *) arg2;  
     if (e1->value < e2->value)  
         return -1;  
     if (e1->value > e2->value)  
         return 1;  
     if (e1->index < e2->index)  
         return -1;  
     if (e1->index > e2->index)  
         return 1;  
     return 0;  
 }  
   
 void Dudley_Util_sortValueAndIndex(dim_t n, Dudley_Util_ValueAndIndex * array)  
 {  
     /* OMP : needs parallelization ! */  
     qsort(array, n, sizeof(Dudley_Util_ValueAndIndex), Dudley_Util_ValueAndIndex_compar);  
 }  
   
 /* calculates the minimum value from a dim X N integer array */  
 index_t Dudley_Util_getMinInt(dim_t dim, dim_t N, index_t * values)  
 {  
     dim_t i, j;  
     index_t out, out_local;  
     out = escript::DataTypes::index_t_max();  
     if (values != NULL && dim * N > 0) {  
         out = values[0];  
 #pragma omp parallel private(out_local)  
         {  
             out_local = out;  
 #pragma omp for private(i,j) schedule(static)  
             for (j = 0; j < N; j++) {  
                 for (i = 0; i < dim; i++)  
                     out_local = std::min(out_local, values[INDEX2(i, j, dim)]);  
111              }              }
112  #pragma omp critical              break;
113              out = std::min(out_local, out);          case 3:
114          }              for (q = 0; q < len; q++) {
115      }                  const double A11 = A[INDEX3(0,0,q,3,dim1)];
116      return out;                  const double A21 = A[INDEX3(1,0,q,3,dim1)];
117  }                  const double A31 = A[INDEX3(2,0,q,3,dim1)];
118                    const double A12 = A[INDEX3(0,1,q,3,dim1)];
119  /* calculates the maximum value from a dim X N integer array */                  const double A22 = A[INDEX3(1,1,q,3,dim1)];
120  index_t Dudley_Util_getMaxInt(dim_t dim, dim_t N, index_t* values)                  const double A32 = A[INDEX3(2,1,q,3,dim1)];
121  {                  const double CO_A13 = A21*A32-A31*A22;
122      dim_t i, j;                  const double CO_A23 = A31*A12-A11*A32;
123      index_t out, out_local;                  const double CO_A33 = A11*A22-A21*A12;
124      out = -escript::DataTypes::index_t_max();                  const double length = sqrt(CO_A13*CO_A13 + CO_A23*CO_A23
125      if (values != NULL && dim * N > 0) {                                             + CO_A33*CO_A33);
126          out = values[0];                  if (length <= 0) {
127  #pragma omp parallel private(out_local)                      throw DudleyException("normalVector: area equals zero.");
128          {                  } else {
129              out_local = out;                      const double invlength = 1./length;
130  #pragma omp for private(i,j) schedule(static)                      Normal[INDEX2(0,q,3)] = CO_A13*invlength;
131              for (j = 0; j < N; j++) {                      Normal[INDEX2(1,q,3)] = CO_A23*invlength;
132                  for (i = 0; i < dim; i++) {                      Normal[INDEX2(2,q,3)] = CO_A33*invlength;
                     out_local = std::max(out_local, values[INDEX2(i, j, dim)]);  
133                  }                  }
134              }              }
135  #pragma omp critical              break;
             out = std::max(out_local, out);  
         }  
136      }      }
     return out;  
137  }  }
138    
139    IndexPair getMinMaxInt(int dim, dim_t N, const index_t* values)
 /* calculates the minimum value from a dim X N integer array */  
 index_t Dudley_Util_getFlaggedMinInt(dim_t dim, dim_t N, index_t* values, index_t ignore)  
140  {  {
141      dim_t i, j;      index_t vmin = escript::DataTypes::index_t_max();
142      index_t out, out_local;      index_t vmax = escript::DataTypes::index_t_min();
143      out = escript::DataTypes::index_t_max();      if (values && dim*N > 0) {
144      if (values != NULL && dim * N > 0) {          vmin = vmax = values[0];
145          out = values[0];  #pragma omp parallel
 #pragma omp parallel private(out_local)  
146          {          {
147              out_local = out;              index_t vmin_local = vmin;
148  #pragma omp for private(i,j) schedule(static)              index_t vmax_local = vmax;
149              for (j = 0; j < N; j++) {  #pragma omp for
150                  for (i = 0; i < dim; i++)              for (index_t j = 0; j < N; j++) {
151                      if (values[INDEX2(i, j, dim)] != ignore)                  for (int i = 0; i < dim; i++) {
152                          out_local = std::min(out_local, values[INDEX2(i, j, dim)]);                      vmin_local = std::min(vmin_local, values[INDEX2(i,j,dim)]);
153                        vmax_local = std::max(vmax_local, values[INDEX2(i,j,dim)]);
154                    }
155              }              }
156  #pragma omp critical  #pragma omp critical
157              out = std::min(out_local, out);              {
158                    vmin = std::min(vmin_local, vmin);
159                    vmax = std::max(vmax_local, vmax);
160                }
161          }          }
162      }      }
163      return out;      return IndexPair(vmin,vmax);
164  }  }
165    
166  /* calculates the maximum value from a dim X N integer array */  IndexPair getFlaggedMinMaxInt(dim_t N, const index_t* values, index_t ignore)
 index_t Dudley_Util_getFlaggedMaxInt(dim_t dim, dim_t N, index_t* values, index_t ignore)  
167  {  {
168      dim_t i, j;      index_t vmin = escript::DataTypes::index_t_max();
169      index_t out, out_local;      index_t vmax = escript::DataTypes::index_t_min();
170      out = -escript::DataTypes::index_t_max();      if (values && N > 0) {
171      if (values != NULL && dim * N > 0) {          vmin = vmax = values[0];
172          out = values[0];  #pragma omp parallel
 #pragma omp parallel private(out_local)  
173          {          {
174              out_local = out;              index_t vmin_local = vmin;
175  #pragma omp for private(i,j) schedule(static)              index_t vmax_local = vmax;
176              for (j = 0; j < N; j++) {  #pragma omp for
177                  for (i = 0; i < dim; i++)              for (index_t i = 0; i < N; i++) {
178                      if (values[INDEX2(i, j, dim)] != ignore)                  if (values[i] != ignore) {
179                          out_local = std::max(out_local, values[INDEX2(i, j, dim)]);                      vmin_local = std::min(vmin_local, values[i]);
180                        vmax_local = std::max(vmax_local, values[i]);
181                    }
182              }              }
183  #pragma omp critical  #pragma omp critical
184              out = std::max(out_local, out);              {
185                    vmin = std::min(vmin_local, vmin);
186                    vmax = std::max(vmax_local, vmax);
187                }
188          }          }
189      }      }
190      return out;      return IndexPair(vmin,vmax);
191  }  }
192    
193  /* set the index of the positive entries in mask. The length of index is returned. */  std::vector<index_t> packMask(const std::vector<short>& mask)
 dim_t Dudley_Util_packMask(dim_t N, index_t* mask, index_t* index)  
194  {  {
195      dim_t out, k;      std::vector<index_t> index;
196      out = 0;      for (index_t k = 0; k < mask.size(); k++) {
     /*OMP */  
     for (k = 0; k < N; k++) {  
197          if (mask[k] >= 0) {          if (mask[k] >= 0) {
198              index[out] = k;              index.push_back(k);
             out++;  
199          }          }
200      }      }
201      return out;      return index;
202  }  }
203    
204  /* returns true if array contains value */  void setValuesInUse(const int* values, dim_t numValues,
205  bool Dudley_Util_isAny(dim_t N, index_t* array, index_t value)                      std::vector<int>& valuesInUse, escript::JMPI mpiinfo)
206  {  {
207      bool out = false;      const int MAX_VALUE = std::numeric_limits<int>::max();
208      dim_t i;      int lastFoundValue = std::numeric_limits<int>::min();
 #pragma omp parallel for private(i) schedule(static) reduction(||:out)  
     for (i = 0; i < N; i++)  
         out = out || (array[i] == value);  
     return out;  
 }  
   
 /* calculates the cummulative sum in array and returns the total sum */  
 index_t Dudley_Util_cumsum(dim_t N, index_t* array)  
 {  
     index_t out = 0, tmp;  
     dim_t i;  
 #ifdef _OPENMP  
     index_t *partial_sums = NULL, sum;  
     partial_sums = new  index_t[omp_get_max_threads()];  
 #pragma omp parallel private(sum,i,tmp)  
     {  
         sum = 0;  
 #pragma omp for schedule(static)  
         for (i = 0; i < N; ++i)  
             sum += array[i];  
         partial_sums[omp_get_thread_num()] = sum;  
 #pragma omp barrier  
 #pragma omp master  
         {  
             out = 0;  
             for (i = 0; i < omp_get_max_threads(); ++i) {  
                 tmp = out;  
                 out += partial_sums[i];  
                 partial_sums[i] = tmp;  
             }  
         }  
 #pragma omp barrier  
         sum = partial_sums[omp_get_thread_num()];  
 #pragma omp for schedule(static)  
         for (i = 0; i < N; ++i) {  
             tmp = sum;  
             sum += array[i];  
             array[i] = tmp;  
         }  
     }  
     delete[] partial_sums;  
 #else  
     for (i = 0; i < N; ++i) {  
         tmp = out;  
         out += array[i];  
         array[i] = tmp;  
     }  
 #endif  
     return out;  
 }  
   
 void Dudley_Util_setValuesInUse(const index_t* values, dim_t numValues, dim_t* numValuesInUse,  
                                 index_t** valuesInUse, escript::JMPI& mpiinfo)  
 {  
     dim_t i;  
     index_t lastFoundValue = escript::DataTypes::index_t_min(), minFoundValue, local_minFoundValue, *newValuesInUse = NULL;  
     index_t itmp;  
209      bool allFound = false;      bool allFound = false;
210      dim_t nv = 0;  
211        valuesInUse.clear();
212    
213      while (!allFound) {      while (!allFound) {
214          /*          // find smallest value bigger than lastFoundValue
215           *  find smallest value bigger than lastFoundValue          int minFoundValue = MAX_VALUE;
216           */  #pragma omp parallel
217          minFoundValue = escript::DataTypes::index_t_max();          {
218  #pragma omp parallel private(local_minFoundValue)              int local_minFoundValue = minFoundValue;
219          {  #pragma omp for
220              local_minFoundValue = minFoundValue;              for (index_t i = 0; i < numValues; i++) {
221  #pragma omp for private(i,itmp) schedule(static)                  const int val = values[i];
222              for (i = 0; i < numValues; i++) {                  if (val > lastFoundValue && val < local_minFoundValue)
223                  itmp = values[i];                      local_minFoundValue = val;
                 if ((itmp > lastFoundValue) && (itmp < local_minFoundValue))  
                     local_minFoundValue = itmp;  
224              }              }
225  #pragma omp critical  #pragma omp critical
226              {              {
227                  if (local_minFoundValue < minFoundValue)                  if (local_minFoundValue < minFoundValue)
228                      minFoundValue = local_minFoundValue;                      minFoundValue = local_minFoundValue;
229              }              }
   
230          }          }
231  #ifdef ESYS_MPI  #ifdef ESYS_MPI
232          local_minFoundValue = minFoundValue;          int local_minFoundValue = minFoundValue;
233          MPI_Allreduce(&local_minFoundValue, &minFoundValue, 1, MPI_INT, MPI_MIN, mpiinfo->comm);          MPI_Allreduce(&local_minFoundValue, &minFoundValue, 1, MPI_INT,
234                          MPI_MIN, mpiinfo->comm);
235  #endif  #endif
         /* if we found a new tag we need to add this too the valuesInUseList */  
236    
237          if (minFoundValue < escript::DataTypes::index_t_max()) {          // if we found a new value we need to add this to valuesInUse
238              newValuesInUse = new index_t[nv + 1];          if (minFoundValue < MAX_VALUE) {
239              if (*valuesInUse != NULL) {              valuesInUse.push_back(minFoundValue);
                 memcpy(newValuesInUse, *valuesInUse, sizeof(index_t) * nv);  
                 delete[] *valuesInUse;  
             }  
             newValuesInUse[nv] = minFoundValue;  
             *valuesInUse = newValuesInUse;  
             newValuesInUse = NULL;  
             nv++;  
240              lastFoundValue = minFoundValue;              lastFoundValue = minFoundValue;
241          } else {          } else {
242              allFound = true;              allFound = true;
243          }          }
244      }      }
     *numValuesInUse = nv;  
 }  
   
 #ifdef ESYS_MPI  
 void Dudley_printDoubleArray(FILE* fid, dim_t n, double *array, char *name)  
 {  
     index_t i;  
   
     if (name)  
         fprintf(fid, "%s [ ", name);  
     else  
         fprintf(fid, "[ ");  
     for (i = 0; i < (n < 60 ? n : 60); i++)  
         fprintf(fid, "%g ", array[i]);  
     if (n >= 30)  
         fprintf(fid, "... ");  
     fprintf(fid, "]\n");  
 }  
   
 void Dudley_printIntArray(FILE* fid, dim_t n, int *array, char *name)  
 {  
     index_t i;  
   
     if (name)  
         fprintf(fid, "%s [ ", name);  
     else  
         fprintf(fid, "[ ");  
     for (i = 0; i < (n < 60 ? n : 60); i++)  
         fprintf(fid, "%d ", array[i]);  
     if (n >= 30)  
         fprintf(fid, "... ");  
     fprintf(fid, "]\n");  
 }  
   
 void Dudley_printMaskArray(FILE* fid, dim_t n, int *array, char *name)  
 {  
     index_t i;  
   
     if (name)  
         fprintf(fid, "%s [ ", name);  
     else  
         fprintf(fid, "[ ");  
     for (i = 0; i < (n < 60 ? n : 60); i++)  
         if (array[i] != -1)  
             fprintf(fid, "%3d ", array[i]);  
         else  
             fprintf(fid, "  * ");  
     if (n >= 30)  
         fprintf(fid, "... ");  
     fprintf(fid, "]\n");  
245  }  }
 #endif // ESYS_MPI  
246    
247    } // namespace util
248  } // namespace dudley  } // namespace dudley
249    

Legend:
Removed from v.6078  
changed lines
  Added in v.6079

  ViewVC Help
Powered by ViewVC 1.1.26