solvers/ida/idalinear.c

#include "idalinear.h"
#include <ida/ida_impl.h>
#include <utilities/error.h>
#include <utilities/ascMalloc.h>

#include <sundials/sundials_math.h>
#define ZERO RCONST(0.0)
#define ONE  RCONST(1.0)
#define TWO  RCONST(2.0)

typedef struct IntegratorIdaAscendMemStruct{
    long                   integ_neq;   /* problem size */
    IntegratorSparseJacFn *integ_jacfn; /* sparse mtx jacobian evaluation function */
    void *                 integ_jac_data;               /* data for use by jacobvian evaluation function */
    int                    integ_lastflag;
    unsigned long          integ_nje;
    unsigned long          integ_nre;
    mtx_matrix_t           integ_sparse_jac_matrix;
} IntegratorIdaAscendMem;

/* readability replacements (see also ida_dense.c from SUNDIALS distro */

#define linit        (IDA_mem->ida_linit)
#define lsetup       (IDA_mem->ida_lsetup)
#define lsolve       (IDA_mem->ida_lsolve)
#define lperf        (IDA_mem->ida_lperf)
#define lfree        (IDA_mem->ida_lfree)
#define lmem         (IDA_mem->ida_lmem)
#define tn           (IDA_mem->ida_tn)
#define cjratio      (IDA_mem->ida_cjratio)
#define cj           (IDA_mem->ida_cj)

#define nje          (iamem->integ_nje)
#define nre          (iamem->integ_nre)
#define lastflag     (iamem->integ_lastflag)
#define neq          (iamem->integ_neq)
#define jacfn        (iamem->integ_jacfn)
#define jacdata      (iamem->integ_jac_data)
#define JJ           (iamem->integ_sparse_jac_matrix)

#define MSGD_IDAMEM_NULL "Integrator memory is NULL."
#define MSGD_MEM_FAIL    "A memory request failed."
#define MSGD_LMEM_NULL   "IDAASCEND memory is NULL."
#define MSGD_JACFN_UNDEF "The sparse jacobian evaluation routine has not been provided."

/*------------------------------------
  Internal setup/evaluation routines required by IDA. As documented in IDA Manual Ch. 8
*/
int integrator_ida_linit(IDAMem ida_mem);

int integrator_ida_lsetup(IDAMem ida_mem, N_Vector yyp, N_Vector ypp,
    N_Vector resp, N_Vector vtemp1, N_Vector vtemp2, N_Vector vtemp3
);

int integrator_ida_lsolve(IDAMem ida_mem, N_Vector b, N_Vector weight,
    N_Vector ycur, N_Vector ypcur, N_Vector rescur
);

int integrator_ida_lfree(IDAMem ida_mem);

/*------------------------------------
  User functions (called from ida.c in this directory)
*/

int IDAASCEND(void *ida_mem, long _neq){
    IDAMem IDA_mem;
    IntegratorIdaAscendMem *iamem;

    if(ida_mem == NULL){
        IDAProcessError(NULL, IDAASCEND_MEM_NULL, "IDAASCEND", __FUNCTION__, MSGD_IDAMEM_NULL);
        return(IDAASCEND_MEM_NULL);
    }
    IDA_mem = (IDAMem)ida_mem;

    iamem = ASC_NEW(IntegratorIdaAscendMem);
    if(iamem == NULL){
        return IDAASCEND_MEM_FAIL;
    }
    lmem = (void *)iamem;

    /* free linsolver memory with the previous lfree fn, if allocated */
    if(lfree != NULL)lfree(ida_mem);

    /* set the internal-use linear solver function pointers for IDA */
    linit  = &integrator_ida_linit;
    lsetup = &integrator_ida_lsetup;
    lsolve = &integrator_ida_lsolve;
    lperf  = NULL;
    lfree  = &integrator_ida_lfree;

    /* no jacobian assigned, initially (we will throw an error if the user doesn't assign it though) */
    jacfn = NULL;
    lastflag = IDAASCEND_SUCCESS;
    neq = _neq;

    /* allocate memory required for Jacobian mtx_matrix_t ?? */

    return IDAASCEND_SUCCESS;
}

