C++ N_VSetArrayPointer函数代码示例

/*-----------------------------------------------------------------
  cvDlsDenseDQJac 
  -----------------------------------------------------------------
  This routine generates a dense difference quotient approximation 
  to the Jacobian of f(t,y). It assumes that a dense SUNMatrix is 
  stored column-wise, and that elements within each column are 
  contiguous. The address of the jth column of J is obtained via
  the accessor function SUNDenseMatrix_Column, and this pointer 
  is associated with an N_Vector using the N_VSetArrayPointer
  function.  Finally, the actual computation of the jth column of 
  the Jacobian is done with a call to N_VLinearSum.
  -----------------------------------------------------------------*/ 
int cvDlsDenseDQJac(realtype t, N_Vector y, N_Vector fy, 
                    SUNMatrix Jac, CVodeMem cv_mem, N_Vector tmp1)
{
  realtype fnorm, minInc, inc, inc_inv, yjsaved, srur;
  realtype *y_data, *ewt_data;
  N_Vector ftemp, jthCol;
  sunindextype j, N;
  int retval = 0;
  CVDlsMem cvdls_mem;

  /* access DlsMem interface structure */
  cvdls_mem = (CVDlsMem) cv_mem->cv_lmem;

  /* access matrix dimension */
  N = SUNDenseMatrix_Rows(Jac);

  /* Rename work vector for readibility */
  ftemp = tmp1;

  /* Create an empty vector for matrix column calculations */
  jthCol = N_VCloneEmpty(tmp1);

  /* Obtain pointers to the data for ewt, y */
  ewt_data = N_VGetArrayPointer(cv_mem->cv_ewt);
  y_data   = N_VGetArrayPointer(y);

  /* Set minimum increment based on uround and norm of f */
  srur = SUNRsqrt(cv_mem->cv_uround);
  fnorm = N_VWrmsNorm(fy, cv_mem->cv_ewt);
  minInc = (fnorm != ZERO) ?
    (MIN_INC_MULT * SUNRabs(cv_mem->cv_h) * cv_mem->cv_uround * N * fnorm) : ONE;

  for (j = 0; j < N; j++) {

    /* Generate the jth col of J(tn,y) */

    N_VSetArrayPointer(SUNDenseMatrix_Column(Jac,j), jthCol);

    yjsaved = y_data[j];
    inc = SUNMAX(srur*SUNRabs(yjsaved), minInc/ewt_data[j]);
    y_data[j] += inc;

    retval = cv_mem->cv_f(t, y, ftemp, cv_mem->cv_user_data);
    cvdls_mem->nfeDQ++;
    if (retval != 0) break;
    
    y_data[j] = yjsaved;

    inc_inv = ONE/inc;
    N_VLinearSum(inc_inv, ftemp, -inc_inv, fy, jthCol);

    /* DENSE_COL(Jac,j) = N_VGetArrayPointer(jthCol);   /\*UNNECESSARY?? *\/ */
  }

  /* Destroy jthCol vector */
  N_VSetArrayPointer(NULL, jthCol);  /* SHOULDN'T BE NEEDED */
  N_VDestroy(jthCol);

  return(retval);
}

