Changeset 13061


Ignore:
Timestamp:
08/16/12 09:24:06 (13 years ago)
Author:
Mathieu Morlighem
Message:

CHG: ISSMized cubic.cpp

Location:
issm/trunk-jpl/src/c/shared/Numerics
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • issm/trunk-jpl/src/c/shared/Numerics/cubic.cpp

    r13050 r13061  
    22#include "./numerics.h"
    33
    4 int signR(double Z){
    5         int ret;
    6         if (Z > 0.0)    ret = 1;
    7         if (Z < 0.0)    ret = -1;
    8         if (Z == 0.0)   ret =0;
     4int cubic(IssmDouble a,IssmDouble b,IssmDouble c,IssmDouble d, IssmDouble x[3], int* num){
     5        /* Find the real roots of linear/quadratic and cubic equations:
     6         *
     7         *   a x^3 + bx^2 + cx + d = 0
     8         *
     9         *   returns the roots in x
     10         *   num is the number of roots */
    911
    10         return ret;
    11 }
     12        /*Some useful constants*/
     13        const IssmDouble pi    = 3.1415926535897932;
     14        const IssmDouble third = 1./3.;
    1215
    13 double CBRT(double Z){
    14         double ret;
    15         const double THIRD = 1./3.;
    16         //define cubic root as statement function
    17         //SIGN has different meanings in both C and Fortran
    18         // Was unable to use the sign command of C, so wrote my own
    19         // that why a new variable needs to be introduced that keeps track of the sign of
    20         // SIGN is supposed to return a 1, -1 or 0 depending on what the sign of the argument is
    21         ret = fabs(pow(fabs(Z),THIRD)) * signR(Z);
    22         return ret;
    23 }
     16        /*Intermediaries*/
     17        IssmDouble U[3],W,P,Q,delta,phi;
    2418
    25 int cubic(double A[4], double X[3], int* num){
    26         const double PI = 3.1415926535897932;
    27         const double THIRD = 1./3.;
    28         double U[3],W, P, Q, DIS, PHI;
    29         int i;
    30 
    31         //define cubic root as statement function
    32         // In C, the function is defined outside of the cubic fct
    33 
    34         // ====determine the degree of the polynomial ====
    35 
    36         if (A[3] != 0.0){
     19        /* determine the degree of the polynomial */
     20        if (a != 0.0){
    3721                //cubic problem
    38                 W = A[2]/A[3]*THIRD;
    39                 P = pow((A[1]/A[3]*THIRD - pow(W,2)),3);
    40                 Q = -.5*(2.0*pow(W,3)-(A[1]*W-A[0])/A[3] );
    41                 DIS = pow(Q,2)+P;
    42                 if ( DIS < 0.0 ){
     22                W     = b/a *third;
     23                P     = pow((c/a *third - pow(W,2)),3);
     24                Q     = -.5 *(2.0*pow(W,3)-(c*W-d)/a );
     25                delta = pow(Q,2)+P;
     26                if ( delta < 0.0 ){
    4327                        //three real solutions!
    44                         //Confine the argument of ACOS to the interval [-1;1]!
    45                         PHI = acos(min(1.0,max(-1.0,Q/sqrt(-P))));
    46                         P=2.0*pow((-P),(5.e-1*THIRD));
    47                         for (i=0;i<3;i++)       U[i] = P*cos((PHI+2*((double)i)*PI)*THIRD)-W;
    48                         X[0] = min(U[0], min(U[1], U[2]));
    49                         X[1] = max(min(U[0], U[1]),max( min(U[0], U[2]), min(U[1], U[2])));
    50                         X[2] = max(U[0], max(U[1], U[2]));
     28                        //Confine the argument of coeffCOS to the interval [-1;1]!
     29                        phi = acos(min(1.0,max(-1.0,Q/sqrt(-P))));
     30                        P   = 2.0*pow((-P),(5.e-1*third));
     31                        for(int i=0;i<3;i++)    U[i] = P*cos((phi+2*((IssmDouble)i)*pi)*third)-W;
     32                        x[0] = min(U[0], min(U[1], U[2]));
     33                        x[1] = max(min(U[0], U[1]),max( min(U[0], U[2]), min(U[1], U[2])));
     34                        x[2] = max(U[0], max(U[1], U[2]));
    5135                        *num = 3;
    5236                }
    5337                else{
    5438                        // only one real solution!
    55                         DIS = sqrt(DIS);
    56                         X[0] = CBRT(Q+DIS)+CBRT(Q-DIS)-W;
     39                        delta = sqrt(delta);
     40                        x[0] = pow(Q+delta,third)+pow(Q-delta,third)-W;
    5741                        *num=1;
    5842                }
    5943        }
    60         else if (A[2] != 0.0){
     44        else if (b != 0.0){
    6145                // quadratic problem
    62                 P = 0.5*A[1]/A[2];
    63                 DIS = pow(P,2)-A[0]/A[2];
    64                 if (DIS > 0.0)
    65                 {
     46                P     = 0.5*c/b;
     47                delta = pow(P,2)-d/b;
     48                if (delta > 0.0){
    6649                        // 2 real solutions
    67                         X[0] = -P - sqrt(DIS);
    68                         X[1] = -P + sqrt(DIS);
    69                         *num=2;
     50                        x[0] = -P - sqrt(delta);
     51                        x[1] = -P + sqrt(delta);
     52                        *num = 2;
    7053                }
    71                 else
    72                 {
     54                else{
    7355                        // no real solution
    74                         *num=0;
     56                        *num = 0;
    7557                }
    7658        }
    77         else if (A[1] != 0.0)
    78         {
     59        else if (c != 0.0){
    7960                //linear equation
    80                 X[0] =A[0]/A[1];
    81                 *num=1;
     61                x[0] = d/c;
     62                *num = 1;
    8263        }
    83         else
    84         {
     64        else{
    8565                //no equation
    86                 *num=0;
     66                *num = 0;
    8767        }
    88  /*
    89   *     ==== perform one step of a newton iteration in order to minimize
    90   *          round-off errors ====
    91   */
    92         for (i=0;i<*num;i++){
    93                 X[i] = X[i] - (A[0]+X[i]*(A[1]+X[i]*(A[2]+X[i]*A[3])))/(A[1]+X[i]*(2.0*A[2]+X[i]*3.0*A[3]));
     68
     69        /* perform one step of a newton iteration in order to minimize round-off errors */
     70        for(int i=0;i<*num;i++){
     71                x[i] = x[i] - (d+x[i]*(c+x[i]*(b+x[i]*a)))/(c+x[i]*(2.0*b+x[i]*3.0*a));
    9472        }
     73
    9574        return 0;
    9675}
  • issm/trunk-jpl/src/c/shared/Numerics/numerics.h

    r13050 r13061  
    1717struct OptPars;
    1818
    19 IssmDouble min(IssmDouble a,IssmDouble b);
    20 IssmDouble max(IssmDouble a,IssmDouble b);
    21 int    min(int a,int b);
    22 int    max(int a,int b);
    23 IssmDouble OptFunc(IssmDouble scalar, OptArgs* optargs);
    24 void   BrentSearch(IssmDouble* psearch_scalar,IssmDouble* pJ,OptPars* optpars,IssmDouble (*f)(IssmDouble,OptArgs*), OptArgs* optargs);
    25 void   OptimalSearch(IssmDouble* psearch_scalar,IssmDouble* pJ,OptPars* optpars,IssmDouble (*f)(IssmDouble,OptArgs*), OptArgs* optargs);
    26 void   cross(IssmDouble* result,IssmDouble* vector1,IssmDouble* vector2);
    27 void   IsInputConverged(IssmDouble* peps, Input** new_inputs,Input** old_inputs,int num_inputs,int criterion_enum);
    28 void   UnitConversion(IssmDouble* values, int numvalues,int direction_enum, int type_enum);
    29 IssmDouble UnitConversion(IssmDouble value, int direction_enum, int type_enum);
    30 char*  OptionsFromAnalysis(Parameters* parameters,int analysis_type);
    31 void   XZvectorsToCoordinateSystem(IssmDouble* T,IssmDouble* xzvectors);
    32 int cubic(double A[4], double X[3], int* num);
     19IssmDouble  min(IssmDouble a,IssmDouble b);
     20IssmDouble  max(IssmDouble a,IssmDouble b);
     21int         min(int a,int b);
     22int         max(int a,int b);
     23IssmDouble  OptFunc(IssmDouble scalar, OptArgs *optargs);
     24void        BrentSearch(IssmDouble *psearch_scalar,IssmDouble*pJ,OptPars*optpars,IssmDouble (*f)(IssmDouble,OptArgs*), OptArgs*optargs);
     25void        OptimalSearch(IssmDouble *psearch_scalar,IssmDouble*pJ,OptPars*optpars,IssmDouble (*f)(IssmDouble,OptArgs*), OptArgs*optargs);
     26void        cross(IssmDouble *result,IssmDouble*vector1,IssmDouble*vector2);
     27void        IsInputConverged(IssmDouble *peps, Input**new_inputs,Input**old_inputs,int num_inputs,int criterion_enum);
     28void        UnitConversion(IssmDouble *values, int numvalues,int direction_enum, int type_enum);
     29IssmDouble  UnitConversion(IssmDouble value, int direction_enum, int type_enum);
     30char       *OptionsFromAnalysis(Parameters *parameters,int analysis_type);
     31void        XZvectorsToCoordinateSystem(IssmDouble *T,IssmDouble*xzvectors);
     32int         cubic(IssmDouble a, IssmDouble b, IssmDouble c, IssmDouble d, double X[3], int *num);
    3333#ifdef _HAVE_PETSC_
    3434void   PetscOptionsFromAnalysis(Parameters* parameters,int analysis_type);
Note: See TracChangeset for help on using the changeset viewer.