int IDAASCENDSetJacFn(void *ida_mem, IntegratorSparseJacFn *_jacfn, void *_jac_data){
    IDAMem IDA_mem;
    IntegratorIdaAscendMem *iamem;

    if(ida_mem == NULL){
        IDAProcessError(NULL, IDAASCEND_MEM_NULL, "IDAASCEND", __FUNCTION__, MSGD_IDAMEM_NULL);
        return(IDAASCEND_MEM_NULL);
    }
    IDA_mem = (IDAMem)ida_mem;

    if(lmem == NULL){
        IDAProcessError(ida_mem, IDAASCEND_LMEM_NULL, "IDAASCEND", __FUNCTION__, MSGD_LMEM_NULL);
        return(IDAASCEND_LMEM_NULL);
    }
    iamem = (IntegratorIdaAscendMem *)lmem;

    jacfn = _jacfn;

    return IDAASCEND_SUCCESS;
}

int IDAASCENDGetLastFlag(void *ida_mem, int *flag){
    IDAMem IDA_mem;
    IntegratorIdaAscendMem *iamem;

    if(ida_mem == NULL){
        IDAProcessError(NULL, IDAASCEND_MEM_NULL, "IDAASCEND", __FUNCTION__, MSGD_IDAMEM_NULL);
        return(IDAASCEND_MEM_NULL);
    }
    IDA_mem = (IDAMem)ida_mem;

    if(lmem == NULL){
        IDAProcessError(ida_mem, IDAASCEND_LMEM_NULL, "IDAASCEND", __FUNCTION__, MSGD_LMEM_NULL);
        return(IDAASCEND_LMEM_NULL);
    }
    iamem = (IntegratorIdaAscendMem *)lmem;

    *flag = lastflag;

    return IDAASCEND_SUCCESS;
}

char *IDAASCENDGetReturnFlagName(int flag){
    char *name;

    name = ASC_NEW_ARRAY(char,30);

#define FLAG(N) case N:sprintf(name,#N);break

    switch(flag) {
        FLAG(IDAASCEND_SUCCESS);
        FLAG(IDAASCEND_MEM_NULL);
        FLAG(IDAASCEND_LMEM_NULL);
        FLAG(IDAASCEND_MEM_FAIL);
        FLAG(IDAASCEND_JACFN_UNDEF);
        FLAG(IDAASCEND_JACFN_UNRECVR);
        FLAG(IDAASCEND_JACFN_RECVR);
        default:
            sprintf(name,"Unknown flag value '%d'",flag);
    }

    return name;
}


/*------------------------------------
  Internal setup/evaluation routines required by IDA
*/

int integrator_ida_linit(IDAMem IDA_mem){
    IntegratorIdaAscendMem *iamem;
    iamem = (IntegratorIdaAscendMem *)lmem;
  
    CONSOLE_DEBUG("Initialising IDA linear solver");
    nje = 0;
    nre = 0;
    jacfn = NULL;

    /* initialise anything else in the IntegratorIdaAscendMem that needs it */

    return 0;
}

int integrator_ida_lsetup(IDAMem IDA_mem
    , N_Vector yyp, N_Vector ypp, N_Vector rrp
    , N_Vector tmp1, N_Vector tmp2, N_Vector tmp3
){
    int retval;

    IntegratorIdaAscendMem *iamem;
    iamem = (IntegratorIdaAscendMem *)lmem;

    CONSOLE_DEBUG("Setting up IDA linear problem");

    if(jacfn==NULL){
        lastflag = IDAASCEND_JACFN_UNDEF;
        return -1; /* unrecoverable */
    }

    /* Increment nje counter. */
    nje++;

    /* clear the jacobian matrix */
    /* ... */

    /* evaluate the jacobian */
    retval = jacfn(neq, tn, yyp, ypp, rrp, cj, jacdata, JJ, 
        tmp1, tmp2, tmp3
    );

    if(retval < 0){
        lastflag = IDAASCEND_JACFN_UNRECVR;
        return -1;
    }
    if (retval > 0) {
        lastflag = IDAASCEND_JACFN_RECVR;
        return +1;
    }

    /* do block decomposition, LU factorisation or whatever, return success or fail flag */
    /* ... */

    ERROR_REPORTER_HERE(ASC_PROG_ERR,"Not implemented");
    return(-1);
}

