ext/cvode/source/cvband.c
void CVBandDQJac(integer N, integer mupper, integer mlower, BandMat J,
		 RhsFn f, void *f_data, real tn, 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)
{
  real    fnorm, minInc, inc, inc_inv, srur;
  N_Vector ftemp, ytemp;
  integer group, i, j, width, ngroups, i1, i2;
  real *col_j, *ewt_data, *fy_data, *ftemp_data, *y_data, *ytemp_data;

  /* Rename work vectors for use as temporary values of y and f */
  ftemp = vtemp1;
  ytemp = vtemp2;

  /* Obtain pointers to the data for ewt, fy, ftemp, y, ytemp */
  ewt_data   = N_VDATA(ewt);
  fy_data    = N_VDATA(fy);
  ftemp_data = N_VDATA(ftemp);
  y_data     = N_VDATA(y);
  ytemp_data = N_VDATA(ytemp);

  /* Load ytemp with y = predicted y vector */
  N_VScale(ONE, y, ytemp);

  /* Set minimum increment based on uround and norm of f */
  srur = RSqrt(uround);
  fnorm = N_VWrmsNorm(fy, ewt);
  minInc = (fnorm != ZERO) ?
           (MIN_INC_MULT * ABS(h) * uround * N * fnorm) : ONE;

  /* Set bandwidth and number of column groups for band differencing */
  width = mlower + mupper + 1;
  ngroups = MIN(width, N);
  
  for (group=1; group <= ngroups; group++) {
    
    /* Increment all y_j in group */
    for(j=group-1; j < N; j+=width) {
      inc = MAX(srur*ABS(y_data[j]), minInc/ewt_data[j]);
      ytemp_data[j] += inc;
    }

    /* Evaluate f with incremented y */
    f(N, tn, ytemp, ftemp, f_data);

    /* Restore ytemp, then form and load difference quotients */
    for (j=group-1; j < N; j+=width) {
      ytemp_data[j] = y_data[j];
      col_j = BAND_COL(J,j);
      inc = MAX(srur*ABS(y_data[j]), minInc/ewt_data[j]);
      inc_inv = ONE/inc;
      i1 = MAX(0, j-mupper);
      i2 = MIN(j+mlower, N-1);
      for (i=i1; i <= i2; i++)
	BAND_COL_ELEM(col_j,i,j) =
	  inc_inv * (ftemp_data[i] - fy_data[i]);
    }
  }
  
  /* Increment counter nfe = *nfePtr */
  *nfePtr += ngroups;
}
/*************** CVBandDQJac *****************************************

 This routine generates a banded difference quotient approximation to
 the Jacobian of f(t,y).  It assumes that a band matrix of type
 BandMat is stored column-wise, and that elements within each column
 are contiguous. This makes it possible to get the address of a column
 of J via the macro BAND_COL and to write a simple for loop to set
 each of the elements of a column in succession.

**********************************************************************/
CVBandDQJac
This routine generates a banded difference quotient approximation to
the Jacobian of f(t,y).  It assumes that a band matrix of type
BandMat is stored column-wise, and that elements within each column
are contiguous. This makes it possible to get the address of a column
of J via the macro BAND_COL and to write a simple for loop to set
each of the elements of a column in succession.
CVBandDQJac This routine generates a banded difference quotient approximation to the Jacobian of f ( t y ) . It assumes that a band matrix of type BandMat is stored column - wise and that elements within each column are contiguous . This makes it possible to get the address of a column of J via the macro BAND_COL and to write a simple for loop to set each of the elements of a column in succession .
void CVBandDQJac(integer N, integer mupper, integer mlower, BandMat J,
		 RhsFn f, void *f_data, real tn, 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)
{
  real    fnorm, minInc, inc, inc_inv, srur;
  N_Vector ftemp, ytemp;
  integer group, i, j, width, ngroups, i1, i2;
  real *col_j, *ewt_data, *fy_data, *ftemp_data, *y_data, *ytemp_data;
  ftemp = vtemp1;
  ytemp = vtemp2;
  ewt_data   = N_VDATA(ewt);
  fy_data    = N_VDATA(fy);
  ftemp_data = N_VDATA(ftemp);
  y_data     = N_VDATA(y);
  ytemp_data = N_VDATA(ytemp);
  N_VScale(ONE, y, ytemp);
  srur = RSqrt(uround);
  fnorm = N_VWrmsNorm(fy, ewt);
  minInc = (fnorm != ZERO) ?
           (MIN_INC_MULT * ABS(h) * uround * N * fnorm) : ONE;
  width = mlower + mupper + 1;
  ngroups = MIN(width, N);
  for (group=1; group <= ngroups; group++) {
    for(j=group-1; j < N; j+=width) {
      inc = MAX(srur*ABS(y_data[j]), minInc/ewt_data[j]);
      ytemp_data[j] += inc;
    }
    f(N, tn, ytemp, ftemp, f_data);
    for (j=group-1; j < N; j+=width) {
      ytemp_data[j] = y_data[j];
      col_j = BAND_COL(J,j);
      inc = MAX(srur*ABS(y_data[j]), minInc/ewt_data[j]);
      inc_inv = ONE/inc;
      i1 = MAX(0, j-mupper);
      i2 = MIN(j+mlower, N-1);
      for (i=i1; i <= i2; i++)
	BAND_COL_ELEM(col_j,i,j) =
	  inc_inv * (ftemp_data[i] - fy_data[i]);
    }
  }
  *nfePtr += ngroups;
}
void CVBandDQJac ( integer N , integer mupper , integer mlower , BandMat J , RhsFn f , void * f_data , real tn , 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 ) { real fnorm , minInc , inc , inc_inv , srur ; N_Vector ftemp , ytemp ; integer group , i , j , width , ngroups , i1 , i2 ; real * col_j , * ewt_data , * fy_data , * ftemp_data , * y_data , * ytemp_data ; ftemp = vtemp1 ; ytemp = vtemp2 ; ewt_data = N_VDATA ( ewt ) ; fy_data = N_VDATA ( fy ) ; ftemp_data = N_VDATA ( ftemp ) ; y_data = N_VDATA ( y ) ; ytemp_data = N_VDATA ( ytemp ) ; N_VScale ( ONE , y , ytemp ) ; srur = RSqrt ( uround ) ; fnorm = N_VWrmsNorm ( fy , ewt ) ; minInc = ( fnorm != ZERO ) ? ( MIN_INC_MULT * ABS ( h ) * uround * N * fnorm ) : ONE ; width = mlower + mupper + 1 ; ngroups = MIN ( width , N ) ; for ( group = 1 ; group <= ngroups ; group ++ ) { for ( j = group - 1 ; j < N ; j += width ) { inc = MAX ( srur * ABS ( y_data [ j ] ) , minInc / ewt_data [ j ] ) ; ytemp_data [ j ] += inc ; } f ( N , tn , ytemp , ftemp , f_data ) ; for ( j = group - 1 ; j < N ; j += width ) { ytemp_data [ j ] = y_data [ j ] ; col_j = BAND_COL ( J , j ) ; inc = MAX ( srur * ABS ( y_data [ j ] ) , minInc / ewt_data [ j ] ) ; inc_inv = ONE / inc ; i1 = MAX ( 0 , j - mupper ) ; i2 = MIN ( j + mlower , N - 1 ) ; for ( i = i1 ; i <= i2 ; i ++ ) BAND_COL_ELEM ( col_j , i , j ) = inc_inv * ( ftemp_data [ i ] - fy_data [ i ] ) ; } } * nfePtr += ngroups ; }
CVBandDQJac 
This routine generates a banded difference quotient approximation to
the Jacobian of f(t,y).
CVBandDQJac This routine generates a banded difference quotient approximation to the Jacobian of f ( t y ) .
/* Rename work vectors for use as temporary values of y and f */ /* Obtain pointers to the data for ewt, fy, ftemp, y, ytemp */ /* Load ytemp with y = predicted y vector */ /* Set minimum increment based on uround and norm of f */ /* Set bandwidth and number of column groups for band differencing */ /* Increment all y_j in group */ /* Evaluate f with incremented y */ /* Restore ytemp, then form and load difference quotients */ /* Increment counter nfe = *nfePtr */