Changeset 13061

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:

• cubic.cpp (modified) (1 diff)
• numerics.h (modified) (1 diff)

issm/trunk-jpl/src/c/shared/Numerics/cubic.cpp

@@ -2 +2 @@
 #include "./numerics.h"
 
-int signR(double Z){
-        int ret;
-        if (Z > 0.0)    ret = 1;
-        if (Z < 0.0)    ret = -1;
-        if (Z == 0.0)   ret =0;
+int cubic(IssmDouble a,IssmDouble b,IssmDouble c,IssmDouble d, IssmDouble x[3], int* num){
+        /* Find the real roots of linear/quadratic and cubic equations:
+         *
+         *   a x^3 + bx^2 + cx + d = 0
+         *
+         *   returns the roots in x
+         *   num is the number of roots */
 
-        return ret;
-}
+        /*Some useful constants*/
+        const IssmDouble pi    = 3.1415926535897932;
+        const IssmDouble third = 1./3.;
 
-double CBRT(double Z){
-        double ret;
-        const double THIRD = 1./3.;
-        //define cubic root as statement function
-        //SIGN has different meanings in both C and Fortran
-        // Was unable to use the sign command of C, so wrote my own
-        // that why a new variable needs to be introduced that keeps track of the sign of
-        // SIGN is supposed to return a 1, -1 or 0 depending on what the sign of the argument is
-        ret = fabs(pow(fabs(Z),THIRD)) * signR(Z);
-        return ret;
-}
+        /*Intermediaries*/
+        IssmDouble U[3],W,P,Q,delta,phi;
 
-int cubic(double A[4], double X[3], int* num){
-        const double PI = 3.1415926535897932;
-        const double THIRD = 1./3.;
-        double U[3],W, P, Q, DIS, PHI;
-        int i;
-
-        //define cubic root as statement function
-        // In C, the function is defined outside of the cubic fct
-
-        // ====determine the degree of the polynomial ====
-
-        if (A[3] != 0.0){
+        /* determine the degree of the polynomial */
+        if (a != 0.0){
                 //cubic problem
-                W = A[2]/A[3]*THIRD;
-                P = pow((A[1]/A[3]*THIRD - pow(W,2)),3);
-                Q = -.5*(2.0*pow(W,3)-(A[1]*W-A[0])/A[3] );
-                DIS = pow(Q,2)+P;
-                if ( DIS < 0.0 ){
+                W     = b/a *third;
+                P     = pow((c/a *third - pow(W,2)),3);
+                Q     = -.5 *(2.0*pow(W,3)-(c*W-d)/a );
+                delta = pow(Q,2)+P;
+                if ( delta < 0.0 ){
                         //three real solutions!
-                        //Confine the argument of ACOS to the interval [-1;1]!
-                        PHI = acos(min(1.0,max(-1.0,Q/sqrt(-P))));
-                        P=2.0*pow((-P),(5.e-1*THIRD));
-                        for (i=0;i<3;i++)       U[i] = P*cos((PHI+2*((double)i)*PI)*THIRD)-W;
-                        X[0] = min(U[0], min(U[1], U[2]));
-                        X[1] = max(min(U[0], U[1]),max( min(U[0], U[2]), min(U[1], U[2])));
-                        X[2] = max(U[0], max(U[1], U[2]));
+                        //Confine the argument of coeffCOS to the interval [-1;1]!
+                        phi = acos(min(1.0,max(-1.0,Q/sqrt(-P))));
+                        P   = 2.0*pow((-P),(5.e-1*third));
+                        for(int i=0;i<3;i++)    U[i] = P*cos((phi+2*((IssmDouble)i)*pi)*third)-W;
+                        x[0] = min(U[0], min(U[1], U[2]));
+                        x[1] = max(min(U[0], U[1]),max( min(U[0], U[2]), min(U[1], U[2])));
+                        x[2] = max(U[0], max(U[1], U[2]));
                         *num = 3;
                 }
                 else{
                         // only one real solution!
-                        DIS = sqrt(DIS);
-                        X[0] = CBRT(Q+DIS)+CBRT(Q-DIS)-W;
+                        delta = sqrt(delta);
+                        x[0] = pow(Q+delta,third)+pow(Q-delta,third)-W;
                         *num=1;
                 }
         }
-        else if (A[2] != 0.0){
+        else if (b != 0.0){
                 // quadratic problem
-                P = 0.5*A[1]/A[2];
-                DIS = pow(P,2)-A[0]/A[2];
-                if (DIS > 0.0)
-                {
+                P     = 0.5*c/b;
+                delta = pow(P,2)-d/b;
+                if (delta > 0.0){
                         // 2 real solutions
-                        X[0] = -P - sqrt(DIS);
-                        X[1] = -P + sqrt(DIS);
-                        *num=2;
+                        x[0] = -P - sqrt(delta);
+                        x[1] = -P + sqrt(delta);
+                        *num = 2;
                 }
-                else
-                {
+                else{
                         // no real solution
-                        *num=0;
+                        *num = 0;
                 }
         }
-        else if (A[1] != 0.0)
-        {
+        else if (c != 0.0){
                 //linear equation
-                X[0] =A[0]/A[1];
-                *num=1;
+                x[0] = d/c;
+                *num = 1;
         }
-        else
-        {
+        else{
                 //no equation
-                *num=0;
+                *num = 0;
         }
- /*
-  *     ==== perform one step of a newton iteration in order to minimize
-  *          round-off errors ====
-  */
-        for (i=0;i<*num;i++){
-                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]));
+
+        /* perform one step of a newton iteration in order to minimize round-off errors */
+        for(int i=0;i<*num;i++){
+                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));
         }
+
         return 0;
 }

issm/trunk-jpl/src/c/shared/Numerics/numerics.h

@@ -17 +17 @@
 struct OptPars;
 
-IssmDouble min(IssmDouble a,IssmDouble b);
-IssmDouble max(IssmDouble a,IssmDouble b);
-int    min(int a,int b);
-int    max(int a,int b);
-IssmDouble OptFunc(IssmDouble scalar, OptArgs* optargs);
-void   BrentSearch(IssmDouble* psearch_scalar,IssmDouble* pJ,OptPars* optpars,IssmDouble (*f)(IssmDouble,OptArgs*), OptArgs* optargs);
-void   OptimalSearch(IssmDouble* psearch_scalar,IssmDouble* pJ,OptPars* optpars,IssmDouble (*f)(IssmDouble,OptArgs*), OptArgs* optargs);
-void   cross(IssmDouble* result,IssmDouble* vector1,IssmDouble* vector2);
-void   IsInputConverged(IssmDouble* peps, Input** new_inputs,Input** old_inputs,int num_inputs,int criterion_enum);
-void   UnitConversion(IssmDouble* values, int numvalues,int direction_enum, int type_enum);
-IssmDouble UnitConversion(IssmDouble value, int direction_enum, int type_enum);
-char*  OptionsFromAnalysis(Parameters* parameters,int analysis_type);
-void   XZvectorsToCoordinateSystem(IssmDouble* T,IssmDouble* xzvectors);
-int cubic(double A[4], double X[3], int* num);
+IssmDouble  min(IssmDouble a,IssmDouble b);
+IssmDouble  max(IssmDouble a,IssmDouble b);
+int         min(int a,int b);
+int         max(int a,int b);
+IssmDouble  OptFunc(IssmDouble scalar, OptArgs *optargs);
+void        BrentSearch(IssmDouble *psearch_scalar,IssmDouble*pJ,OptPars*optpars,IssmDouble (*f)(IssmDouble,OptArgs*), OptArgs*optargs);
+void        OptimalSearch(IssmDouble *psearch_scalar,IssmDouble*pJ,OptPars*optpars,IssmDouble (*f)(IssmDouble,OptArgs*), OptArgs*optargs);
+void        cross(IssmDouble *result,IssmDouble*vector1,IssmDouble*vector2);
+void        IsInputConverged(IssmDouble *peps, Input**new_inputs,Input**old_inputs,int num_inputs,int criterion_enum);
+void        UnitConversion(IssmDouble *values, int numvalues,int direction_enum, int type_enum);
+IssmDouble  UnitConversion(IssmDouble value, int direction_enum, int type_enum);
+char       *OptionsFromAnalysis(Parameters *parameters,int analysis_type);
+void        XZvectorsToCoordinateSystem(IssmDouble *T,IssmDouble*xzvectors);
+int         cubic(IssmDouble a, IssmDouble b, IssmDouble c, IssmDouble d, double X[3], int *num);
 #ifdef _HAVE_PETSC_
 void   PetscOptionsFromAnalysis(Parameters* parameters,int analysis_type);