AMROC - Source code documentation

00001 #ifndef _included_GridFunctionOpsRed1_h
00002 #define _included_GridFunctionOpsRed1_h
00003 
00009 /*************************************************************************/
00010 /* Maxval */
00011 /*************************************************************************/
00012 template <class Type>
00013 Type GridFunction(1)<Type>::GF_maxval(const int time, const int level)
00014   {
00015    register int t = dagh_timeindex(time,level);
00016    register int l = level;
00017    Type maxval = (Type) DAGHSmall, tmax = (Type) DAGHSmall;
00018    register int i;
00019    short flag = DAGHTrue;
00020    for (i=0; i<length[l]; i++) if (gdb[t][l][i]) {
00021      tmax = (gdb[t][l][i]->griddata()).maxval(gdb[t][l][i]->interiorbox());
00022      if (flag == DAGHTrue) { maxval = tmax; flag = DAGHFalse; }
00023      else if (maxval < tmax) maxval = tmax;
00024      tmax = (Type) DAGHSmall;
00025    }
00026 
00027 #ifdef DAGH_NO_MPI
00028 #else
00029    if (comm_service::dce() && comm_service::proc_num() > 1) { 
00030 
00031      int num = comm_service::proc_num();
00032  
00033      int R;
00034      int size = sizeof(Type);
00035      void *values = (void *) new char[size*num];
00036 
00037      R = MPI_Allgather(&maxval, size, MPI_BYTE, values, size, MPI_BYTE, 
00038                        comm_service::comm(gfid));
00039      if ( MPI_SUCCESS != R ) 
00040         comm_service::error_die( "GridFunction::GF_maxval", "MPI_Allgather", R );
00041 
00042 #ifdef DEBUG_PRINT_GF_COMM
00043    ( comm_service::log() << "GridFunction::GF_maxval" << " " << comm_service::proc_me() << " "
00044                          << "MPI_Allgather { "
00045                          << size
00046                          << " }"
00047                          << endl ).flush();
00048 #endif
00049 
00050      maxval = (Type) DAGHSmall;
00051      flag = DAGHTrue;
00052      for (i=0;i<num;i++) {
00053        tmax = *((Type *)values + i);
00054        if (flag == DAGHTrue) { maxval = tmax; flag = DAGHFalse; }
00055        else if (maxval < tmax) maxval = tmax;
00056      }
00057      if (values) delete [] (char *)values;
00058    }
00059 #endif
00060    
00061    return (maxval);
00062   }
00063 
00064 template <class Type>
00065 Type GridFunction(1)<Type>::GF_maxval(const int time, const int level, const BBox &where)
00066   {
00067    assert (where.rank == dagh.rank);
00068 
00069    BBox alignedwhere = where;
00070    gdbAlignBBox(gfrank,alignedwhere,alignment);
00071 
00072    register int t = dagh_timeindex(time,level);
00073    register int l = level;
00074    Type maxval = (Type) DAGHSmall, tmax = (Type) DAGHSmall;
00075    register int i;
00076    short flag = DAGHTrue;
00077    for (i=0; i<length[l]; i++) if (gdb[t][l][i]) {
00078      tmax = (gdb[t][l][i]->griddata()).maxval(
00079                            alignedwhere*gdb[t][l][i]->interiorbox());
00080      if (flag == DAGHTrue) { maxval = tmax; flag = DAGHFalse; }
00081      else if (maxval < tmax) maxval = tmax;
00082      tmax = (Type) DAGHSmall;
00083    }
00084 
00085 #ifdef DAGH_NO_MPI
00086 #else
00087    if (comm_service::dce() && comm_service::proc_num() > 1) { 
00088 
00089      int num = comm_service::proc_num();
00090  
00091      int R;
00092      int size = sizeof(Type);
00093      void *values = (void *) new char[size*num];
00094 
00095      R = MPI_Allgather(&maxval, size, MPI_BYTE, values, size, MPI_BYTE, 
00096                        comm_service::comm(gfid));
00097      if ( MPI_SUCCESS != R ) 
00098         comm_service::error_die( "GridFunction::GF_maxval", "MPI_Allgather", R );
00099 
00100 #ifdef DEBUG_PRINT_GF_COMM
00101    ( comm_service::log() << "GridFunction::GF_maxval" << " " << comm_service::proc_me() << " "
00102                          << "MPI_Allgather { "
00103                          << size
00104                          << " }"
00105                          << endl ).flush();
00106 #endif
00107 
00108      maxval = (Type) DAGHSmall;
00109      flag = DAGHTrue;
00110      for (i=0;i<num;i++) {
00111        tmax = *((Type *)values + i);
00112        if (flag == DAGHTrue) { maxval = tmax; flag = DAGHFalse; }
00113        else if (maxval < tmax) maxval = tmax;
00114      }
00115      if (values) delete [] (char *)values;
00116    }
00117 #endif
00118 
00119    return (maxval);
00120   }
00121 
00122 /*************************************************************************/
00123 /* Minval */
00124 /*************************************************************************/
00125 template <class Type>
00126 Type GridFunction(1)<Type>::GF_minval(const int time, const int level)
00127   {
00128    register int t = dagh_timeindex(time,level);
00129    register int l = level;
00130    Type minval = MAXLIM(Type), tmin = MAXLIM(Type);
00131    register int i;
00132    short flag = DAGHTrue;
00133    for (i=0; i<length[l]; i++) if (gdb[t][l][i]) {  
00134      tmin = (gdb[t][l][i]->griddata()).minval(gdb[t][l][i]->interiorbox());
00135      if (flag == DAGHTrue) { minval = tmin; flag = DAGHFalse; }
00136      else if (minval > tmin) minval = tmin;
00137      tmin = MAXLIM(Type);
00138    }
00139 
00140 #ifdef DAGH_NO_MPI
00141 #else
00142    if (comm_service::dce() && comm_service::proc_num() > 1) { 
00143 
00144      int num = comm_service::proc_num();
00145  
00146      int R;
00147      int size = sizeof(Type);
00148      void *values = (void *) new char[size*num];
00149 
00150      R = MPI_Allgather(&minval, size, MPI_BYTE, values, size, MPI_BYTE, 
00151                        comm_service::comm(gfid));
00152      if ( MPI_SUCCESS != R ) 
00153         comm_service::error_die( "GridFunction::GF_minval", "MPI_Allgather", R );
00154 
00155 #ifdef DEBUG_PRINT_GF_COMM
00156    ( comm_service::log() << "GridFunction::GF_minval" << " " << comm_service::proc_me() << " "
00157                          << "MPI_Allgather { "
00158                          << size
00159                          << " }"
00160                          << endl ).flush();
00161 #endif
00162 
00163      minval = MAXLIM(Type);
00164      flag = DAGHTrue;
00165      for (i=0;i<num;i++) {
00166        tmin = *((Type *)values + i);
00167        if (flag == DAGHTrue) { minval = tmin; flag = DAGHFalse; }
00168        else if (minval > tmin) minval = tmin;
00169      }
00170      if (values) delete [] (char *)values;
00171    }
00172 #endif
00173 
00174    return (minval);
00175   }
00176 
00177 template <class Type>
00178 Type GridFunction(1)<Type>::GF_minval(const int time, const int level, const BBox &where)
00179   {
00180    assert (where.rank == dagh.rank);
00181 
00182    BBox alignedwhere = where;
00183    gdbAlignBBox(gfrank,alignedwhere,alignment);
00184 
00185    register int t = dagh_timeindex(time,level);
00186    register int l = level;
00187    Type minval = MAXLIM(Type), tmin = MAXLIM(Type);
00188    register int i;
00189    short flag = DAGHTrue;
00190    for (i=0; i<length[l]; i++) if (gdb[t][l][i]) { 
00191      tmin = (gdb[t][l][i]->griddata()).minval(
00192                            alignedwhere*gdb[t][l][i]->interiorbox());
00193      if (flag == DAGHTrue) { minval = tmin; flag = DAGHFalse; }     
00194      else if (minval > tmin) minval = tmin;
00195      tmin = MAXLIM(Type);
00196    }
00197 
00198 #ifdef DAGH_NO_MPI
00199 #else
00200 
00201    if (comm_service::dce() && comm_service::proc_num() > 1) { 
00202 
00203      int num = comm_service::proc_num();
00204  
00205      int R;
00206      int size = sizeof(Type);
00207      void *values = (void *) new char[size*num];
00208 
00209      R = MPI_Allgather(&minval, size, MPI_BYTE, values, size, MPI_BYTE, 
00210                        comm_service::comm(gfid));
00211      if ( MPI_SUCCESS != R ) 
00212         comm_service::error_die( "GridFunction::GF_minval", "MPI_Allgather", R );
00213 
00214 #ifdef DEBUG_PRINT_GF_COMM
00215    ( comm_service::log() << "GridFunction::GF_minval" << " " << comm_service::proc_me() << " "
00216                          << "MPI_Allgather { "
00217                          << size
00218                          << " }"
00219                          << endl ).flush();
00220 #endif
00221 
00222      minval = MAXLIM(Type);
00223      flag = DAGHTrue;
00224      for (i=0;i<num;i++) {
00225        tmin = *((Type *)values + i);
00226        if (flag == DAGHTrue) { minval = tmin; flag = DAGHFalse; }
00227        else if (minval > tmin) minval = tmin;
00228      }
00229      if (values) delete [] (char *)values;
00230    }
00231 #endif
00232 
00233    return (minval);
00234   }
00235 
00236 /*************************************************************************/
00237 /* Sum */
00238 /*************************************************************************/
00239 template <class Type>
00240 Type GridFunction(1)<Type>::GF_sum(const int time, const int level)
00241   {
00242    register int t = dagh_timeindex(time,level);
00243    register int l = level;
00244    Type sum = (Type) 0;
00245    register int i;
00246    for (i=0; i<length[l]; i++) if (gdb[t][l][i]) { 
00247      sum += (gdb[t][l][i]->griddata()).sum(gdb[t][l][i]->interiorbox());
00248    }
00249 
00250 #ifdef DAGH_NO_MPI
00251 #else
00252    if (comm_service::dce() && comm_service::proc_num() > 1) { 
00253 
00254      int num = comm_service::proc_num();
00255  
00256      int R;
00257      int size = sizeof(Type);
00258      void *values = (void *) new char[size*num];
00259 
00260      R = MPI_Allgather(&sum, size, MPI_BYTE, values, size, MPI_BYTE, 
00261                        comm_service::comm(gfid));
00262      if ( MPI_SUCCESS != R ) 
00263         comm_service::error_die( "GridFunction::GF_sum", "MPI_Allgather", R );
00264 
00265 #ifdef DEBUG_PRINT_GF_COMM
00266    ( comm_service::log() << "GridFunction::GF_sum" << " " << comm_service::proc_me() << " "
00267                          << "MPI_Allgather { "
00268                          << size
00269                          << " }"
00270                          << endl ).flush();
00271 #endif
00272 
00273      sum = (Type) 0;
00274      for (i=0;i<num;i++)
00275        sum += *((Type *)values + i);
00276 
00277      if (values) delete [] (char *)values;
00278    }
00279 #endif
00280 
00281    return (sum);
00282   }
00283 
00284 template <class Type>
00285 Type GridFunction(1)<Type>::GF_sum(const int time, const int level, const BBox &where)
00286   {
00287    assert (where.rank == dagh.rank);
00288 
00289    BBox alignedwhere = where;
00290    gdbAlignBBox(gfrank,alignedwhere,alignment);
00291 
00292    register int t = dagh_timeindex(time,level);
00293    register int l = level;
00294    Type sum = (Type) 0;
00295    register int i;
00296    for (i=0; i<length[l]; i++) if (gdb[t][l][i]) {  
00297      sum += (gdb[t][l][i]->griddata()).sum(
00298                            alignedwhere*gdb[t][l][i]->interiorbox());
00299    }
00300 
00301 #ifdef DAGH_NO_MPI
00302 #else
00303    if (comm_service::dce() && comm_service::proc_num() > 1) { 
00304 
00305      int num = comm_service::proc_num();
00306  
00307      int R;
00308      int size = sizeof(Type);
00309      void *values = (void *) new char[size*num];
00310 
00311      R = MPI_Allgather(&sum, size, MPI_BYTE, values, size, MPI_BYTE, 
00312                        comm_service::comm(gfid));
00313      if ( MPI_SUCCESS != R ) 
00314         comm_service::error_die( "GridFunction::GF_sum", "MPI_Allgather", R );
00315 
00316 #ifdef DEBUG_PRINT_GF_COMM
00317    ( comm_service::log() << "GridFunction::GF_sum" << " " << comm_service::proc_me() << " "
00318                          << "MPI_Allgather { "
00319                          << size
00320                          << " }"
00321                          << endl ).flush();
00322 #endif
00323 
00324      sum = (Type) 0;
00325      for (i=0;i<num;i++)
00326        sum += *((Type *)values + i);
00327    
00328      if (values) delete [] (char *)values;
00329    }
00330 #endif
00331 
00332    return (sum);
00333   }
00334 
00335 /*************************************************************************/
00336 /* Product */
00337 /*************************************************************************/
00338 template <class Type>
00339 Type GridFunction(1)<Type>::GF_product(const int time, const int level)
00340   {
00341    register int t = dagh_timeindex(time,level);
00342    register int l = level;
00343    Type prod = (Type) 1;
00344    register int i;
00345    for (i=0; i<length[l]; i++) {
00346      if (gdb[t][l][i]) 
00347        prod *= (gdb[t][l][i]->griddata()).product(gdb[t][l][i]->interiorbox());
00348    }
00349 
00350 #ifdef DAGH_NO_MPI
00351 #else
00352    if (comm_service::dce() && comm_service::proc_num() > 1) { 
00353 
00354      int num = comm_service::proc_num();
00355  
00356      int R;
00357      int size = sizeof(Type);
00358      void *values = (void *) new char[size*num];
00359 
00360      R = MPI_Allgather(&prod, size, MPI_BYTE, values, size, MPI_BYTE, 
00361                        comm_service::comm(gfid));
00362      if ( MPI_SUCCESS != R ) 
00363         comm_service::error_die( "GridFunction::GF_product", "MPI_Allgather", R );
00364 
00365 #ifdef DEBUG_PRINT_GF_COMM
00366    ( comm_service::log() << "GridFunction::GF_product" << " " << comm_service::proc_me() << " "
00367                          << "MPI_Allgather { "
00368                          << size
00369                          << " }"
00370                          << endl ).flush();
00371 #endif
00372 
00373      prod = (Type) 1;
00374      for (i=0;i<num;i++)
00375        prod *= *((Type *)values + i);
00376    
00377      if (values) delete [] (char *)values;
00378    }
00379 #endif
00380 
00381    return (prod);
00382   }
00383 
00384 template <class Type>
00385 Type GridFunction(1)<Type>::GF_product(const int time, const int level, const BBox &where)
00386   {
00387    assert (where.rank == dagh.rank);
00388 
00389    BBox alignedwhere = where;
00390    gdbAlignBBox(gfrank,alignedwhere,alignment);
00391 
00392    register int t = dagh_timeindex(time,level);
00393    register int l = level;
00394    Type prod = (Type) 1;
00395    register int i;
00396    for (i=0; i<length[l]; i++) {
00397      if (gdb[t][l][i]) 
00398        prod *= (gdb[t][l][i]->griddata()).product(
00399                               alignedwhere*gdb[t][l][i]->interiorbox());
00400    }
00401 
00402 #ifdef DAGH_NO_MPI
00403 #else
00404    if (comm_service::dce() && comm_service::proc_num() > 1) { 
00405 
00406      int num = comm_service::proc_num();
00407  
00408      int R;
00409      int size = sizeof(Type);
00410      void *values = (void *) new char[size*num];
00411 
00412      R = MPI_Allgather(&prod, size, MPI_BYTE, values, size, MPI_BYTE, 
00413                        comm_service::comm(gfid));
00414      if ( MPI_SUCCESS != R ) 
00415         comm_service::error_die( "GridFunction::GF_product", "MPI_Allgather", R );
00416 
00417 #ifdef DEBUG_PRINT_GF_COMM
00418    ( comm_service::log() << "GridFunction::GF_product" << " " << comm_service::proc_me() << " "
00419                          << "MPI_Allgather { "
00420                          << size
00421                          << " }"
00422                          << endl ).flush();
00423 #endif
00424 
00425      prod = (Type) 1;
00426      for (i=0;i<num;i++)
00427        prod *= *((Type *)values + i);
00428    
00429      if (values) delete [] (char *)values;
00430    }
00431 #endif
00432 
00433    return (prod);
00434   }
00435 
00436 /*************************************************************************/
00437 /* Norms */
00438 /* Exclusion of higher levels would make sense */
00439 /*************************************************************************/
00440 template <class Type>
00441 double GridFunction(1)<Type>::GF_norm(const int time, const int level, const int param) {
00442 
00443   register int t = dagh_timeindex(time,level);
00444   register int l = level;
00445   double s = 0.0;
00446   register int i;
00447   for (i=0; i<length[l]; i++) 
00448     if (gdb[t][l][i]) {
00449       switch (param) {
00450       case DAGHNormL1:
00451         s += (gdb[t][l][i]->griddata()).sumabs(gdb[t][l][i]->interiorbox());
00452         break;
00453       case DAGHNormL2:
00454         s += (gdb[t][l][i]->griddata()).sumsqrd(gdb[t][l][i]->interiorbox());
00455         break;
00456       case DAGHNormMax: 
00457         {
00458           double tm = (gdb[t][l][i]->griddata()).maxabs(gdb[t][l][i]->interiorbox());
00459           s = Max(tm, s);
00460           break;
00461         }
00462       default:
00463         break;
00464       }
00465     }
00466 
00467 #ifdef DAGH_NO_MPI
00468 #else
00469   if (comm_service::dce() && comm_service::proc_num() > 1) { 
00470 
00471     int num = comm_service::proc_num();
00472     
00473     int R;
00474     int size = sizeof(double);
00475     void *values = (void *) new char[size*num];
00476     
00477     R = MPI_Allgather(&s, size, MPI_BYTE, values, size, MPI_BYTE, 
00478                       comm_service::comm(gfid));
00479     if ( MPI_SUCCESS != R ) 
00480       comm_service::error_die( "GridFunction::GF_norm", "MPI_Allgather", R );
00481      
00482 #ifdef DEBUG_PRINT_GF_COMM
00483     ( comm_service::log() << "GridFunction::GF_norm" << " " << comm_service::proc_me() << " "
00484                           << "MPI_Allgather { "
00485                           << size
00486                           << " }"
00487                           << endl ).flush();
00488 #endif
00489     s = 0.0;
00490     for (i=0;i<num;i++) {
00491       switch (param) {
00492       case DAGHNormL1:
00493       case DAGHNormL2:
00494         s += *((double *)values + i);
00495         break;
00496       case DAGHNormMax: 
00497         s = Max(s, *((double *)values + i));
00498         break;
00499       default:
00500         break;
00501       }
00502     }
00503    
00504     if (values) delete [] (char *)values;
00505   }
00506 #endif
00507    
00508   int d;
00509   switch (param) {
00510   case DAGHNormL1:
00511     for(d=0; d<gfrank; d++)
00512       s *= DeltaX(dagh,d,level);
00513     break;
00514   case DAGHNormL2:
00515     for(d=0; d<gfrank; d++)
00516       s *= DeltaX(dagh,d,level);
00517     s = sqrt(s);
00518     break;
00519   default:
00520     break;
00521   }
00522 
00523   return s;
00524 }
00525 
00526 template <class Type>
00527 double GridFunction(1)<Type>::GF_norm(const int time, const int level,
00528                                      const BBox &where, const int param) {
00529   assert (where.rank == dagh.rank);
00530 
00531   BBox alignedwhere = where;
00532   gdbAlignBBox(gfrank,alignedwhere,alignment);
00533   
00534   register int t = dagh_timeindex(time,level);
00535   register int l = level;
00536   double s = 0.0; 
00537   register int i;
00538   for (i=0; i<length[l]; i++) 
00539     if (gdb[t][l][i]) {
00540       switch (param) {
00541       case DAGHNormL1:
00542         s += (gdb[t][l][i]->griddata()).sumabs(
00543               alignedwhere*gdb[t][l][i]->interiorbox());
00544         break;
00545       case DAGHNormL2:
00546         s += (gdb[t][l][i]->griddata()).sumsqrd(
00547               alignedwhere*gdb[t][l][i]->interiorbox());
00548         break;
00549       case DAGHNormMax: 
00550         {
00551           double tm = (gdb[t][l][i]->griddata()).maxabs(
00552                       alignedwhere*gdb[t][l][i]->interiorbox());
00553           s = Max(tm, s);
00554           break;
00555         }
00556       default:
00557         break;
00558       }
00559     }
00560 
00561 #ifdef DAGH_NO_MPI
00562 #else
00563   if (comm_service::dce() && comm_service::proc_num() > 1) { 
00564     
00565     int num = comm_service::proc_num();
00566     
00567     int R;
00568     int size = sizeof(double);
00569     void *values = (void *) new char[size*num];
00570     
00571     R = MPI_Allgather(&s, size, MPI_BYTE, values, size, MPI_BYTE, 
00572                       comm_service::comm(gfid));
00573     if ( MPI_SUCCESS != R ) 
00574       comm_service::error_die( "GridFunction::GF_norm", "MPI_Allgather", R );
00575 
00576 #ifdef DEBUG_PRINT_GF_COMM
00577     ( comm_service::log() << "GridFunction::GF_norm" << " " << comm_service::proc_me() << " "
00578                           << "MPI_Allgather { "
00579                           << size
00580                           << " }"
00581                           << endl ).flush();
00582 #endif
00583     s = 0.0;
00584     for (i=0;i<num;i++) {
00585       switch (param) {
00586       case DAGHNormL1:
00587       case DAGHNormL2:
00588         s += *((double *)values + i);
00589         break;
00590       case DAGHNormMax: 
00591         s = Max(s, *((double *)values + i));
00592         break;
00593       default:
00594         break;
00595       }
00596     }
00597     
00598     if (values) delete [] (char *)values;
00599   }
00600 #endif
00601    
00602   int d;
00603   switch (param) {
00604   case DAGHNormL1:
00605     for(d=0; d<gfrank; d++)
00606       s *= DeltaX(dagh,d,level);
00607     break;
00608   case DAGHNormL2:
00609     for(d=0; d<gfrank; d++)
00610       s *= DeltaX(dagh,d,level);
00611     s = sqrt(s);
00612     break;
00613   default:
00614     break;
00615   }
00616   
00617   return s;
00618 }
00619 
00620 
00621 #endif
00622
GridFunctionOpsRed1.h

last update: 06/01/04