/**
    This routines handles the linear solve operation for the IDAASCEND linear
    solver. It interfaces to the appropriate mtx_matrix routines for this,
    but also scales solution vector according to cjratio.

    @return IDAASCEND_SUCESS on success
*/
int integrator_ida_lsolve(IDAMem IDA_mem
    , N_Vector b, N_Vector weight
    , N_Vector ycur, N_Vector ypcur, N_Vector rrcur
){
    realtype *bd;
  
    IntegratorIdaAscendMem *iamem;
    iamem = (IntegratorIdaAscendMem *)lmem;
  
    /* retrieve the data array for the RHS vector, 'b' */
    bd = N_VGetArrayPointer(b);

    /* call the necessary linsolqr routine here */
    /* ... */
    /* For IDADENSE, it's: DenseGETRS(JJ, pivots, bd); */

    /* not sure what this is doing yet */
    /* For IDADENSE, it's: Scale the correction to account for change in cj. */
    if(cjratio != ONE){
        N_VScale(TWO/(ONE + cjratio), b, b);
    }

    ERROR_REPORTER_HERE(ASC_PROG_ERR,"Not implemented");
    return -1;
}

int integrator_ida_lfree(IDAMem IDA_mem){
    CONSOLE_DEBUG("Freeing IDA linear solver data");

    /* free jacobian mtx_matrix_t data ?? */

    if(lmem!=NULL){
        ASC_FREE(lmem);
        lmem=NULL;
    }
    return 0;
}
1	#include "idalinear.h"
2	#include <ida/ida_impl.h>
3	#include <utilities/error.h>
4	#include <utilities/ascMalloc.h>
5
6	#include <sundials/sundials_math.h>
7	#define ZERO RCONST(0.0)
8	#define ONE RCONST(1.0)
9	#define TWO RCONST(2.0)
10
11	typedef struct IntegratorIdaAscendMemStruct{
12	long integ_neq; /* problem size */
13	IntegratorSparseJacFn integ_jacfn; / sparse mtx jacobian evaluation function */
14	void * integ_jac_data; /* data for use by jacobvian evaluation function */
15	int integ_lastflag;
16	unsigned long integ_nje;
17	unsigned long integ_nre;
18	mtx_matrix_t integ_sparse_jac_matrix;
19	} IntegratorIdaAscendMem;
20
21	/* readability replacements (see also ida_dense.c from SUNDIALS distro */
22
23	#define linit (IDA_mem->ida_linit)
24	#define lsetup (IDA_mem->ida_lsetup)
25	#define lsolve (IDA_mem->ida_lsolve)
26	#define lperf (IDA_mem->ida_lperf)
27	#define lfree (IDA_mem->ida_lfree)
28	#define lmem (IDA_mem->ida_lmem)
29	#define tn (IDA_mem->ida_tn)
30	#define cjratio (IDA_mem->ida_cjratio)
31	#define cj (IDA_mem->ida_cj)
32
33	#define nje (iamem->integ_nje)
34	#define nre (iamem->integ_nre)
35	#define lastflag (iamem->integ_lastflag)
36	#define neq (iamem->integ_neq)
37	#define jacfn (iamem->integ_jacfn)
38	#define jacdata (iamem->integ_jac_data)
39	#define JJ (iamem->integ_sparse_jac_matrix)
40
41	#define MSGD_IDAMEM_NULL "Integrator memory is NULL."
42	#define MSGD_MEM_FAIL "A memory request failed."
43	#define MSGD_LMEM_NULL "IDAASCEND memory is NULL."
44	#define MSGD_JACFN_UNDEF "The sparse jacobian evaluation routine has not been provided."
45
46	/*------------------------------------
47	Internal setup/evaluation routines required by IDA. As documented in IDA Manual Ch. 8
48	*/
49	int integrator_ida_linit(IDAMem ida_mem);
50
51	int integrator_ida_lsetup(IDAMem ida_mem, N_Vector yyp, N_Vector ypp,
52	N_Vector resp, N_Vector vtemp1, N_Vector vtemp2, N_Vector vtemp3
53	);
54
55	int integrator_ida_lsolve(IDAMem ida_mem, N_Vector b, N_Vector weight,
56	N_Vector ycur, N_Vector ypcur, N_Vector rescur
57	);
58
59	int integrator_ida_lfree(IDAMem ida_mem);
60
61	/*------------------------------------
62	User functions (called from ida.c in this directory)
63	*/
64
65	int IDAASCEND(void *ida_mem, long _neq){
66	IDAMem IDA_mem;
67	IntegratorIdaAscendMem *iamem;
68
69	if(ida_mem == NULL){
70	IDAProcessError(NULL, IDAASCEND_MEM_NULL, "IDAASCEND", __FUNCTION__, MSGD_IDAMEM_NULL);
71	return(IDAASCEND_MEM_NULL);
72	}
73	IDA_mem = (IDAMem)ida_mem;
74
75	iamem = ASC_NEW(IntegratorIdaAscendMem);
76	if(iamem == NULL){
77	return IDAASCEND_MEM_FAIL;
78	}
79	lmem = (void *)iamem;
80
81	/* free linsolver memory with the previous lfree fn, if allocated */
82	if(lfree != NULL)lfree(ida_mem);
83
84	/* set the internal-use linear solver function pointers for IDA */
85	linit = &integrator_ida_linit;
86	lsetup = &integrator_ida_lsetup;
87	lsolve = &integrator_ida_lsolve;
88	lperf = NULL;
89	lfree = &integrator_ida_lfree;
90
91	/* no jacobian assigned, initially (we will throw an error if the user doesn't assign it though) */
92	jacfn = NULL;
93	lastflag = IDAASCEND_SUCCESS;
94	neq = _neq;
95
96	/* allocate memory required for Jacobian mtx_matrix_t ?? */
97
98	return IDAASCEND_SUCCESS;
99	}
100
101	int IDAASCENDSetJacFn(void ida_mem, IntegratorSparseJacFn _jacfn, void *_jac_data){
102	IDAMem IDA_mem;
103	IntegratorIdaAscendMem *iamem;
104
105	if(ida_mem == NULL){
106	IDAProcessError(NULL, IDAASCEND_MEM_NULL, "IDAASCEND", __FUNCTION__, MSGD_IDAMEM_NULL);
107	return(IDAASCEND_MEM_NULL);
108	}
109	IDA_mem = (IDAMem)ida_mem;
110
111	if(lmem == NULL){
112	IDAProcessError(ida_mem, IDAASCEND_LMEM_NULL, "IDAASCEND", __FUNCTION__, MSGD_LMEM_NULL);
113	return(IDAASCEND_LMEM_NULL);
114	}
115	iamem = (IntegratorIdaAscendMem *)lmem;
116
117	jacfn = _jacfn;
118
119	return IDAASCEND_SUCCESS;
120	}
121
122	int IDAASCENDGetLastFlag(void ida_mem, int flag){
123	IDAMem IDA_mem;
124	IntegratorIdaAscendMem *iamem;
125
126	if(ida_mem == NULL){
127	IDAProcessError(NULL, IDAASCEND_MEM_NULL, "IDAASCEND", __FUNCTION__, MSGD_IDAMEM_NULL);
128	return(IDAASCEND_MEM_NULL);
129	}
130	IDA_mem = (IDAMem)ida_mem;
131
132	if(lmem == NULL){
133	IDAProcessError(ida_mem, IDAASCEND_LMEM_NULL, "IDAASCEND", __FUNCTION__, MSGD_LMEM_NULL);
134	return(IDAASCEND_LMEM_NULL);
135	}
136	iamem = (IntegratorIdaAscendMem *)lmem;
137
138	*flag = lastflag;
139
140	return IDAASCEND_SUCCESS;
141	}
142
143	char *IDAASCENDGetReturnFlagName(int flag){
144	char *name;
145
146	name = ASC_NEW_ARRAY(char,30);
147
148	#define FLAG(N) case N:sprintf(name,#N);break
149
150	switch(flag) {
151	FLAG(IDAASCEND_SUCCESS);
152	FLAG(IDAASCEND_MEM_NULL);
153	FLAG(IDAASCEND_LMEM_NULL);
154	FLAG(IDAASCEND_MEM_FAIL);
155	FLAG(IDAASCEND_JACFN_UNDEF);
156	FLAG(IDAASCEND_JACFN_UNRECVR);
157	FLAG(IDAASCEND_JACFN_RECVR);
158	default:
159	sprintf(name,"Unknown flag value '%d'",flag);
160	}
161
162	return name;
163	}
164
165
166	/*------------------------------------
167	Internal setup/evaluation routines required by IDA
168	*/
169
170	int integrator_ida_linit(IDAMem IDA_mem){
171	IntegratorIdaAscendMem *iamem;
172	iamem = (IntegratorIdaAscendMem *)lmem;
173
174	CONSOLE_DEBUG("Initialising IDA linear solver");
175	nje = 0;
176	nre = 0;
177	jacfn = NULL;
178
179	/* initialise anything else in the IntegratorIdaAscendMem that needs it */
180
181	return 0;
182	}
183
184	int integrator_ida_lsetup(IDAMem IDA_mem
185	, N_Vector yyp, N_Vector ypp, N_Vector rrp
186	, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3
187	){
188	int retval;
189
190	IntegratorIdaAscendMem *iamem;
191	iamem = (IntegratorIdaAscendMem *)lmem;
192
193	CONSOLE_DEBUG("Setting up IDA linear problem");
194
195	if(jacfn==NULL){
196	lastflag = IDAASCEND_JACFN_UNDEF;
197	return -1; /* unrecoverable */
198	}
199
200	/* Increment nje counter. */
201	nje++;
202
203	/* clear the jacobian matrix */
204	/* ... */
205
206	/* evaluate the jacobian */
207	retval = jacfn(neq, tn, yyp, ypp, rrp, cj, jacdata, JJ,
208	tmp1, tmp2, tmp3
209	);
210
211	if(retval < 0){
212	lastflag = IDAASCEND_JACFN_UNRECVR;
213	return -1;
214	}
215	if (retval > 0) {
216	lastflag = IDAASCEND_JACFN_RECVR;
217	return +1;
218	}
219
220	/* do block decomposition, LU factorisation or whatever, return success or fail flag */
221	/* ... */
222
223	ERROR_REPORTER_HERE(ASC_PROG_ERR,"Not implemented");
224	return(-1);
225	}
226
227	/**
228	This routines handles the linear solve operation for the IDAASCEND linear
229	solver. It interfaces to the appropriate mtx_matrix routines for this,
230	but also scales solution vector according to cjratio.
231
232	@return IDAASCEND_SUCESS on success
233	*/
234	int integrator_ida_lsolve(IDAMem IDA_mem
235	, N_Vector b, N_Vector weight
236	, N_Vector ycur, N_Vector ypcur, N_Vector rrcur
237	){
238	realtype *bd;
239
240	IntegratorIdaAscendMem *iamem;
241	iamem = (IntegratorIdaAscendMem *)lmem;
242
243	/* retrieve the data array for the RHS vector, 'b' */
244	bd = N_VGetArrayPointer(b);
245
246	/* call the necessary linsolqr routine here */
247	/* ... */
248	/* For IDADENSE, it's: DenseGETRS(JJ, pivots, bd); */
249
250	/* not sure what this is doing yet */
251	/* For IDADENSE, it's: Scale the correction to account for change in cj. */
252	if(cjratio != ONE){
253	N_VScale(TWO/(ONE + cjratio), b, b);
254	}
255
256	ERROR_REPORTER_HERE(ASC_PROG_ERR,"Not implemented");
257	return -1;
258	}
259
260	int integrator_ida_lfree(IDAMem IDA_mem){
261	CONSOLE_DEBUG("Freeing IDA linear solver data");
262
263	/* free jacobian mtx_matrix_t data ?? */
264
265	if(lmem!=NULL){
266	ASC_FREE(lmem);
267	lmem=NULL;
268	}
269	return 0;
270	}