/* Interface to CVODE for solution of an S-d ODE.
Last modified on 09/21/01 by Sergej V. Aksenov. */

/* Return codes: 
0 - ok, 
1-8 - CVODE error, 
9 - CVODE produces negative solution. */

#include <math.h>
#include <stdio.h>
#include "llnltyps.h" 
#include "llnlmath.h"
#include "cvode.h"    
#include "cvdense.h"  
#include "nvector.h"  
#include "dense.h"   

/* S-d parameters. */

typedef struct
    {
      	double g, h, alpha;
    } *UserData;

#define Ith(v,i) N_VIth(v,i-1)         
#define IJth(A,i,j) DENSE_ELEM(A,i-1,j-1)

/* Rhs of an S-d ODE. */

static void f(integer N,real t,N_Vector y,N_Vector ydot,void *f_data)
    {
      	double y1, g, h, alpha;
      	UserData data;

      	data = (UserData) f_data;
      	g = data->g;
      	h = data->h;
      	alpha = data->alpha;

      	y1 = Ith(y,1);
      	Ith(ydot,1) = alpha*pow(y1,g)-alpha*pow(y1,h);
    }
   
/* Jacobian for the rhs of an S-d ODE. */

static void Jac(integer N,DenseMat J,RhsFn f,void *f_data,real t,N_Vector y,N_Vector fy,N_Vector ewt,real h,real uround,void *jac_data,long int *nfePtr,N_Vector vtemp1,N_Vector vtemp2,N_Vector vtemp3)
    {
      	double y1, g, h1, alpha;
      	UserData data;

      	data = (UserData) f_data;
      	g = data->g;
      	h1 = data->h;
      	alpha = data->alpha;

      	y1 = Ith(y,1);
      	IJth(J,1,1) = alpha*g*pow(y1,g-1.0)-alpha*h1*pow(y1,h1-1.0);
    }

/* Solve an S-d ODE. */

int sdode(double *sdpar, double *f0, double *xend, int *nsolp, double *sol1, double *rtol, double *atol, int *meth, int *steps)
    {
      	real ropt[OPT_SIZE], reltol, abstol, t;
      	long int iopt[OPT_SIZE];
      	N_Vector y;
      	void *cvode_mem;
      	int flag, k, returnstatus, mm, itr, sign;
        CVDenseJacFn jcc;
      	UserData data;
        FILE *myfp;
        double x0, tround, uround, tdist, tdiff;
		
      	/* Initialize CVODE. */
		
        switch(*meth)
            {
                case 1000: mm = BDF; itr = NEWTON; jcc = Jac; break;
                case 1010: mm = BDF; itr = FUNCTIONAL; jcc = Jac; break;
                case 1001: mm = BDF; itr = NEWTON; jcc = NULL; break;
                case 1011: mm = BDF; itr = FUNCTIONAL; jcc = NULL; break;
                case 1100: mm = ADAMS; itr = NEWTON; jcc = Jac; break;
                case 1110: mm = ADAMS; itr = FUNCTIONAL; jcc = Jac; break;
                case 1101: mm = ADAMS; itr = NEWTON; jcc = NULL; break;
                case 1111: mm = ADAMS; itr = FUNCTIONAL; jcc = NULL; break;
                default: mm = BDF; itr = NEWTON; jcc = Jac; break;
            }

      	data = (UserData) malloc(sizeof *data);
      	data->g = sdpar[0];
      	data->h = sdpar[1];
      	data->alpha = sdpar[2];
        x0 = sdpar[3];
      	y = N_VNew(1, NULL);
      	reltol = *rtol;
      	abstol = *atol;
        Ith(y,1) = *f0;
        uround = UnitRoundoff();
        for (k = 0; k <= OPT_SIZE-1; k++)
            {
                iopt[k] = 0;
                ropt[k] = 0.0;
            }
        iopt[MXSTEP] = *steps;
        /* ropt[HMIN] = DBL_EPSILON; */
        /* myfp = fopen("cvode.err","w");
        if (myfp == NULL) myfp = stderr; */
        myfp = stderr;
        
        cvode_mem = CVodeMalloc(1,f,x0,y,mm,itr,SS,&reltol,&abstol,data,myfp,TRUE,iopt,ropt,NULL);
        CVDense(cvode_mem,jcc,data);
		
        /* Compute solution at all points. */

      	for (k = 0; k <= *nsolp; k++)
            {
                tdiff = *(xend+k) - x0;
                sign = (tdiff > 0.0) ? 1 : -1;
                tdist = ABS(tdiff);
                tround = uround * MAX(ABS(x0), ABS(*(xend+k)));
                if ( tdiff == 0.0 || tdist < 2.0*tround) 
                    {
                        *(sol1+k) = *f0;
                        returnstatus = 0;
                    }
  
             /*   if (fabs(*(xend+k)-x0) <= 8*DBL_EPSILON) 
                    {
                        *(sol1+k) = *f0;
                        returnstatus = 0;
                    } */
                else
                    {
                        flag = CVode(cvode_mem,*(xend+k),y,&t,NORMAL);
                        if (flag)
                            {
                                *nsolp = k - 1;
                                returnstatus = flag+9;
                                break;
                            }
                        else
                            {
                                if (Ith(y,1) < 0)
                                    {
                                        *nsolp = k - 1;
                                        returnstatus = 9;
                                        break;
                                    }
                                else
                                    {
                                        *(sol1+k) = Ith(y,1);
                                        returnstatus = 0;
                                    }
                            }
                    }
            }

        /* Clean up. */

        /* fclose(myfp); */
      	CVodeFree(cvode_mem);
      	N_VFree(y);
      	free(data);
      	return returnstatus;
    }

#undef Ith
#undef IJth