text: /*
   * -----------------------------------------------------------------
   * $Revision: 4272 $
   * $Date: 2014-12-02 11:19:41 -0800 (Tue, 02 Dec 2014) $
   * ----------------------------------------------------------------- 
   * Programmer(s): Radu Serban @ LLNL
   * -----------------------------------------------------------------
   * LLNS Copyright Start
   * Copyright (c) 2014, Lawrence Livermore National Security
   * This work was performed under the auspices of the U.S. Department 
   * of Energy by Lawrence Livermore National Laboratory in part under 
   * Contract W-7405-Eng-48 and in part under Contract DE-AC52-07NA27344.
   * Produced at the Lawrence Livermore National Laboratory.
   * All rights reserved.
   * For details, see the LICENSE file.
   * LLNS Copyright End
   * -----------------------------------------------------------------
   * This is the implementation file for the CVSBAND linear solver.
   * -----------------------------------------------------------------
   */
  
  #include <stdio.h>
  #include <stdlib.h>
  
  #include <cvodes/cvodes_band.h>
  #include "cvodes_direct_impl.h"
  #include "cvodes_impl.h"
  
  #include <sundials/sundials_math.h>
  
  /* Constants */
  
  #define ZERO         RCONST(0.0)
  #define ONE          RCONST(1.0)
  #define TWO          RCONST(2.0)
  
  /* CVSBAND linit, lsetup, lsolve, and lfree routines */
  static int cvBandInit(CVodeMem cv_mem);
  static int cvBandSetup(CVodeMem cv_mem, int convfail, N_Vector ypred,
                         N_Vector fpred, booleantype *jcurPtr, N_Vector vtemp1,
                         N_Vector vtemp2, N_Vector vtemp3);
  static int cvBandSolve(CVodeMem cv_mem, N_Vector b, N_Vector weight,
                         N_Vector ycur, N_Vector fcur);
  static void cvBandFree(CVodeMem cv_mem);
  
  /* CVSBAND lfreeB function */
  static void cvBandFreeB(CVodeBMem cvB_mem);
  
  /* 
   * ================================================================
   *
   *                   PART I - forward problems
   *
   * ================================================================
   */
  
  
  /* Readability Replacements */
  
  #define lmm       (cv_mem->cv_lmm)
  #define f         (cv_mem->cv_f)
  #define nst       (cv_mem->cv_nst)
  #define tn        (cv_mem->cv_tn)
  #define h         (cv_mem->cv_h)
  #define gamma     (cv_mem->cv_gamma)
  #define gammap    (cv_mem->cv_gammap)
  #define gamrat    (cv_mem->cv_gamrat)
  #define ewt       (cv_mem->cv_ewt)
  #define nfe       (cv_mem->cv_nfe)
  #define linit     (cv_mem->cv_linit)
  #define lsetup    (cv_mem->cv_lsetup)
  #define lsolve    (cv_mem->cv_lsolve)
  #define lfree     (cv_mem->cv_lfree)
  #define lmem      (cv_mem->cv_lmem)
  #define vec_tmpl     (cv_mem->cv_tempv)
  #define setupNonNull (cv_mem->cv_setupNonNull)
  
  #define mtype      (cvdls_mem->d_type)
  #define n          (cvdls_mem->d_n)
  #define jacDQ      (cvdls_mem->d_jacDQ)
  #define jac        (cvdls_mem->d_bjac)
  #define M          (cvdls_mem->d_M)
  #define mu         (cvdls_mem->d_mu)
  #define ml         (cvdls_mem->d_ml)
  #define smu        (cvdls_mem->d_smu)
  #define lpivots    (cvdls_mem->d_lpivots)
  #define savedJ     (cvdls_mem->d_savedJ)
  #define nstlj      (cvdls_mem->d_nstlj)
  #define nje        (cvdls_mem->d_nje)
  #define nfeDQ      (cvdls_mem->d_nfeDQ)
  #define J_data     (cvdls_mem->d_J_data)
  #define last_flag  (cvdls_mem->d_last_flag)
  
  /*
   * -----------------------------------------------------------------
   * CVBand
   * -----------------------------------------------------------------
   * This routine initializes the memory record and sets various function
   * fields specific to the band linear solver module.  CVBand first calls
   * the existing lfree routine if this is not NULL.  It then sets the
   * cv_linit, cv_lsetup, cv_lsolve, and cv_lfree fields in (*cvode_mem)
   * to be cvBandInit, cvBandSetup, cvBandSolve, and cvBandFree,
   * respectively.  It allocates memory for a structure of type
   * CVDlsMemRec and sets the cv_lmem field in (*cvode_mem) to the
   * address of this structure.  It sets setupNonNull in (*cvode_mem) to be
   * TRUE, b_mu to be mupper, b_ml to be mlower, and the b_jac field to be 
   * CVBandDQJac.
   * Finally, it allocates memory for M, savedJ, and lpivots.  The CVBand
   * return value is SUCCESS = 0, LMEM_FAIL = -1, or LIN_ILL_INPUT = -2.
   *
   * NOTE: The band linear solver assumes a serial implementation
   *       of the NVECTOR package. Therefore, CVBand will first 
   *       test for compatible a compatible N_Vector internal
   *       representation by checking that the function 
   *       N_VGetArrayPointer exists.
   * -----------------------------------------------------------------
   */
                    
  int CVBand(void *cvode_mem, long int N, long int mupper, long int mlower)
  {
    CVodeMem cv_mem;
    CVDlsMem cvdls_mem;
  
    /* Return immediately if cvode_mem is NULL */
    if (cvode_mem == NULL) {
      cvProcessError(NULL, CVDLS_MEM_NULL, "CVSBAND", "CVBand", MSGD_CVMEM_NULL);
      return(CVDLS_MEM_NULL);
    }
    cv_mem = (CVodeMem) cvode_mem;
  
    /* Test if the NVECTOR package is compatible with the BAND solver */
    if (vec_tmpl->ops->nvgetarraypointer == NULL) {
      cvProcessError(cv_mem, CVDLS_ILL_INPUT, "CVSBAND", "CVBand", MSGD_BAD_NVECTOR);
      return(CVDLS_ILL_INPUT);
    }
  
    if (lfree != NULL) lfree(cv_mem);
  
    /* Set four main function fields in cv_mem */  
    linit  = cvBandInit;
    lsetup = cvBandSetup;
    lsolve = cvBandSolve;
    lfree  = cvBandFree;
    
    /* Get memory for CVDlsMemRec */
    cvdls_mem = NULL;
    cvdls_mem = (CVDlsMem) malloc(sizeof(struct CVDlsMemRec));
    if (cvdls_mem == NULL) {
      cvProcessError(cv_mem, CVDLS_MEM_FAIL, "CVSBAND", "CVBand", MSGD_MEM_FAIL);
      return(CVDLS_MEM_FAIL);
    }
  
    /* Set matrix type */
    mtype = SUNDIALS_BAND;
  
    /* Initialize Jacobian-related data */
    jacDQ = TRUE;
    jac = NULL;
    J_data = NULL;
  
    last_flag = CVDLS_SUCCESS;
  
    setupNonNull = TRUE;
    
    /* Load problem dimension */
    n = N;
  
    /* Load half-bandwiths in cvdls_mem */
    ml = mlower;
    mu = mupper;
  
    /* Test ml and mu for legality */
    if ((ml < 0) || (mu < 0) || (ml >= N) || (mu >= N)) {
      cvProcessError(cv_mem, CVDLS_ILL_INPUT, "CVSBAND", "CVBand", MSGD_BAD_SIZES);
      free(cvdls_mem); cvdls_mem = NULL;
      return(CVDLS_ILL_INPUT);
    }
  
    /* Set extended upper half-bandwith for M (required for pivoting) */
    smu = SUNMIN(N-1, mu + ml);
  
    /* Allocate memory for M, savedJ, and pivot arrays */
    M = NULL;
    M = NewBandMat(N, mu, ml, smu);
    if (M == NULL) {
      cvProcessError(cv_mem, CVDLS_MEM_FAIL, "CVSBAND", "CVBand", MSGD_MEM_FAIL);
      free(cvdls_mem); cvdls_mem = NULL;
      return(CVDLS_MEM_FAIL);
    }
    savedJ = NULL;
    savedJ = NewBandMat(N, mu, ml, mu);
    if (savedJ == NULL) {
      cvProcessError(cv_mem, CVDLS_MEM_FAIL, "CVSBAND", "CVBand", MSGD_MEM_FAIL);
      DestroyMat(M);
      free(cvdls_mem); cvdls_mem = NULL;
      return(CVDLS_MEM_FAIL);
    }
    lpivots = NULL;
    lpivots = NewLintArray(N);
    if (lpivots == NULL) {
      cvProcessError(cv_mem, CVDLS_MEM_FAIL, "CVSBAND", "CVBand", MSGD_MEM_FAIL);
      DestroyMat(M);
      DestroyMat(savedJ);
      free(cvdls_mem); cvdls_mem = NULL;
      return(CVDLS_MEM_FAIL);
    }
  
    /* Attach linear solver memory to integrator memory */
    lmem = cvdls_mem;
  
    return(CVDLS_SUCCESS);
  }
  
  /*
   * -----------------------------------------------------------------
   * cvBandInit
   * -----------------------------------------------------------------
   * This routine does remaining initializations specific to the band
   * linear solver.
   * -----------------------------------------------------------------
   */
  
  static int cvBandInit(CVodeMem cv_mem)
  {
    CVDlsMem cvdls_mem;
  
    cvdls_mem = (CVDlsMem) lmem;
  
    nje   = 0;
    nfeDQ = 0;
    nstlj = 0;
  
    /* Set Jacobian function and data, depending on jacDQ */
    if (jacDQ) {
      jac = cvDlsBandDQJac;
      J_data = cv_mem;
    } else {
      J_data = cv_mem->cv_user_data;
    }
  
    last_flag = CVDLS_SUCCESS;
    return(0);
  }
  
  /*
   * -----------------------------------------------------------------
   * cvBandSetup
   * -----------------------------------------------------------------
   * This routine does the setup operations for the band linear solver.
   * It makes a decision whether or not to call the Jacobian evaluation
   * routine based on various state variables, and if not it uses the 
   * saved copy.  In any case, it constructs the Newton matrix 
   * M = I - gamma*J, updates counters, and calls the band LU 
   * factorization routine.
   * -----------------------------------------------------------------
   */
  
  static int cvBandSetup(CVodeMem cv_mem, int convfail, N_Vector ypred,
                         N_Vector fpred, booleantype *jcurPtr, N_Vector vtemp1,
                         N_Vector vtemp2, N_Vector vtemp3)
  {
    CVDlsMem cvdls_mem;
    booleantype jbad, jok;
    realtype dgamma;
    int retval;
    long int ier;
  
    cvdls_mem = (CVDlsMem) lmem;
  
    /* Use nst, gamma/gammap, and convfail to set J eval. flag jok */
  
    dgamma = SUNRabs((gamma/gammap) - ONE);
    jbad = (nst == 0) || (nst > nstlj + CVD_MSBJ) ||
           ((convfail == CV_FAIL_BAD_J) && (dgamma < CVD_DGMAX)) ||
           (convfail == CV_FAIL_OTHER);
    jok = !jbad;
    
    if (jok) {
  
      /* If jok = TRUE, use saved copy of J */
      *jcurPtr = FALSE;
      BandCopy(savedJ, M, mu, ml);
  
    } else {
  
      /* If jok = FALSE, call jac routine for new J value */
      nje++;
      nstlj = nst;
      *jcurPtr = TRUE;
      SetToZero(M); 
  
      retval = jac(n, mu, ml, tn, ypred, fpred, M, J_data, vtemp1, vtemp2, vtemp3);
      if (retval < 0) {
        cvProcessError(cv_mem, CVDLS_JACFUNC_UNRECVR, "CVSBAND", "cvBandSetup", MSGD_JACFUNC_FAILED);
        last_flag = CVDLS_JACFUNC_UNRECVR;
        return(-1);
      }
      if (retval > 0) {
        last_flag = CVDLS_JACFUNC_RECVR;
        return(1);
      }
  
      BandCopy(M, savedJ, mu, ml);
  
    }
    
    /* Scale and add I to get M = I - gamma*J */
    BandScale(-gamma, M);
    AddIdentity(M);
  
    /* Do LU factorization of M */
    ier = BandGBTRF(M, lpivots);
  
    /* Return 0 if the LU was complete; otherwise return 1 */
    if (ier > 0) {
      last_flag = ier;
      return(1);
    }
    last_flag = CVDLS_SUCCESS;
    return(0);
  }
  
  /*
   * -----------------------------------------------------------------
   * cvBandSolve
   * -----------------------------------------------------------------
   * This routine handles the solve operation for the band linear solver
   * by calling the band backsolve routine.  The return value is 0.
   * -----------------------------------------------------------------
   */
  
  static int cvBandSolve(CVodeMem cv_mem, N_Vector b, N_Vector weight,
                         N_Vector ycur, N_Vector fcur)
  {
    CVDlsMem cvdls_mem;
    realtype *bd;
  
    cvdls_mem = (CVDlsMem) lmem;
  
    bd = N_VGetArrayPointer(b);
  
    BandGBTRS(M, lpivots, bd);
  
    /* If CV_BDF, scale the correction to account for change in gamma */
    if ((lmm == CV_BDF) && (gamrat != ONE)) {
      N_VScale(TWO/(ONE + gamrat), b, b);
    }
  
    last_flag = CVDLS_SUCCESS;
    return(0);
  }
  
  /*
   * -----------------------------------------------------------------
   * cvBandFree
   * -----------------------------------------------------------------
   * This routine frees memory specific to the band linear solver.
   * -----------------------------------------------------------------
   */
  
  static void cvBandFree(CVodeMem cv_mem)
  {
    CVDlsMem cvdls_mem;
  
    cvdls_mem = (CVDlsMem) lmem;
  
    DestroyMat(M);
    DestroyMat(savedJ);
    DestroyArray(lpivots);
    free(cvdls_mem);
    cv_mem->cv_lmem = NULL;
  }
  
  /* 
   * ================================================================
   *
   *                   PART II - backward problems
   *
   * ================================================================
   */
  
  /*
   * CVBandB is a wrapper around CVBand. It attaches the CVSBAND linear solver
   * to the backward problem memory block.
   */
  
  int CVBandB(void *cvode_mem, int which,
              long int nB, long int mupperB, long int mlowerB)
  {
    CVodeMem cv_mem;
    CVadjMem ca_mem;
    CVodeBMem cvB_mem;
    void *cvodeB_mem;
    CVDlsMemB cvdlsB_mem;
    int flag;
  
    /* Check if cvode_mem exists */
    if (cvode_mem == NULL) {
      cvProcessError(NULL, CVDLS_MEM_NULL, "CVSBAND", "CVBandB", MSGD_CVMEM_NULL);
      return(CVDLS_MEM_NULL);
    }
    cv_mem = (CVodeMem) cvode_mem;
  
    /* Was ASA initialized? */
    if (cv_mem->cv_adjMallocDone == FALSE) {
      cvProcessError(cv_mem, CVDLS_NO_ADJ, "CVSBAND", "CVBandB", MSGD_NO_ADJ);
      return(CVDLS_NO_ADJ);
    } 
    ca_mem = cv_mem->cv_adj_mem;
  
    /* Check which */
    if ( which >= ca_mem->ca_nbckpbs ) {
      cvProcessError(cv_mem, CVDLS_ILL_INPUT, "CVSBAND", "CVBandB", MSGCV_BAD_WHICH);
      return(CVDLS_ILL_INPUT);
    }
  
    /* Find the CVodeBMem entry in the linked list corresponding to which */
    cvB_mem = ca_mem->cvB_mem;
    while (cvB_mem != NULL) {
      if ( which == cvB_mem->cv_index ) break;
      cvB_mem = cvB_mem->cv_next;
    }
  
    cvodeB_mem = (void *) (cvB_mem->cv_mem);
  
    /* Get memory for CVDlsMemRecB */
    cvdlsB_mem = (CVDlsMemB) malloc(sizeof(struct CVDlsMemRecB));
    if (cvdlsB_mem == NULL) {
      cvProcessError(cv_mem, CVDLS_MEM_FAIL, "CVSBAND", "CVBandB", MSGD_MEM_FAIL);
      return(CVDLS_MEM_FAIL);
    }
  
    /* set matrix type */
    cvdlsB_mem->d_typeB = SUNDIALS_BAND;
  
    /* initialize Jacobian function */
    cvdlsB_mem->d_bjacB = NULL;
  
    /* attach lmemB and lfreeB */
    cvB_mem->cv_lmem = cvdlsB_mem;
    cvB_mem->cv_lfree = cvBandFreeB;
  
    flag = CVBand(cvodeB_mem, nB, mupperB, mlowerB);
  
    if (flag != CVDLS_SUCCESS) {
      free(cvdlsB_mem);
      cvdlsB_mem = NULL;
    }
  
    return(flag);
  }
  
  /*
   * cvBandFreeB frees the memory associated with the CVSBAND linear
   * solver for backward integration.
   */
  
  static void cvBandFreeB(CVodeBMem cvB_mem)
  {
    CVDlsMemB cvdlsB_mem;
  
    cvdlsB_mem = (CVDlsMemB) (cvB_mem->cv_lmem);
  
    free(cvdlsB_mem);
  }
  
  