Changeset 15741


Ignore:
Timestamp:
08/07/13 11:25:05 (12 years ago)
Author:
Mathieu Morlighem
Message:

CHG: keep the exponent as an integer may execute faster (compiler dependent)

Location:
issm/trunk-jpl/src/c
Files:
13 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/classes/Elements/Penta.cpp

    r15738 r15741  
    39383938                        /*Build K: */
    39393939                        GetElementSizes(&hx,&hy,&hz);
    3940                         vel=sqrt(pow(vx,2.)+pow(vy,2.)+pow(vz,2.))+1.e-14;
    3941                         h=sqrt( pow(hx*vx/vel,2.) + pow(hy*vy/vel,2.) + pow(hz*vz/vel,2.));
     3940                        vel=sqrt(vx*vx + vy*vy + vz*vz)+1.e-14;
     3941                        h=sqrt( pow(hx*vx/vel,2) + pow(hy*vy/vel,2) + pow(hz*vz/vel,2));
    39423942                        K[0][0]=h/(2*vel)*vx*vx;  K[0][1]=h/(2*vel)*vx*vy; K[0][2]=h/(2*vel)*vx*vz;
    39433943                        K[1][0]=h/(2*vel)*vy*vx;  K[1][1]=h/(2*vel)*vy*vy; K[1][2]=h/(2*vel)*vy*vz;
    … …  
    41674167                        /*Build K: */
    41684168                        GetElementSizes(&hx,&hy,&hz);
    4169                         vel=sqrt(pow(vx,2.)+pow(vy,2.)+pow(vz,2.))+1.e-14;
    4170                         h=sqrt( pow(hx*vx/vel,2.) + pow(hy*vy/vel,2.) + pow(hz*vz/vel,2.));
     4169                        vel=sqrt(vx*vx + vy*vy + vz*vz)+1.e-14;
     4170                        h=sqrt( pow(hx*vx/vel,2) + pow(hy*vy/vel,2) + pow(hz*vz/vel,2));
    41714171
    41724172                        K[0][0]=h/(2*vel)*fabs(vx*vx);  K[0][1]=h/(2*vel)*fabs(vx*vy); K[0][2]=h/(2*vel)*fabs(vx*vz);
    … …  
    84968496                /*Deal with lower surface*/
    84978497                if(IsOnBed()){
    8498                         constant_part=-1.58*pow(10.,-10.)*rho_ice*gravity*thickness;
     8498                        constant_part=-1.58*1.e-10*rho_ice*gravity*thickness;
    84998499                        ub=constant_part*slope[0];
    85008500                        vb=constant_part*slope[1];

  • issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp

    r15694 r15741  
    18741874        /*Jdet = (Area of the trapezoid)/(Area trapezoid ref) with AreaRef = 4*/
    18751875        /*Area of a trabezoid = altitude * (base1 + base2)/2 */
    1876         *Jdet= pow(pow(x2-x1,2.) + pow(y2-y1,2.),0.5) * (z4-z1 + z3-z2)/8;
    1877         if(*Jdet<0) _error_("negative jacobian determinant!");
     1876        *Jdet= sqrt((x2-x1)*(x2-x1) + (y2-y1)*(y2-y1)) * (z4-z1 + z3-z2)/8.;
     1877        if(*Jdet<0.) _error_("negative jacobian determinant!");
    18781878
    18791879}

  • issm/trunk-jpl/src/c/classes/Elements/Tria.cpp

    r15732 r15741  
    33213321                constant_part=-2.*pow(rho_ice*gravity,n)*pow(slope2,((n-1.)/2.));
    33223322
    3323                 ub=-1.58*pow(10.,-10.)*rho_ice*gravity*thickness*slope[0];
    3324                 vb=-1.58*pow(10.,-10.)*rho_ice*gravity*thickness*slope[1];
     3323                ub=-1.58*1.e-10*rho_ice*gravity*thickness*slope[0];
     3324                vb=-1.58*1.e-10*rho_ice*gravity*thickness*slope[1];
    33253325
    33263326                pe->values[2*i+0]=(ub-2.*pow(rho_ice*gravity,n)*pow(slope2,((n-1.)/2.))*pow(thickness,n)/(pow(B,n)*(n+1))*slope[0])/connectivity;
    … …  
    44294429
    44304430                /*Tikhonov regularization: J = 1/2 ((dp/dx)^2 + (dp/dy)^2) */
    4431                 Jelem+=weight*1/2*(pow(dp[0],2.)+pow(dp[1],2.))*Jdet*gauss->weight;
     4431                Jelem+=weight*1/2*(dp[0]*dp[0] + dp[1]*dp[1])*Jdet*gauss->weight;
    44324432        }
    44334433

  • issm/trunk-jpl/src/c/classes/kriging/ExponentialVariogram.cpp

    r15104 r15741  
    6565
    6666        /*Calculate length*/
    67         h=sqrt(pow(deltax,2.)+pow(deltay,2.));
     67        h=sqrt(deltax*deltax + deltay*deltay);
    6868
    6969        /*If h is too small, return sill*/
    … …  
    8282
    8383        /*Calculate length*/
    84         h=sqrt(pow(deltax,2.)+pow(deltay,2.));
     84        h=sqrt(deltax*deltax + deltay*deltay);
    8585
    8686        /*return semi-variogram*/

  • issm/trunk-jpl/src/c/classes/kriging/GaussianVariogram.cpp

    r15104 r15741  
    6565
    6666        /*Calculate length square*/
    67         h2=pow(deltax,2.)+pow(deltay,2.);
     67        h2=deltax*deltax + deltay*deltay;
    6868
    6969        /*If h is too small, return sill*/
    … …  
    8383
    8484        /*Calculate length square*/
    85         h2=pow(deltax,2.)+pow(deltay,2.);
     85        h2=deltax*deltax + deltay*deltay;
    8686
    8787        /*return semi-variogram*/

  • issm/trunk-jpl/src/c/classes/kriging/Observations.cpp

    r15104 r15741  
    563563                        observation2=dynamic_cast<Observation*>(this->GetObjectByOffset(j));
    564564
    565                         distance=sqrt(pow(observation1->x - observation2->x,2.) + pow(observation1->y - observation2->y,2.));
     565                        distance=sqrt(pow(observation1->x - observation2->x,2) + pow(observation1->y - observation2->y,2));
    566566                        if(distance>x[n-1]) continue;
    567567
    … …  
    570570                        if(index<0)   index = 0;
    571571
    572                         gamma[index]   += 1./2.*pow(observation1->value - observation2->value,2.);
     572                        gamma[index]   += 1./2.*pow(observation1->value - observation2->value,2);
    573573                        counter[index] += 1.;
    574574                }

  • issm/trunk-jpl/src/c/classes/kriging/PowerVariogram.cpp

    r15104 r15741  
    6666
    6767        /*Calculate length square*/
    68         h=sqrt(pow(deltax,2.)+pow(deltay,2.));
     68        h=sqrt(deltax*deltax + deltay*deltay);
    6969
    7070        /*return covariance*/
    … …  
    8080
    8181        /*Calculate length square*/
    82         h=sqrt(pow(deltax,2.)+pow(deltay,2.));
     82        h=sqrt(deltax*deltax + deltay*deltay);
    8383
    8484        /*return semi-variogram*/

  • issm/trunk-jpl/src/c/classes/kriging/Quadtree.cpp

    r15104 r15741  
    230230        if(box && box->nbitems==4){
    231231                index  = 0;
    232                 length = pow(box->obs[0]->x - x,2.) + pow(box->obs[0]->y - y,2.);
     232                length = pow(box->obs[0]->x - x,2) + pow(box->obs[0]->y - y,2);
    233233                for(int i=1;i<4;i++){
    234                         length2 = pow(box->obs[i]->x - x,2.) + pow(box->obs[i]->y - y,2.);
     234                        length2 = pow(box->obs[i]->x - x,2) + pow(box->obs[i]->y - y,2);
    235235                        if(length2<length){
    236236                                index  = i;
    … …  
    282282        if(box && box->nbitems>0){
    283283                index  = box->obs[0]->index;
    284                 length = pow(box->obs[0]->x - x,2.) + pow(box->obs[0]->y - y,2.);
     284                length = pow(box->obs[0]->x - x,2) + pow(box->obs[0]->y - y,2);
    285285                for(int i=1;i<box->nbitems;i++){
    286                         length2 = pow(box->obs[i]->x - x,2.) + pow(box->obs[i]->y - y,2.);
     286                        length2 = pow(box->obs[i]->x - x,2) + pow(box->obs[i]->y - y,2);
    287287                        if(length2<length){
    288288                                index  = box->obs[i]->index;

  • issm/trunk-jpl/src/c/classes/kriging/SphericalVariogram.cpp

    r15104 r15741  
    6565
    6666        /*Calculate length square*/
    67         h=sqrt(pow(deltax,2.)+pow(deltay,2.));
     67        h=sqrt(deltax*deltax + deltay*deltay);
    6868
    6969        /*return covariance*/
    7070        if(h<=range)
    71          cova = (sill-nugget)*(1 - (3*h)/(2*range) + pow(h,3.)/(2*pow(range,3.)) );
     71         cova = (sill-nugget)*(1 - (3*h)/(2*range) + pow(h,3)/(2*pow(range,3)) );
    7272        else
    7373         cova = 0.;
    … …  
    8282
    8383        /*Calculate length square*/
    84         h=sqrt(pow(deltax,2.)+pow(deltay,2.));
     84        h=sqrt(deltax*deltax + deltay*deltay);
    8585
    8686        /*return semi-variogram*/
    8787        if(h<=range)
    88          gamma = (sill-nugget)*( (3*h)/(2*range) - pow(h,3.)/(2*pow(range,3.)) ) + nugget;
     88         gamma = (sill-nugget)*( (3*h)/(2*range) - pow(h,3)/(2*pow(range,3)) ) + nugget;
    8989        else
    9090         gamma = sill;

  • issm/trunk-jpl/src/c/shared/Elements/Arrhenius.cpp

    r14975 r15741  
    2626
    2727        /*Some physical constants (Payne2000)*/
    28         IssmDouble beta=8.66*pow(10.,-4.);
     28        IssmDouble beta=8.66*1.e-4;
    2929        IssmDouble R=8.314;
    3030
    … …  
    3939        /*Get A*/
    4040        if(Tstar<263.15){
    41                 A=3.61*pow(10.,-13.) * exp(  -6.*pow(10.,4.)/(R*Tstar));
     41                A=3.61e-13*exp(  -6.e+4/(R*Tstar));
    4242        }
    4343        else{
    44                 A=1.73*pow(10.,  3.) * exp(-13.9*pow(10.,4.)/(R*Tstar));
     44                A=1.73e+3 *exp(-13.9e+4/(R*Tstar));
    4545        }
    4646

  • issm/trunk-jpl/src/c/shared/LatLong/Ll2xyx.cpp

    r15104 r15741  
    5858        cde = 57.29577951;
    5959        /*  Radius of the earth in meters  */
    60         re  = 6378.273*pow(10.,3.);
     60        re  = 6378.273*1.e3;
    6161        /*  Eccentricity of the Hughes ellipsoid squared  */
    6262        ex2 = 0.006693883;
    … …  
    7777                        sl  = slat*PI/180.;
    7878                        tc  = tan(PI/4.-sl/2.)/pow(((1.-ex*sin(sl))/(1.+ex*sin(sl))),(ex/2.));
    79                         mc  = cos(sl)/sqrt(1.0-ex2*(pow(sin(sl),2.)));
     79                        mc  = cos(sl)/sqrt(1.0-ex2*(pow(sin(sl),2)));
    8080                        rho = re*mc*T/tc;
    8181                }

  • issm/trunk-jpl/src/c/shared/LatLong/Xy2llx.cpp

    r15104 r15741  
    5959        cde  = 57.29577951;
    6060        /*  Radius of the earth in meters  */
    61         re   = 6378.273*pow(10.,3.);
     61        re   = 6378.273e+3;
    6262        /*  Eccentricity of the Hughes ellipsoid squared  */
    6363        ex2  = 0.006693883;
    … …  
    6868        for(i=0; i<ncoord; i++){
    6969                sl = slat*PI/180.;
    70                 cm = cos(sl)/sqrt(1.0-ex2*(pow(sin(sl),2.)));
    71                 rho= sqrt(pow(x[i],2.) + pow(y[i],2.));
     70                cm = cos(sl)/sqrt(1.0-ex2*(pow(sin(sl),2)));
     71                rho= sqrt(pow(x[i],2) + pow(y[i],2));
    7272                T  = tan((PI/4.0) - (sl/2.0))/pow(((1.0-ex*sin(sl))/(1.0+ex*sin(sl))),(ex/2.0));
    7373

  • issm/trunk-jpl/src/c/shared/Numerics/GaussPoints.cpp

    r15104 r15741  
    14521452                                          *(5.+3.*egaus[j]-zgaus[k]-3.*egaus[j]*zgaus[k]);
    14531453                                        zeta      =sqrt(2./3.)*(1.+zgaus[k]);
    1454                                         (*pwgt)[ipt]=xwgt[i]*ewgt[j]*zwgt[k]
    1455                                           *(SQRT2/16.
    1456                                                                   *(1.-egaus[j])*pow(1.-zgaus[k],2.));
     1454                                        (*pwgt)[ipt]=xwgt[i]*ewgt[j]*zwgt[k]*(SQRT2/16.*(1.-egaus[j])*pow(1.-zgaus[k],2));
    14571455
    14581456                                        (*pl1 )[ipt]=(1.-xi-eta/SQRT3-zeta/sqrt(6.))/2.;
Note: See TracChangeset for help on using the changeset viewer.