static int KINDenseDQJac(long int n, DenseMat J,
                         N_Vector u, N_Vector fu, void *jac_data,
                         N_Vector tmp1, N_Vector tmp2)
{
  realtype inc, inc_inv, ujsaved, ujscale, sign;
  realtype *tmp2_data, *u_data, *uscale_data;
  N_Vector ftemp, jthCol;
  long int j;
  int retval;

  KINMem kin_mem;
  KINDenseMem  kindense_mem;

  /* jac_data points to kin_mem */
  kin_mem = (KINMem) jac_data;
  kindense_mem = (KINDenseMem) lmem;

  /* Save pointer to the array in tmp2 */
  tmp2_data = N_VGetArrayPointer(tmp2);

  /* Rename work vectors for readibility */
  ftemp = tmp1; 
  jthCol = tmp2;

  /* Obtain pointers to the data for u and uscale */
  u_data   = N_VGetArrayPointer(u);
  uscale_data = N_VGetArrayPointer(uscale);

  /* This is the only for loop for 0..N-1 in KINSOL */

  for (j = 0; j < n; j++) {

    /* Generate the jth col of J(u) */

    N_VSetArrayPointer(DENSE_COL(J,j), jthCol);

    ujsaved = u_data[j];
    ujscale = ONE/uscale_data[j];
    sign = (ujsaved >= ZERO) ? ONE : -ONE;
    inc = sqrt_relfunc*MAX(ABS(ujsaved), ujscale)*sign;
    u_data[j] += inc;

    retval = func(u, ftemp, f_data);
    if (retval != 0) return(-1); 

    u_data[j] = ujsaved;

    inc_inv = ONE/inc;
    N_VLinearSum(inc_inv, ftemp, -inc_inv, fu, jthCol);

  }

  /* Restore original array pointer in tmp2 */
  N_VSetArrayPointer(tmp2_data, tmp2);

  /* Increment counter nfeD */
  nfeD += n;

  return(0);
}

void FIDA_SETVIN(char key_name[], realtype *vval, int *ier)
{
  N_Vector Vec;

  *ier = 0;

  if (!strncmp(key_name,"ID_VEC",6)) {
    Vec = NULL;
    Vec = N_VCloneEmpty(F2C_IDA_vec);
    if (Vec == NULL) {
      *ier = -1;
      return;
    }
    N_VSetArrayPointer(vval, Vec);
    IDASetId(IDA_idamem, Vec);
    N_VDestroy(Vec);
  } else if (!strncmp(key_name,"CONSTR_VEC",10)) {
    Vec = NULL;
    Vec = N_VCloneEmpty(F2C_IDA_vec);
    if (Vec == NULL) {
      *ier = -1;
      return;
    }
    N_VSetArrayPointer(vval, Vec);
    IDASetConstraints(IDA_idamem, Vec);
    N_VDestroy(Vec);
  } else {
    *ier = -99;
    printf("FIDASETVIN: Unrecognized key.\n\n");
  }

}

