Index: /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/c/shared/Numerics/cubic.cpp =================================================================== --- /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/c/shared/Numerics/cubic.cpp (revision 13060) +++ /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/c/shared/Numerics/cubic.cpp (revision 13061) @@ -1,96 +1,75 @@ #include #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; } Index: /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/c/shared/Numerics/numerics.h =================================================================== --- /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/c/shared/Numerics/numerics.h (revision 13060) +++ /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/c/shared/Numerics/numerics.h (revision 13061) @@ -16,20 +16,20 @@ struct OptArgs; 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); #endif