static void CVDenseDQJac(long int N, DenseMat J, realtype t,
                         N_Vector y, N_Vector fy, void *jac_data,
                         N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
{
    realtype fnorm, minInc, inc, inc_inv, yjsaved, srur;
    realtype *tmp2_data, *y_data, *ewt_data;
    N_Vector ftemp, jthCol;
    long int j;

    CVodeMem cv_mem;
    CVDenseMem  cvdense_mem;

    /* jac_data points to cvode_mem */
    cv_mem = (CVodeMem) jac_data;
    cvdense_mem = (CVDenseMem) lmem;

    /* Save pointer to the array in tmp2 */
    tmp2_data = N_VGetArrayPointer(tmp2);

    /* Rename work vectors for readibility */
    ftemp = tmp1;
    jthCol = tmp2;

    /* Obtain pointers to the data for ewt, y */
    ewt_data = N_VGetArrayPointer(ewt);
    y_data   = N_VGetArrayPointer(y);

    /* 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;

    /* This is the only for loop for 0..N-1 in CVODE */

    for (j = 0; j < N; j++) {

        /* Generate the jth col of J(tn,y) */

        N_VSetArrayPointer(DENSE_COL(J,j), jthCol);

        yjsaved = y_data[j];
        inc = MAX(srur*ABS(yjsaved), minInc/ewt_data[j]);
        y_data[j] += inc;
        f(tn, y, ftemp, f_data);
        y_data[j] = yjsaved;

        inc_inv = ONE/inc;
        N_VLinearSum(inc_inv, ftemp, -inc_inv, fy, jthCol);

        DENSE_COL(J,j) = N_VGetArrayPointer(jthCol);
    }

    /* Restore original array pointer in tmp2 */
    N_VSetArrayPointer(tmp2_data, tmp2);

    /* Increment counter nfeD */
    nfeD += N;
}

void FIDA_REINIT(realtype *t0, realtype *yy0, realtype *yp0,
                 int *iatol, realtype *rtol, realtype *atol,
                 int *ier)
{
  N_Vector Vatol;

  *ier = 0;

  /* Initialize all pointers to NULL */
  Vatol = NULL;

  /* Attach user's yy0 to F2C_IDA_vec */
  N_VSetArrayPointer(yy0, F2C_IDA_vec);

  /* Attach user's yp0 to F2C_IDA_ypvec */
  N_VSetArrayPointer(yp0, F2C_IDA_ypvec);

  /* Call IDAReInit */
  *ier = IDAReInit(IDA_idamem, *t0, F2C_IDA_vec, F2C_IDA_ypvec);

  /* Reset data pointers */
  N_VSetArrayPointer(NULL, F2C_IDA_vec);
  N_VSetArrayPointer(NULL, F2C_IDA_ypvec);

  /* On failure, exit */
  if (*ier != IDA_SUCCESS) {
    *ier = -1;
    return;
  }

  /* Set tolerances */
  switch (*iatol) {
  case 1:
    *ier = IDASStolerances(IDA_idamem, *rtol, *atol);
    break;
  case 2:
    Vatol = NULL;
    Vatol= N_VCloneEmpty(F2C_IDA_vec);
    if (Vatol == NULL) {
      *ier = -1;
      return;
    }
    N_VSetArrayPointer(atol, Vatol);
    *ier = IDASVtolerances(IDA_idamem, *rtol, Vatol);
    N_VDestroy(Vatol);
    break;
  }

  /* On failure, exit */
  if (*ier != IDA_SUCCESS) {
    *ier = -1;
    return;
  }

  return;
}

void FCV_REINIT(realtype *t0, realtype *y0, 
                int *iatol, realtype *rtol, realtype *atol, 
                int *ier)
{
  N_Vector Vatol;

  *ier = 0;

  /* Initialize all pointers to NULL */
  Vatol = NULL;

  /* Set data in F2C_CVODE_vec to y0 */
  N_VSetArrayPointer(y0, F2C_CVODE_vec);

  /* Call CVReInit */
  *ier = CVodeReInit(CV_cvodemem, *t0, F2C_CVODE_vec);

  /* Reset data pointers */
  N_VSetArrayPointer(NULL, F2C_CVODE_vec);

  /* On failure, exit */
  if (*ier != CV_SUCCESS) {
    *ier = -1;
    return;
  }

  /* Set tolerances */
  switch (*iatol) {
  case 1:
    *ier = CVodeSStolerances(CV_cvodemem, *rtol, *atol); 
    break;
  case 2:
    Vatol = NULL;
    Vatol = N_VCloneEmpty(F2C_CVODE_vec);
    if (Vatol == NULL) {
      *ier = -1;
      return;
    }
    N_VSetArrayPointer(atol, Vatol);
    *ier = CVodeSVtolerances(CV_cvodemem, *rtol, Vatol);
    N_VDestroy(Vatol);
    break;
  }

  /* On failure, exit */
  if (*ier != CV_SUCCESS) {
    *ier = -1;
    return;
  }

  return;
}

void FIDA_GETESTLOCALERR(realtype *ele, int *ier)
{
  /* Attach user data to vector */
  N_VSetArrayPointer(ele, F2C_IDA_vec);

  *ier = 0;
  *ier = IDAGetEstLocalErrors(IDA_idamem, F2C_IDA_vec);

  /* Reset data pointers */
  N_VSetArrayPointer(NULL, F2C_IDA_vec);

  return;
}

void FIDA_GETERRWEIGHTS(realtype *eweight, int *ier)
{
  /* Attach user data to vector */
  N_VSetArrayPointer(eweight, F2C_IDA_vec);

  *ier = 0;
  *ier = IDAGetErrWeights(IDA_idamem, F2C_IDA_vec);

  /* Reset data pointer */
  N_VSetArrayPointer(NULL, F2C_IDA_vec);

  return;
}

void FIDA_GETDKY(realtype *t, int *k, realtype *dky, int *ier)
{
  /* Attach user data to vectors */
  N_VSetArrayPointer(dky, F2C_IDA_vec);

  *ier = 0;
  *ier = IDAGetDky(IDA_idamem, *t, *k, F2C_IDA_vec);

  /* Reset data pointers */
  N_VSetArrayPointer(NULL, F2C_IDA_vec);

  return;
}

void FCV_GETERRWEIGHTS(realtype *eweight, int *ier)
{
  /* Attach user data to vector */
  realtype *f2c_data = N_VGetArrayPointer(F2C_CVODE_vec);
  N_VSetArrayPointer(eweight, F2C_CVODE_vec);

  *ier = 0;
  *ier = CVodeGetErrWeights(CV_cvodemem, F2C_CVODE_vec);

  /* Reset data pointers */
  N_VSetArrayPointer(f2c_data, F2C_CVODE_vec);

  return;
}

void FCV_GETESTLOCALERR(realtype *ele, int *ier)
{
  /* Attach user data to vector */
  realtype *f2c_data = N_VGetArrayPointer(F2C_CVODE_vec);
  N_VSetArrayPointer(ele, F2C_CVODE_vec);

  *ier = 0;
  *ier = CVodeGetEstLocalErrors(CV_cvodemem, F2C_CVODE_vec);

  /* Reset data pointers */
  N_VSetArrayPointer(f2c_data, F2C_CVODE_vec);

  return;
}

void FIDA_GETSOL(realtype *t, realtype *yret, realtype *ypret, int *ier)
{
  /* Attach user data to vectors */
  N_VSetArrayPointer(yret, F2C_IDA_vec);
  N_VSetArrayPointer(ypret, F2C_IDA_ypvec);

  *ier = 0;
  *ier = IDAGetSolution(IDA_idamem, *t, F2C_IDA_vec, F2C_IDA_ypvec);

  /* Reset data pointers */
  N_VSetArrayPointer(NULL, F2C_IDA_vec);
  N_VSetArrayPointer(NULL, F2C_IDA_ypvec);

  return;
}

void FIDA_GETDKY(realtype *t, int *k, realtype *dky, int *ier)
{
  /* Store existing F2C_IDA_vec data pointer */
  realtype *f2c_data = N_VGetArrayPointer(F2C_IDA_vec);

  /* Attach user data to vectors */
  N_VSetArrayPointer(dky, F2C_IDA_vec);

  *ier = 0;
  *ier = IDAGetDky(IDA_idamem, *t, *k, F2C_IDA_vec);

  /* Reset data pointers */
  N_VSetArrayPointer(f2c_data, F2C_IDA_vec);

  return;
}

void FIDA_GETESTLOCALERR(realtype *ele, int *ier)
{
  /* Store existing F2C_IDA_vec data pointer */
  realtype *f2c_data = N_VGetArrayPointer(F2C_IDA_vec);

  /* Attach user data to vector */
  N_VSetArrayPointer(ele, F2C_IDA_vec);

  *ier = 0;
  *ier = IDAGetEstLocalErrors(IDA_idamem, F2C_IDA_vec);

  /* Reset data pointers */
  N_VSetArrayPointer(f2c_data, F2C_IDA_vec);

  return;
}

/* Fortran interface routine to set residual tolerance 
   scalar/array; functions as an all-in-one interface to the C 
   routines ARKodeResStolerance and ARKodeResVtolerance; 
   see farkode.h for further details */
void FARK_SETRESTOLERANCE(int *itol, realtype *atol, int *ier) {

  N_Vector Vatol;
  realtype abstol;

  *ier = 0;

  /* Set tolerance, based on itol argument */
  abstol=1.e-9;
  switch (*itol) {
  case 1:
    if (*atol > 0.0)  abstol = *atol;
    *ier = ARKodeResStolerance(ARK_arkodemem, abstol); 
    break;
  case 2:
    Vatol = NULL;
    Vatol = N_VCloneEmpty(F2C_ARKODE_vec);
    if (Vatol == NULL) {
      *ier = -1;
      return;
    }
    N_VSetArrayPointer(atol, Vatol);
    if (N_VMin(Vatol) <= 0.0)  N_VConst(abstol, Vatol);
    *ier = ARKodeResVtolerance(ARK_arkodemem, Vatol);
    N_VDestroy(Vatol);
    break;
  }

  return;
}