1 /*
2 * -----------------------------------------------------------------
3 * Programmer(s): Allan Taylor, Alan Hindmarsh, Radu Serban, Carol Woodward,
4 * John Loffeld, and Aaron Collier @ LLNL
5 * -----------------------------------------------------------------
6 * SUNDIALS Copyright Start
7 * Copyright (c) 2002-2021, Lawrence Livermore National Security
8 * and Southern Methodist University.
9 * All rights reserved.
10 *
11 * See the top-level LICENSE and NOTICE files for details.
12 *
13 * SPDX-License-Identifier: BSD-3-Clause
14 * SUNDIALS Copyright End
15 * -----------------------------------------------------------------
16 * This is the implementation file for the main KINSol solver.
17 * It is independent of the KINSol linear solver in use.
18 * -----------------------------------------------------------------
19 *
20 * EXPORTED FUNCTIONS
21 * ------------------
22 * Creation and allocation functions
23 * KINCreate
24 * KINInit
25 * Main solver function
26 * KINSol
27 * Deallocation function
28 * KINFree
29 *
30 * PRIVATE FUNCTIONS
31 * -----------------
32 * KINCheckNvector
33 * Memory allocation/deallocation
34 * KINAllocVectors
35 * KINFreeVectors
36 * Initial setup
37 * KINSolInit
38 * Step functions
39 * KINLinSolDrv
40 * KINFullNewton
41 * KINLineSearch
42 * KINConstraint
43 * KINFP
44 * KINPicardAA
45 * Stopping tests
46 * KINStop
47 * KINForcingTerm
48 * Norm functions
49 * KINScFNorm
50 * KINScSNorm
51 * KINSOL Verbose output functions
52 * KINPrintInfo
53 * KINInfoHandler
54 * KINSOL Error Handling functions
55 * KINProcessError
56 * KINErrHandler
57 * -----------------------------------------------------------------
58 */
59
60 /*
61 * =================================================================
62 * IMPORTED HEADER FILES
63 * =================================================================
64 */
65
66 #include <stdio.h>
67 #include <stdlib.h>
68 #include <stdarg.h>
69 #include <string.h>
70
71 #include <math.h>
72
73 #include "kinsol_impl.h"
74 #include <sundials/sundials_math.h>
75
76 /*
77 * =================================================================
78 * KINSOL PRIVATE CONSTANTS
79 * =================================================================
80 */
81
82 #define HALF RCONST(0.5)
83 #define ZERO RCONST(0.0)
84 #define ONE RCONST(1.0)
85 #define ONEPT5 RCONST(1.5)
86 #define TWO RCONST(2.0)
87 #define THREE RCONST(3.0)
88 #define FIVE RCONST(5.0)
89 #define TWELVE RCONST(12.0)
90 #define POINT1 RCONST(0.1)
91 #define POINT01 RCONST(0.01)
92 #define POINT99 RCONST(0.99)
93 #define THOUSAND RCONST(1000.0)
94 #define ONETHIRD RCONST(0.3333333333333333)
95 #define TWOTHIRDS RCONST(0.6666666666666667)
96 #define POINT9 RCONST(0.9)
97 #define POINT0001 RCONST(0.0001)
98
99 /*
100 * =================================================================
101 * KINSOL ROUTINE-SPECIFIC CONSTANTS
102 * =================================================================
103 */
104
105 /*
106 * Control constants for lower-level functions used by KINSol
107 * ----------------------------------------------------------
108 *
109 * KINStop return value requesting more iterations
110 * RETRY_ITERATION
111 * CONTINUE_ITERATIONS
112 *
113 * KINFullNewton, KINLineSearch, KINFP, and KINPicardAA return values:
114 * KIN_SUCCESS
115 * KIN_SYSFUNC_FAIL
116 * STEP_TOO_SMALL
117 *
118 * KINConstraint return values:
119 * KIN_SUCCESS
120 * CONSTR_VIOLATED
121 */
122
123 #define RETRY_ITERATION -998
124 #define CONTINUE_ITERATIONS -999
125 #define STEP_TOO_SMALL -997
126 #define CONSTR_VIOLATED -996
127
128 /*
129 * Algorithmic constants
130 * ---------------------
131 *
132 * MAX_RECVR max. no. of attempts to correct a recoverable func error
133 */
134
135 #define MAX_RECVR 5
136
137 /*
138 * Keys for KINPrintInfo
139 * ---------------------
140 */
141
142 #define PRNT_RETVAL 1
143 #define PRNT_NNI 2
144 #define PRNT_TOL 3
145 #define PRNT_FMAX 4
146 #define PRNT_PNORM 5
147 #define PRNT_PNORM1 6
148 #define PRNT_FNORM 7
149 #define PRNT_LAM 8
150 #define PRNT_ALPHA 9
151 #define PRNT_BETA 10
152 #define PRNT_ALPHABETA 11
153 #define PRNT_ADJ 12
154
155 /*
156 * =================================================================
157 * PRIVATE FUNCTION PROTOTYPES
158 * =================================================================
159 */
160
161 static booleantype KINCheckNvector(N_Vector tmpl);
162 static booleantype KINAllocVectors(KINMem kin_mem, N_Vector tmpl);
163 static int KINSolInit(KINMem kin_mem);
164 static int KINConstraint(KINMem kin_mem );
165 static void KINForcingTerm(KINMem kin_mem, realtype fnormp);
166 static void KINFreeVectors(KINMem kin_mem);
167
168 static int KINFullNewton(KINMem kin_mem, realtype *fnormp,
169 realtype *f1normp, booleantype *maxStepTaken);
170 static int KINLineSearch(KINMem kin_mem, realtype *fnormp,
171 realtype *f1normp, booleantype *maxStepTaken);
172 static int KINPicardAA(KINMem kin_mem, long int *iter, realtype *R,
173 realtype *gamma, realtype *fmax);
174 static int KINFP(KINMem kin_mem);
175
176 static int KINLinSolDrv(KINMem kinmem);
177 static int KINPicardFcnEval(KINMem kin_mem, N_Vector gval, N_Vector uval,
178 N_Vector fval1);
179 static realtype KINScFNorm(KINMem kin_mem, N_Vector v, N_Vector scale);
180 static realtype KINScSNorm(KINMem kin_mem, N_Vector v, N_Vector u);
181 static int KINStop(KINMem kin_mem, booleantype maxStepTaken,
182 int sflag);
183 static int AndersonAcc(KINMem kin_mem, N_Vector gval, N_Vector fv, N_Vector x,
184 N_Vector x_old, long int iter, realtype *R, realtype *gamma);
185
186 /*
187 * =================================================================
188 * EXPORTED FUNCTIONS IMPLEMENTATION
189 * =================================================================
190 */
191
192 /*
193 * -----------------------------------------------------------------
194 * Creation and allocation functions
195 * -----------------------------------------------------------------
196 */
197
198 /*
199 * Function : KINCreate
200 *
201 * KINCreate creates an internal memory block for a problem to
202 * be solved by KINSOL. If successful, KINCreate returns a pointer
203 * to the problem memory. This pointer should be passed to
204 * KINInit. If an initialization error occurs, KINCreate prints
205 * an error message to standard error and returns NULL.
206 */
207
KINCreate(void)208 void *KINCreate(void)
209 {
210 KINMem kin_mem;
211 realtype uround;
212
213 kin_mem = NULL;
214 kin_mem = (KINMem) malloc(sizeof(struct KINMemRec));
215 if (kin_mem == NULL) {
216 KINProcessError(kin_mem, 0, "KINSOL", "KINCreate", MSG_MEM_FAIL);
217 return(NULL);
218 }
219
220 /* Zero out kin_mem */
221 memset(kin_mem, 0, sizeof(struct KINMemRec));
222
223 /* set uround (unit roundoff) */
224
225 kin_mem->kin_uround = uround = UNIT_ROUNDOFF;
226
227 /* set default values for solver optional inputs */
228
229 kin_mem->kin_func = NULL;
230 kin_mem->kin_user_data = NULL;
231 kin_mem->kin_uu = NULL;
232 kin_mem->kin_unew = NULL;
233 kin_mem->kin_fval = NULL;
234 kin_mem->kin_gval = NULL;
235 kin_mem->kin_uscale = NULL;
236 kin_mem->kin_fscale = NULL;
237 kin_mem->kin_pp = NULL;
238 kin_mem->kin_constraints = NULL;
239 kin_mem->kin_vtemp1 = NULL;
240 kin_mem->kin_vtemp2 = NULL;
241 kin_mem->kin_fold_aa = NULL;
242 kin_mem->kin_gold_aa = NULL;
243 kin_mem->kin_df_aa = NULL;
244 kin_mem->kin_dg_aa = NULL;
245 kin_mem->kin_q_aa = NULL;
246 kin_mem->kin_gamma_aa = NULL;
247 kin_mem->kin_R_aa = NULL;
248 kin_mem->kin_ipt_map = NULL;
249 kin_mem->kin_cv = NULL;
250 kin_mem->kin_Xv = NULL;
251 kin_mem->kin_lmem = NULL;
252 kin_mem->kin_m_aa = 0;
253 kin_mem->kin_aamem_aa = 0;
254 kin_mem->kin_setstop_aa = 0;
255 kin_mem->kin_beta_aa = ONE;
256 kin_mem->kin_damping_aa = SUNFALSE;
257 kin_mem->kin_constraintsSet = SUNFALSE;
258 kin_mem->kin_ehfun = KINErrHandler;
259 kin_mem->kin_eh_data = kin_mem;
260 kin_mem->kin_errfp = stderr;
261 kin_mem->kin_ihfun = KINInfoHandler;
262 kin_mem->kin_ih_data = kin_mem;
263 kin_mem->kin_infofp = stdout;
264 kin_mem->kin_printfl = PRINTFL_DEFAULT;
265 kin_mem->kin_mxiter = MXITER_DEFAULT;
266 kin_mem->kin_noInitSetup = SUNFALSE;
267 kin_mem->kin_msbset = MSBSET_DEFAULT;
268 kin_mem->kin_noResMon = SUNFALSE;
269 kin_mem->kin_msbset_sub = MSBSET_SUB_DEFAULT;
270 kin_mem->kin_update_fnorm_sub = SUNFALSE;
271 kin_mem->kin_mxnbcf = MXNBCF_DEFAULT;
272 kin_mem->kin_sthrsh = TWO;
273 kin_mem->kin_noMinEps = SUNFALSE;
274 kin_mem->kin_mxnstepin = ZERO;
275 kin_mem->kin_sqrt_relfunc = SUNRsqrt(uround);
276 kin_mem->kin_scsteptol = SUNRpowerR(uround,TWOTHIRDS);
277 kin_mem->kin_fnormtol = SUNRpowerR(uround,ONETHIRD);
278 kin_mem->kin_etaflag = KIN_ETACHOICE1;
279 kin_mem->kin_eta = POINT1; /* default for KIN_ETACONSTANT */
280 kin_mem->kin_eta_alpha = TWO; /* default for KIN_ETACHOICE2 */
281 kin_mem->kin_eta_gamma = POINT9; /* default for KIN_ETACHOICE2 */
282 kin_mem->kin_MallocDone = SUNFALSE;
283 kin_mem->kin_eval_omega = SUNTRUE;
284 kin_mem->kin_omega = ZERO; /* default to using min/max */
285 kin_mem->kin_omega_min = OMEGA_MIN;
286 kin_mem->kin_omega_max = OMEGA_MAX;
287
288 /* initialize lrw and liw */
289
290 kin_mem->kin_lrw = 17;
291 kin_mem->kin_liw = 22;
292
293 /* NOTE: needed since KINInit could be called after KINSetConstraints */
294
295 kin_mem->kin_lrw1 = 0;
296 kin_mem->kin_liw1 = 0;
297
298 return((void *) kin_mem);
299 }
300
301 /*
302 * Function : KINInit
303 *
304 * KINInit allocates memory for a problem or execution of KINSol.
305 * If memory is successfully allocated, KIN_SUCCESS is returned.
306 * Otherwise, an error message is printed and an error flag
307 * returned.
308 */
309
KINInit(void * kinmem,KINSysFn func,N_Vector tmpl)310 int KINInit(void *kinmem, KINSysFn func, N_Vector tmpl)
311 {
312 sunindextype liw1, lrw1;
313 KINMem kin_mem;
314 booleantype allocOK, nvectorOK;
315
316 /* check kinmem */
317
318 if (kinmem == NULL) {
319 KINProcessError(NULL, KIN_MEM_NULL, "KINSOL", "KINInit", MSG_NO_MEM);
320 return(KIN_MEM_NULL);
321 }
322 kin_mem = (KINMem) kinmem;
323
324 if (func == NULL) {
325 KINProcessError(kin_mem, KIN_ILL_INPUT, "KINSOL", "KINInit", MSG_FUNC_NULL);
326 return(KIN_ILL_INPUT);
327 }
328
329 /* check if all required vector operations are implemented */
330
331 nvectorOK = KINCheckNvector(tmpl);
332 if (!nvectorOK) {
333 KINProcessError(kin_mem, KIN_ILL_INPUT, "KINSOL", "KINInit", MSG_BAD_NVECTOR);
334 return(KIN_ILL_INPUT);
335 }
336
337 /* set space requirements for one N_Vector */
338
339 if (tmpl->ops->nvspace != NULL) {
340 N_VSpace(tmpl, &lrw1, &liw1);
341 kin_mem->kin_lrw1 = lrw1;
342 kin_mem->kin_liw1 = liw1;
343 }
344 else {
345 kin_mem->kin_lrw1 = 0;
346 kin_mem->kin_liw1 = 0;
347 }
348
349 /* allocate necessary vectors */
350
351 allocOK = KINAllocVectors(kin_mem, tmpl);
352 if (!allocOK) {
353 KINProcessError(kin_mem, KIN_MEM_FAIL, "KINSOL", "KINInit", MSG_MEM_FAIL);
354 free(kin_mem); kin_mem = NULL;
355 return(KIN_MEM_FAIL);
356 }
357
358 /* copy the input parameter into KINSol state */
359
360 kin_mem->kin_func = func;
361
362 /* set the linear solver addresses to NULL */
363
364 kin_mem->kin_linit = NULL;
365 kin_mem->kin_lsetup = NULL;
366 kin_mem->kin_lsolve = NULL;
367 kin_mem->kin_lfree = NULL;
368 kin_mem->kin_lmem = NULL;
369
370 /* problem memory has been successfully allocated */
371
372 kin_mem->kin_MallocDone = SUNTRUE;
373
374 return(KIN_SUCCESS);
375 }
376
377 /*
378 * -----------------------------------------------------------------
379 * Main solver function
380 * -----------------------------------------------------------------
381 */
382
383 /*
384 * Function : KINSol
385 *
386 * KINSol (main KINSOL driver routine) manages the computational
387 * process of computing an approximate solution of the nonlinear
388 * system F(uu) = 0. The KINSol routine calls the following
389 * subroutines:
390 *
391 * KINSolInit checks if initial guess satisfies user-supplied
392 * constraints and initializes linear solver
393 *
394 * KINLinSolDrv interfaces with linear solver to find a
395 * solution of the system J(uu)*x = b (calculate
396 * Newton step)
397 *
398 * KINFullNewton/KINLineSearch implement the global strategy
399 *
400 * KINForcingTerm computes the forcing term (eta)
401 *
402 * KINStop determines if an approximate solution has been found
403 */
404
KINSol(void * kinmem,N_Vector u,int strategy_in,N_Vector u_scale,N_Vector f_scale)405 int KINSol(void *kinmem, N_Vector u, int strategy_in,
406 N_Vector u_scale, N_Vector f_scale)
407 {
408 realtype fnormp, f1normp, epsmin, fmax=ZERO;
409 KINMem kin_mem;
410 int ret, sflag;
411 booleantype maxStepTaken;
412
413 /* intialize to avoid compiler warning messages */
414
415 maxStepTaken = SUNFALSE;
416 f1normp = fnormp = -ONE;
417
418 /* initialize epsmin to avoid compiler warning message */
419
420 epsmin = ZERO;
421
422 /* check for kinmem non-NULL */
423
424 if (kinmem == NULL) {
425 KINProcessError(NULL, KIN_MEM_NULL, "KINSOL", "KINSol", MSG_NO_MEM);
426 return(KIN_MEM_NULL);
427 }
428 kin_mem = (KINMem) kinmem;
429
430 if(kin_mem->kin_MallocDone == SUNFALSE) {
431 KINProcessError(NULL, KIN_NO_MALLOC, "KINSOL", "KINSol", MSG_NO_MALLOC);
432 return(KIN_NO_MALLOC);
433 }
434
435 /* load input arguments */
436
437 kin_mem->kin_uu = u;
438 kin_mem->kin_uscale = u_scale;
439 kin_mem->kin_fscale = f_scale;
440 kin_mem->kin_globalstrategy = strategy_in;
441
442 /* CSW:
443 Call fixed point solver if requested. Note that this should probably
444 be forked off to a FPSOL solver instead of kinsol in the future. */
445 if ( kin_mem->kin_globalstrategy == KIN_FP ) {
446 if (kin_mem->kin_uu == NULL) {
447 KINProcessError(kin_mem, KIN_ILL_INPUT, "KINSOL", "KINSol", MSG_UU_NULL);
448 return(KIN_ILL_INPUT);
449 }
450
451 if (kin_mem->kin_constraintsSet != SUNFALSE) {
452 KINProcessError(kin_mem, KIN_ILL_INPUT, "KINSOL", "KINSol", MSG_CONSTRAINTS_NOTOK);
453 return(KIN_ILL_INPUT);
454 }
455
456 if (kin_mem->kin_printfl > 0)
457 KINPrintInfo(kin_mem, PRNT_TOL, "KINSOL", "KINSol", INFO_TOL, kin_mem->kin_scsteptol, kin_mem->kin_fnormtol);
458
459 kin_mem->kin_nfe = kin_mem->kin_nnilset = kin_mem->kin_nnilset_sub = kin_mem->kin_nni = kin_mem->kin_nbcf = kin_mem->kin_nbktrk = 0;
460 ret = KINFP(kin_mem);
461
462 switch(ret) {
463 case KIN_SYSFUNC_FAIL:
464 KINProcessError(kin_mem, KIN_SYSFUNC_FAIL, "KINSOL", "KINSol", MSG_SYSFUNC_FAILED);
465 break;
466 case KIN_MAXITER_REACHED:
467 KINProcessError(kin_mem, KIN_MAXITER_REACHED, "KINSOL", "KINSol", MSG_MAXITER_REACHED);
468 break;
469 }
470
471 return(ret);
472 }
473
474 /* initialize solver */
475 ret = KINSolInit(kin_mem);
476 if (ret != KIN_SUCCESS) return(ret);
477
478 kin_mem->kin_ncscmx = 0;
479
480 /* Note: The following logic allows the choice of whether or not
481 to force a call to the linear solver setup upon a given call to
482 KINSol */
483
484 if (kin_mem->kin_noInitSetup) kin_mem->kin_sthrsh = ONE;
485 else kin_mem->kin_sthrsh = TWO;
486
487 /* if eps is to be bounded from below, set the bound */
488
489 if (kin_mem->kin_inexact_ls && !(kin_mem->kin_noMinEps))
490 epsmin = POINT01 * kin_mem->kin_fnormtol;
491
492
493 /* if omega is zero at this point, make sure it will be evaluated
494 at each iteration based on the provided min/max bounds and the
495 current function norm. */
496 if (kin_mem->kin_omega == ZERO) kin_mem->kin_eval_omega = SUNTRUE;
497 else kin_mem->kin_eval_omega = SUNFALSE;
498
499
500 /* CSW:
501 Call fixed point solver for Picard method if requested.
502 Note that this should probably be forked off to a part of an
503 FPSOL solver instead of kinsol in the future. */
504 if ( kin_mem->kin_globalstrategy == KIN_PICARD ) {
505
506 if (kin_mem->kin_gval == NULL) {
507 kin_mem->kin_gval = N_VClone(kin_mem->kin_unew);
508 if (kin_mem->kin_gval == NULL) {
509 KINProcessError(kin_mem, KIN_MEM_FAIL, "KINSOL", "KINSol", MSG_MEM_FAIL);
510 return(KIN_MEM_FAIL);
511 }
512 kin_mem->kin_liw += kin_mem->kin_liw1;
513 kin_mem->kin_lrw += kin_mem->kin_lrw1;
514 }
515 ret = KINPicardAA(kin_mem, &(kin_mem->kin_nni), kin_mem->kin_R_aa, kin_mem->kin_gamma_aa, &fmax);
516
517 return(ret);
518 }
519
520
521 for(;;){
522
523 kin_mem->kin_retry_nni = SUNFALSE;
524
525 kin_mem->kin_nni++;
526
527 /* calculate the epsilon (stopping criteria for iterative linear solver)
528 for this iteration based on eta from the routine KINForcingTerm */
529
530 if (kin_mem->kin_inexact_ls) {
531 kin_mem->kin_eps = (kin_mem->kin_eta + kin_mem->kin_uround) * kin_mem->kin_fnorm;
532 if(!(kin_mem->kin_noMinEps)) kin_mem->kin_eps = SUNMAX(epsmin, kin_mem->kin_eps);
533 }
534
535 repeat_nni:
536
537 /* call the appropriate routine to calculate an acceptable step pp */
538
539 sflag = 0;
540
541 if (kin_mem->kin_globalstrategy == KIN_NONE) {
542
543 /* Full Newton Step*/
544
545 /* call KINLinSolDrv to calculate the (approximate) Newton step, pp */
546 ret = KINLinSolDrv(kin_mem);
547 if (ret != KIN_SUCCESS) break;
548
549 sflag = KINFullNewton(kin_mem, &fnormp, &f1normp, &maxStepTaken);
550
551 /* if sysfunc failed unrecoverably, stop */
552 if ((sflag == KIN_SYSFUNC_FAIL) || (sflag == KIN_REPTD_SYSFUNC_ERR)) {
553 ret = sflag;
554 break;
555 }
556
557 } else if (kin_mem->kin_globalstrategy == KIN_LINESEARCH) {
558
559 /* Line Search */
560
561 /* call KINLinSolDrv to calculate the (approximate) Newton step, pp */
562 ret = KINLinSolDrv(kin_mem);
563 if (ret != KIN_SUCCESS) break;
564
565 sflag = KINLineSearch(kin_mem, &fnormp, &f1normp, &maxStepTaken);
566
567 /* if sysfunc failed unrecoverably, stop */
568 if ((sflag == KIN_SYSFUNC_FAIL) || (sflag == KIN_REPTD_SYSFUNC_ERR)) {
569 ret = sflag;
570 break;
571 }
572
573 /* if too many beta condition failures, then stop iteration */
574 if (kin_mem->kin_nbcf > kin_mem->kin_mxnbcf) {
575 ret = KIN_LINESEARCH_BCFAIL;
576 break;
577 }
578
579 }
580
581 if ( (kin_mem->kin_globalstrategy != KIN_PICARD) &&
582 (kin_mem->kin_globalstrategy != KIN_FP) ) {
583
584 /* evaluate eta by calling the forcing term routine */
585 if (kin_mem->kin_callForcingTerm) KINForcingTerm(kin_mem, fnormp);
586
587 kin_mem->kin_fnorm = fnormp;
588
589 /* call KINStop to check if tolerances where met by this iteration */
590 ret = KINStop(kin_mem, maxStepTaken, sflag);
591
592 if (ret == RETRY_ITERATION) {
593 kin_mem->kin_retry_nni = SUNTRUE;
594 goto repeat_nni;
595 }
596 }
597
598 /* update uu after the iteration */
599 N_VScale(ONE, kin_mem->kin_unew, kin_mem->kin_uu);
600
601 kin_mem->kin_f1norm = f1normp;
602
603 /* print the current nni, fnorm, and nfe values if printfl > 0 */
604
605 if (kin_mem->kin_printfl > 0)
606 KINPrintInfo(kin_mem, PRNT_NNI, "KINSOL", "KINSol", INFO_NNI, kin_mem->kin_nni, kin_mem->kin_nfe, kin_mem->kin_fnorm);
607
608 if (ret != CONTINUE_ITERATIONS) break;
609
610 fflush(kin_mem->kin_errfp);
611
612 } /* end of loop; return */
613
614
615
616 if (kin_mem->kin_printfl > 0)
617 KINPrintInfo(kin_mem, PRNT_RETVAL, "KINSOL", "KINSol", INFO_RETVAL, ret);
618
619 switch(ret) {
620 case KIN_SYSFUNC_FAIL:
621 KINProcessError(kin_mem, KIN_SYSFUNC_FAIL, "KINSOL", "KINSol", MSG_SYSFUNC_FAILED);
622 break;
623 case KIN_REPTD_SYSFUNC_ERR:
624 KINProcessError(kin_mem, KIN_REPTD_SYSFUNC_ERR, "KINSOL", "KINSol", MSG_SYSFUNC_REPTD);
625 break;
626 case KIN_LSETUP_FAIL:
627 KINProcessError(kin_mem, KIN_LSETUP_FAIL, "KINSOL", "KINSol", MSG_LSETUP_FAILED);
628 break;
629 case KIN_LSOLVE_FAIL:
630 KINProcessError(kin_mem, KIN_LSOLVE_FAIL, "KINSOL", "KINSol", MSG_LSOLVE_FAILED);
631 break;
632 case KIN_LINSOLV_NO_RECOVERY:
633 KINProcessError(kin_mem, KIN_LINSOLV_NO_RECOVERY, "KINSOL", "KINSol", MSG_LINSOLV_NO_RECOVERY);
634 break;
635 case KIN_LINESEARCH_NONCONV:
636 KINProcessError(kin_mem, KIN_LINESEARCH_NONCONV, "KINSOL", "KINSol", MSG_LINESEARCH_NONCONV);
637 break;
638 case KIN_LINESEARCH_BCFAIL:
639 KINProcessError(kin_mem, KIN_LINESEARCH_BCFAIL, "KINSOL", "KINSol", MSG_LINESEARCH_BCFAIL);
640 break;
641 case KIN_MAXITER_REACHED:
642 KINProcessError(kin_mem, KIN_MAXITER_REACHED, "KINSOL", "KINSol", MSG_MAXITER_REACHED);
643 break;
644 case KIN_MXNEWT_5X_EXCEEDED:
645 KINProcessError(kin_mem, KIN_MXNEWT_5X_EXCEEDED, "KINSOL", "KINSol", MSG_MXNEWT_5X_EXCEEDED);
646 break;
647 }
648
649 return(ret);
650 }
651
652 /*
653 * -----------------------------------------------------------------
654 * Deallocation function
655 * -----------------------------------------------------------------
656 */
657
658 /*
659 * Function : KINFree
660 *
661 * This routine frees the problem memory allocated by KINInit.
662 * Such memory includes all the vectors allocated by
663 * KINAllocVectors, and the memory lmem for the linear solver
664 * (deallocated by a call to lfree).
665 */
666
KINFree(void ** kinmem)667 void KINFree(void **kinmem)
668 {
669 KINMem kin_mem;
670
671 if (*kinmem == NULL) return;
672
673 kin_mem = (KINMem) (*kinmem);
674 KINFreeVectors(kin_mem);
675
676 /* call lfree if non-NULL */
677
678 if (kin_mem->kin_lfree != NULL) kin_mem->kin_lfree(kin_mem);
679
680 free(*kinmem);
681 *kinmem = NULL;
682 }
683
684 /*
685 * =================================================================
686 * PRIVATE FUNCTIONS
687 * =================================================================
688 */
689
690 /*
691 * Function : KINCheckNvector
692 *
693 * This routine checks if all required vector operations are
694 * implemented (excluding those required by KINConstraint). If all
695 * necessary operations are present, then KINCheckNvector returns
696 * SUNTRUE. Otherwise, SUNFALSE is returned.
697 */
698
KINCheckNvector(N_Vector tmpl)699 static booleantype KINCheckNvector(N_Vector tmpl)
700 {
701 if ((tmpl->ops->nvclone == NULL) ||
702 (tmpl->ops->nvdestroy == NULL) ||
703 (tmpl->ops->nvlinearsum == NULL) ||
704 (tmpl->ops->nvprod == NULL) ||
705 (tmpl->ops->nvdiv == NULL) ||
706 (tmpl->ops->nvscale == NULL) ||
707 (tmpl->ops->nvabs == NULL) ||
708 (tmpl->ops->nvinv == NULL) ||
709 (tmpl->ops->nvmaxnorm == NULL) ||
710 (tmpl->ops->nvmin == NULL) ||
711 (tmpl->ops->nvwl2norm == NULL)) return(SUNFALSE);
712 else return(SUNTRUE);
713 }
714
715 /*
716 * -----------------------------------------------------------------
717 * Memory allocation/deallocation
718 * -----------------------------------------------------------------
719 */
720
721 /*
722 * Function : KINAllocVectors
723 *
724 * This routine allocates the KINSol vectors. If all memory
725 * allocations are successful, KINAllocVectors returns SUNTRUE.
726 * Otherwise all allocated memory is freed and KINAllocVectors
727 * returns SUNFALSE.
728 */
729
KINAllocVectors(KINMem kin_mem,N_Vector tmpl)730 static booleantype KINAllocVectors(KINMem kin_mem, N_Vector tmpl)
731 {
732 /* allocate unew, fval, pp, vtemp1 and vtemp2. */
733 /* allocate df, dg, q, for Anderson Acceleration, Broyden and EN */
734
735 if (kin_mem->kin_unew == NULL) {
736 kin_mem->kin_unew = N_VClone(tmpl);
737 if (kin_mem->kin_unew == NULL) return(SUNFALSE);
738 kin_mem->kin_liw += kin_mem->kin_liw1;
739 kin_mem->kin_lrw += kin_mem->kin_lrw1;
740 }
741
742 if (kin_mem->kin_fval == NULL) {
743 kin_mem->kin_fval = N_VClone(tmpl);
744 if (kin_mem->kin_fval == NULL) {
745 N_VDestroy(kin_mem->kin_unew);
746 kin_mem->kin_liw -= kin_mem->kin_liw1;
747 kin_mem->kin_lrw -= kin_mem->kin_lrw1;
748 return(SUNFALSE);
749 }
750 kin_mem->kin_liw += kin_mem->kin_liw1;
751 kin_mem->kin_lrw += kin_mem->kin_lrw1;
752 }
753
754 if (kin_mem->kin_pp == NULL) {
755 kin_mem->kin_pp = N_VClone(tmpl);
756 if (kin_mem->kin_pp == NULL) {
757 N_VDestroy(kin_mem->kin_unew);
758 N_VDestroy(kin_mem->kin_fval);
759 kin_mem->kin_liw -= 2*kin_mem->kin_liw1;
760 kin_mem->kin_lrw -= 2*kin_mem->kin_lrw1;
761 return(SUNFALSE);
762 }
763 kin_mem->kin_liw += kin_mem->kin_liw1;
764 kin_mem->kin_lrw += kin_mem->kin_lrw1;
765 }
766
767 if (kin_mem->kin_vtemp1 == NULL) {
768 kin_mem->kin_vtemp1 = N_VClone(tmpl);
769 if (kin_mem->kin_vtemp1 == NULL) {
770 N_VDestroy(kin_mem->kin_unew);
771 N_VDestroy(kin_mem->kin_fval);
772 N_VDestroy(kin_mem->kin_pp);
773 kin_mem->kin_liw -= 3*kin_mem->kin_liw1;
774 kin_mem->kin_lrw -= 3*kin_mem->kin_lrw1;
775 return(SUNFALSE);
776 }
777 kin_mem->kin_liw += kin_mem->kin_liw1;
778 kin_mem->kin_lrw += kin_mem->kin_lrw1;
779 }
780
781 if (kin_mem->kin_vtemp2 == NULL) {
782 kin_mem->kin_vtemp2 = N_VClone(tmpl);
783 if (kin_mem->kin_vtemp2 == NULL) {
784 N_VDestroy(kin_mem->kin_unew);
785 N_VDestroy(kin_mem->kin_fval);
786 N_VDestroy(kin_mem->kin_pp);
787 N_VDestroy(kin_mem->kin_vtemp1);
788 kin_mem->kin_liw -= 4*kin_mem->kin_liw1;
789 kin_mem->kin_lrw -= 4*kin_mem->kin_lrw1;
790 return(SUNFALSE);
791 }
792 kin_mem->kin_liw += kin_mem->kin_liw1;
793 kin_mem->kin_lrw += kin_mem->kin_lrw1;
794 }
795
796 /* Vectors for Anderson acceleration */
797
798 if (kin_mem->kin_m_aa) {
799
800 if (kin_mem->kin_R_aa == NULL) {
801 kin_mem->kin_R_aa = (realtype *) malloc((kin_mem->kin_m_aa*kin_mem->kin_m_aa) * sizeof(realtype));
802 if (kin_mem->kin_R_aa == NULL) {
803 KINProcessError(kin_mem, 0, "KINSOL", "KINAllocVectors", MSG_MEM_FAIL);
804 N_VDestroy(kin_mem->kin_unew);
805 N_VDestroy(kin_mem->kin_fval);
806 N_VDestroy(kin_mem->kin_pp);
807 N_VDestroy(kin_mem->kin_vtemp1);
808 N_VDestroy(kin_mem->kin_vtemp2);
809 kin_mem->kin_liw -= 5*kin_mem->kin_liw1;
810 kin_mem->kin_lrw -= 5*kin_mem->kin_lrw1;
811 return(KIN_MEM_FAIL);
812 }
813 }
814
815 if (kin_mem->kin_gamma_aa == NULL) {
816 kin_mem->kin_gamma_aa = (realtype *) malloc(kin_mem->kin_m_aa * sizeof(realtype));
817 if (kin_mem->kin_gamma_aa == NULL) {
818 KINProcessError(kin_mem, 0, "KINSOL", "KINAllocVectors", MSG_MEM_FAIL);
819 N_VDestroy(kin_mem->kin_unew);
820 N_VDestroy(kin_mem->kin_fval);
821 N_VDestroy(kin_mem->kin_pp);
822 N_VDestroy(kin_mem->kin_vtemp1);
823 N_VDestroy(kin_mem->kin_vtemp2);
824 free(kin_mem->kin_R_aa);
825 kin_mem->kin_liw -= 5*kin_mem->kin_liw1;
826 kin_mem->kin_lrw -= 5*kin_mem->kin_lrw1;
827 return(KIN_MEM_FAIL);
828 }
829 }
830
831 if (kin_mem->kin_ipt_map == NULL) {
832 kin_mem->kin_ipt_map = (long int *) malloc(kin_mem->kin_m_aa * sizeof(long int));
833 if (kin_mem->kin_ipt_map == NULL) {
834 KINProcessError(kin_mem, 0, "KINSOL", "KINAllocVectors", MSG_MEM_FAIL);
835 N_VDestroy(kin_mem->kin_unew);
836 N_VDestroy(kin_mem->kin_fval);
837 N_VDestroy(kin_mem->kin_pp);
838 N_VDestroy(kin_mem->kin_vtemp1);
839 N_VDestroy(kin_mem->kin_vtemp2);
840 free(kin_mem->kin_R_aa);
841 free(kin_mem->kin_gamma_aa);
842 kin_mem->kin_liw -= 5*kin_mem->kin_liw1;
843 kin_mem->kin_lrw -= 5*kin_mem->kin_lrw1;
844 return(KIN_MEM_FAIL);
845 }
846 }
847
848 if (kin_mem->kin_cv == NULL) {
849 kin_mem->kin_cv = (realtype *) malloc(2 * (kin_mem->kin_m_aa+1) * sizeof(realtype));
850 if (kin_mem->kin_cv == NULL) {
851 KINProcessError(kin_mem, 0, "KINSOL", "KINAllocVectors", MSG_MEM_FAIL);
852 N_VDestroy(kin_mem->kin_unew);
853 N_VDestroy(kin_mem->kin_fval);
854 N_VDestroy(kin_mem->kin_pp);
855 N_VDestroy(kin_mem->kin_vtemp1);
856 N_VDestroy(kin_mem->kin_vtemp2);
857 free(kin_mem->kin_R_aa);
858 free(kin_mem->kin_gamma_aa);
859 free(kin_mem->kin_ipt_map);
860 kin_mem->kin_liw -= 5*kin_mem->kin_liw1;
861 kin_mem->kin_lrw -= 5*kin_mem->kin_lrw1;
862 return(KIN_MEM_FAIL);
863 }
864 }
865
866 if (kin_mem->kin_Xv == NULL) {
867 kin_mem->kin_Xv = (N_Vector *) malloc(2 * (kin_mem->kin_m_aa+1) * sizeof(N_Vector));
868 if (kin_mem->kin_Xv == NULL) {
869 KINProcessError(kin_mem, 0, "KINSOL", "KINAllocVectors", MSG_MEM_FAIL);
870 N_VDestroy(kin_mem->kin_unew);
871 N_VDestroy(kin_mem->kin_fval);
872 N_VDestroy(kin_mem->kin_pp);
873 N_VDestroy(kin_mem->kin_vtemp1);
874 N_VDestroy(kin_mem->kin_vtemp2);
875 free(kin_mem->kin_R_aa);
876 free(kin_mem->kin_gamma_aa);
877 free(kin_mem->kin_ipt_map);
878 free(kin_mem->kin_cv);
879 kin_mem->kin_liw -= 5*kin_mem->kin_liw1;
880 kin_mem->kin_lrw -= 5*kin_mem->kin_lrw1;
881 return(KIN_MEM_FAIL);
882 }
883 }
884
885 if (kin_mem->kin_fold_aa == NULL) {
886 kin_mem->kin_fold_aa = N_VClone(tmpl);
887 if (kin_mem->kin_fold_aa == NULL) {
888 N_VDestroy(kin_mem->kin_unew);
889 N_VDestroy(kin_mem->kin_fval);
890 N_VDestroy(kin_mem->kin_pp);
891 N_VDestroy(kin_mem->kin_vtemp1);
892 N_VDestroy(kin_mem->kin_vtemp2);
893 free(kin_mem->kin_R_aa);
894 free(kin_mem->kin_gamma_aa);
895 free(kin_mem->kin_ipt_map);
896 free(kin_mem->kin_cv);
897 free(kin_mem->kin_Xv);
898 kin_mem->kin_liw -= 5*kin_mem->kin_liw1;
899 kin_mem->kin_lrw -= 5*kin_mem->kin_lrw1;
900 return(SUNFALSE);
901 }
902 kin_mem->kin_liw += kin_mem->kin_liw1;
903 kin_mem->kin_lrw += kin_mem->kin_lrw1;
904 }
905
906 if (kin_mem->kin_gold_aa == NULL) {
907 kin_mem->kin_gold_aa = N_VClone(tmpl);
908 if (kin_mem->kin_gold_aa == NULL) {
909 N_VDestroy(kin_mem->kin_unew);
910 N_VDestroy(kin_mem->kin_fval);
911 N_VDestroy(kin_mem->kin_pp);
912 N_VDestroy(kin_mem->kin_vtemp1);
913 N_VDestroy(kin_mem->kin_vtemp2);
914 free(kin_mem->kin_R_aa);
915 free(kin_mem->kin_gamma_aa);
916 free(kin_mem->kin_ipt_map);
917 free(kin_mem->kin_cv);
918 free(kin_mem->kin_Xv);
919 N_VDestroy(kin_mem->kin_fold_aa);
920 kin_mem->kin_liw -= 6*kin_mem->kin_liw1;
921 kin_mem->kin_lrw -= 6*kin_mem->kin_lrw1;
922 return(SUNFALSE);
923 }
924 kin_mem->kin_liw += kin_mem->kin_liw1;
925 kin_mem->kin_lrw += kin_mem->kin_lrw1;
926 }
927
928 if (kin_mem->kin_df_aa == NULL) {
929 kin_mem->kin_df_aa = N_VCloneVectorArray((int) kin_mem->kin_m_aa,tmpl);
930 if (kin_mem->kin_df_aa == NULL) {
931 N_VDestroy(kin_mem->kin_unew);
932 N_VDestroy(kin_mem->kin_fval);
933 N_VDestroy(kin_mem->kin_pp);
934 N_VDestroy(kin_mem->kin_vtemp1);
935 N_VDestroy(kin_mem->kin_vtemp2);
936 free(kin_mem->kin_R_aa);
937 free(kin_mem->kin_gamma_aa);
938 free(kin_mem->kin_ipt_map);
939 free(kin_mem->kin_cv);
940 free(kin_mem->kin_Xv);
941 N_VDestroy(kin_mem->kin_fold_aa);
942 N_VDestroy(kin_mem->kin_gold_aa);
943 kin_mem->kin_liw -= 7*kin_mem->kin_liw1;
944 kin_mem->kin_lrw -= 7*kin_mem->kin_lrw1;
945 return(SUNFALSE);
946 }
947 kin_mem->kin_liw += kin_mem->kin_m_aa * kin_mem->kin_liw1;
948 kin_mem->kin_lrw += kin_mem->kin_m_aa * kin_mem->kin_lrw1;
949 }
950
951 if (kin_mem->kin_dg_aa == NULL) {
952 kin_mem->kin_dg_aa = N_VCloneVectorArray((int) kin_mem->kin_m_aa,tmpl);
953 if (kin_mem->kin_dg_aa == NULL) {
954 N_VDestroy(kin_mem->kin_unew);
955 N_VDestroy(kin_mem->kin_fval);
956 N_VDestroy(kin_mem->kin_pp);
957 N_VDestroy(kin_mem->kin_vtemp1);
958 N_VDestroy(kin_mem->kin_vtemp2);
959 free(kin_mem->kin_R_aa);
960 free(kin_mem->kin_gamma_aa);
961 free(kin_mem->kin_ipt_map);
962 free(kin_mem->kin_cv);
963 free(kin_mem->kin_Xv);
964 N_VDestroy(kin_mem->kin_fold_aa);
965 N_VDestroy(kin_mem->kin_gold_aa);
966 N_VDestroyVectorArray(kin_mem->kin_df_aa, (int) kin_mem->kin_m_aa);
967 kin_mem->kin_liw -= (7 + kin_mem->kin_m_aa) * kin_mem->kin_liw1;
968 kin_mem->kin_lrw -= (7 + kin_mem->kin_m_aa) * kin_mem->kin_lrw1;
969 return(SUNFALSE);
970 }
971 kin_mem->kin_liw += kin_mem->kin_m_aa * kin_mem->kin_liw1;
972 kin_mem->kin_lrw += kin_mem->kin_m_aa * kin_mem->kin_lrw1;
973 }
974
975 if (kin_mem->kin_q_aa == NULL) {
976 kin_mem->kin_q_aa = N_VCloneVectorArray((int) kin_mem->kin_m_aa,tmpl);
977 if (kin_mem->kin_q_aa == NULL) {
978 N_VDestroy(kin_mem->kin_unew);
979 N_VDestroy(kin_mem->kin_fval);
980 N_VDestroy(kin_mem->kin_pp);
981 N_VDestroy(kin_mem->kin_vtemp1);
982 N_VDestroy(kin_mem->kin_vtemp2);
983 free(kin_mem->kin_R_aa);
984 free(kin_mem->kin_gamma_aa);
985 free(kin_mem->kin_ipt_map);
986 free(kin_mem->kin_cv);
987 free(kin_mem->kin_Xv);
988 N_VDestroy(kin_mem->kin_fold_aa);
989 N_VDestroy(kin_mem->kin_gold_aa);
990 N_VDestroyVectorArray(kin_mem->kin_df_aa, (int) kin_mem->kin_m_aa);
991 N_VDestroyVectorArray(kin_mem->kin_dg_aa, (int) kin_mem->kin_m_aa);
992 kin_mem->kin_liw -= (7 + 2 * kin_mem->kin_m_aa) * kin_mem->kin_liw1;
993 kin_mem->kin_lrw -= (7 + 2 * kin_mem->kin_m_aa) * kin_mem->kin_lrw1;
994 return(SUNFALSE);
995 }
996 kin_mem->kin_liw += kin_mem->kin_m_aa * kin_mem->kin_liw1;
997 kin_mem->kin_lrw += kin_mem->kin_m_aa * kin_mem->kin_lrw1;
998 }
999 }
1000
1001 return(SUNTRUE);
1002 }
1003
1004 /*
1005 * KINFreeVectors
1006 *
1007 * This routine frees the KINSol vectors allocated by
1008 * KINAllocVectors.
1009 */
1010
KINFreeVectors(KINMem kin_mem)1011 static void KINFreeVectors(KINMem kin_mem)
1012 {
1013 if (kin_mem->kin_unew != NULL) {
1014 N_VDestroy(kin_mem->kin_unew);
1015 kin_mem->kin_unew = NULL;
1016 kin_mem->kin_lrw -= kin_mem->kin_lrw1;
1017 kin_mem->kin_liw -= kin_mem->kin_liw1;
1018 }
1019
1020 if (kin_mem->kin_fval != NULL) {
1021 N_VDestroy(kin_mem->kin_fval);
1022 kin_mem->kin_fval = NULL;
1023 kin_mem->kin_lrw -= kin_mem->kin_lrw1;
1024 kin_mem->kin_liw -= kin_mem->kin_liw1;
1025 }
1026
1027 if (kin_mem->kin_pp != NULL) {
1028 N_VDestroy(kin_mem->kin_pp);
1029 kin_mem->kin_pp = NULL;
1030 kin_mem->kin_lrw -= kin_mem->kin_lrw1;
1031 kin_mem->kin_liw -= kin_mem->kin_liw1;
1032 }
1033
1034 if (kin_mem->kin_vtemp1 != NULL) {
1035 N_VDestroy(kin_mem->kin_vtemp1);
1036 kin_mem->kin_vtemp1 = NULL;
1037 kin_mem->kin_lrw -= kin_mem->kin_lrw1;
1038 kin_mem->kin_liw -= kin_mem->kin_liw1;
1039 }
1040
1041 if (kin_mem->kin_vtemp2 != NULL) {
1042 N_VDestroy(kin_mem->kin_vtemp2);
1043 kin_mem->kin_vtemp2 = NULL;
1044 kin_mem->kin_lrw -= kin_mem->kin_lrw1;
1045 kin_mem->kin_liw -= kin_mem->kin_liw1;
1046 }
1047
1048 if (kin_mem->kin_gval != NULL) {
1049 N_VDestroy(kin_mem->kin_gval);
1050 kin_mem->kin_gval = NULL;
1051 kin_mem->kin_lrw -= kin_mem->kin_lrw1;
1052 kin_mem->kin_liw -= kin_mem->kin_liw1;
1053 }
1054
1055 if (kin_mem->kin_R_aa != NULL) {
1056 free(kin_mem->kin_R_aa);
1057 kin_mem->kin_R_aa = NULL;
1058 }
1059
1060 if (kin_mem->kin_gamma_aa != NULL) {
1061 free(kin_mem->kin_gamma_aa);
1062 kin_mem->kin_gamma_aa = NULL;
1063 }
1064
1065 if (kin_mem->kin_ipt_map != NULL) {
1066 free(kin_mem->kin_ipt_map);
1067 kin_mem->kin_ipt_map = NULL;
1068 }
1069
1070 if (kin_mem->kin_cv != NULL) {
1071 free(kin_mem->kin_cv);
1072 kin_mem->kin_cv = NULL;
1073 }
1074
1075 if (kin_mem->kin_Xv != NULL) {
1076 free(kin_mem->kin_Xv);
1077 kin_mem->kin_Xv = NULL;
1078 }
1079
1080 if (kin_mem->kin_fold_aa != NULL) {
1081 N_VDestroy(kin_mem->kin_fold_aa);
1082 kin_mem->kin_fold_aa = NULL;
1083 kin_mem->kin_lrw -= kin_mem->kin_lrw1;
1084 kin_mem->kin_liw -= kin_mem->kin_liw1;
1085 }
1086
1087 if (kin_mem->kin_gold_aa != NULL) {
1088 N_VDestroy(kin_mem->kin_gold_aa);
1089 kin_mem->kin_gold_aa = NULL;
1090 kin_mem->kin_lrw -= kin_mem->kin_lrw1;
1091 kin_mem->kin_liw -= kin_mem->kin_liw1;
1092 }
1093
1094 if (kin_mem->kin_df_aa != NULL) {
1095 N_VDestroyVectorArray(kin_mem->kin_df_aa, (int) kin_mem->kin_m_aa);
1096 kin_mem->kin_df_aa = NULL;
1097 kin_mem->kin_lrw -= kin_mem->kin_m_aa * kin_mem->kin_lrw1;
1098 kin_mem->kin_liw -= kin_mem->kin_m_aa * kin_mem->kin_liw1;
1099 }
1100
1101 if (kin_mem->kin_dg_aa != NULL) {
1102 N_VDestroyVectorArray(kin_mem->kin_dg_aa, (int) kin_mem->kin_m_aa);
1103 kin_mem->kin_dg_aa = NULL;
1104 kin_mem->kin_lrw -= kin_mem->kin_m_aa * kin_mem->kin_lrw1;
1105 kin_mem->kin_liw -= kin_mem->kin_m_aa * kin_mem->kin_liw1;
1106 }
1107
1108 if (kin_mem->kin_q_aa != NULL) {
1109 N_VDestroyVectorArray(kin_mem->kin_q_aa, (int) kin_mem->kin_m_aa);
1110 kin_mem->kin_q_aa = NULL;
1111 kin_mem->kin_lrw -= kin_mem->kin_m_aa * kin_mem->kin_lrw1;
1112 kin_mem->kin_liw -= kin_mem->kin_m_aa * kin_mem->kin_liw1;
1113 }
1114
1115 if (kin_mem->kin_constraints != NULL) {
1116 N_VDestroy(kin_mem->kin_constraints);
1117 kin_mem->kin_constraints = NULL;
1118 kin_mem->kin_lrw -= kin_mem->kin_lrw1;
1119 kin_mem->kin_liw -= kin_mem->kin_liw1;
1120 }
1121
1122 return;
1123 }
1124
1125 /*
1126 * -----------------------------------------------------------------
1127 * Initial setup
1128 * -----------------------------------------------------------------
1129 */
1130
1131 /*
1132 * KINSolInit
1133 *
1134 * KINSolInit initializes the problem for the specific input
1135 * received in this call to KINSol (which calls KINSolInit). All
1136 * problem specification inputs are checked for errors. If any error
1137 * occurs during initialization, it is reported to the file whose
1138 * file pointer is errfp.
1139 *
1140 * The possible return values for KINSolInit are:
1141 * KIN_SUCCESS : indicates a normal initialization
1142 *
1143 * KIN_ILL_INPUT : indicates that an input error has been found
1144 *
1145 * KIN_INITIAL_GUESS_OK : indicates that the guess uu
1146 * satisfied the system func(uu) = 0
1147 * within the tolerances specified
1148 */
1149
KINSolInit(KINMem kin_mem)1150 static int KINSolInit(KINMem kin_mem)
1151 {
1152 int retval;
1153 realtype fmax;
1154
1155 /* check for illegal input parameters */
1156
1157 if (kin_mem->kin_uu == NULL) {
1158 KINProcessError(kin_mem, KIN_ILL_INPUT, "KINSOL", "KINSolInit", MSG_UU_NULL);
1159 return(KIN_ILL_INPUT);
1160 }
1161
1162 /* check for valid strategy */
1163
1164 if ( (kin_mem->kin_globalstrategy != KIN_NONE) &&
1165 (kin_mem->kin_globalstrategy != KIN_LINESEARCH) &&
1166 (kin_mem->kin_globalstrategy != KIN_PICARD) &&
1167 (kin_mem->kin_globalstrategy != KIN_FP) ) {
1168 KINProcessError(kin_mem, KIN_ILL_INPUT, "KINSOL", "KINSolInit", MSG_BAD_GLSTRAT);
1169 return(KIN_ILL_INPUT);
1170 }
1171
1172 if (kin_mem->kin_uscale == NULL) {
1173 KINProcessError(kin_mem, KIN_ILL_INPUT, "KINSOL", "KINSolInit", MSG_BAD_USCALE);
1174 return(KIN_ILL_INPUT);
1175 }
1176
1177 if (N_VMin(kin_mem->kin_uscale) <= ZERO){
1178 KINProcessError(kin_mem, KIN_ILL_INPUT, "KINSOL", "KINSolInit", MSG_USCALE_NONPOSITIVE);
1179 return(KIN_ILL_INPUT);
1180 }
1181
1182 if (kin_mem->kin_fscale == NULL) {
1183 KINProcessError(kin_mem, KIN_ILL_INPUT, "KINSOL", "KINSolInit", MSG_BAD_FSCALE);
1184 return(KIN_ILL_INPUT);
1185 }
1186
1187 if (N_VMin(kin_mem->kin_fscale) <= ZERO){
1188 KINProcessError(kin_mem, KIN_ILL_INPUT, "KINSOL", "KINSolInit", MSG_FSCALE_NONPOSITIVE);
1189 return(KIN_ILL_INPUT);
1190 }
1191
1192 if ( (kin_mem->kin_constraints != NULL) &&
1193 ( (kin_mem->kin_globalstrategy == KIN_PICARD) ||
1194 (kin_mem->kin_globalstrategy == KIN_FP) ) ) {
1195 KINProcessError(kin_mem, KIN_ILL_INPUT, "KINSOL", "KINSolInit", MSG_CONSTRAINTS_NOTOK);
1196 return(KIN_ILL_INPUT);
1197 }
1198
1199
1200 /* set the constraints flag */
1201
1202 if (kin_mem->kin_constraints == NULL)
1203 kin_mem->kin_constraintsSet = SUNFALSE;
1204 else {
1205 kin_mem->kin_constraintsSet = SUNTRUE;
1206 if ((kin_mem->kin_constraints->ops->nvconstrmask == NULL) ||
1207 (kin_mem->kin_constraints->ops->nvminquotient == NULL)) {
1208 KINProcessError(kin_mem, KIN_ILL_INPUT, "KINSOL", "KINSolInit", MSG_BAD_NVECTOR);
1209 return(KIN_ILL_INPUT);
1210 }
1211 }
1212
1213 /* check the initial guess uu against the constraints */
1214
1215 if (kin_mem->kin_constraintsSet) {
1216 if (!N_VConstrMask(kin_mem->kin_constraints, kin_mem->kin_uu, kin_mem->kin_vtemp1)) {
1217 KINProcessError(kin_mem, KIN_ILL_INPUT, "KINSOL", "KINSolInit", MSG_INITIAL_CNSTRNT);
1218 return(KIN_ILL_INPUT);
1219 }
1220 }
1221
1222 /* all error checking is complete at this point */
1223
1224 if (kin_mem->kin_printfl > 0)
1225 KINPrintInfo(kin_mem, PRNT_TOL, "KINSOL", "KINSolInit", INFO_TOL, kin_mem->kin_scsteptol, kin_mem->kin_fnormtol);
1226
1227 /* calculate the default value for mxnewtstep (maximum Newton step) */
1228
1229 if (kin_mem->kin_mxnstepin == ZERO) kin_mem->kin_mxnewtstep = THOUSAND * N_VWL2Norm(kin_mem->kin_uu, kin_mem->kin_uscale);
1230 else kin_mem->kin_mxnewtstep = kin_mem->kin_mxnstepin;
1231
1232 if (kin_mem->kin_mxnewtstep < ONE) kin_mem->kin_mxnewtstep = ONE;
1233
1234 /* additional set-up for inexact linear solvers */
1235
1236 if (kin_mem->kin_inexact_ls) {
1237
1238 /* set up the coefficients for the eta calculation */
1239
1240 kin_mem->kin_callForcingTerm = (kin_mem->kin_etaflag != KIN_ETACONSTANT);
1241
1242 /* this value is always used for choice #1 */
1243
1244 if (kin_mem->kin_etaflag == KIN_ETACHOICE1) kin_mem->kin_eta_alpha = (ONE + SUNRsqrt(FIVE)) * HALF;
1245
1246 /* initial value for eta set to 0.5 for other than the
1247 KIN_ETACONSTANT option */
1248
1249 if (kin_mem->kin_etaflag != KIN_ETACONSTANT) kin_mem->kin_eta = HALF;
1250
1251 /* disable residual monitoring if using an inexact linear solver */
1252
1253 kin_mem->kin_noResMon = SUNTRUE;
1254
1255 } else {
1256
1257 kin_mem->kin_callForcingTerm = SUNFALSE;
1258
1259 }
1260
1261 /* initialize counters */
1262
1263 kin_mem->kin_nfe = kin_mem->kin_nnilset = kin_mem->kin_nnilset_sub = kin_mem->kin_nni = kin_mem->kin_nbcf = kin_mem->kin_nbktrk = 0;
1264
1265 /* see if the initial guess uu satisfies the nonlinear system */
1266 retval = kin_mem->kin_func(kin_mem->kin_uu, kin_mem->kin_fval, kin_mem->kin_user_data); kin_mem->kin_nfe++;
1267
1268 if (retval < 0) {
1269 KINProcessError(kin_mem, KIN_SYSFUNC_FAIL, "KINSOL", "KINSolInit",
1270 MSG_SYSFUNC_FAILED);
1271 return(KIN_SYSFUNC_FAIL);
1272 } else if (retval > 0) {
1273 KINProcessError(kin_mem, KIN_FIRST_SYSFUNC_ERR, "KINSOL", "KINSolInit",
1274 MSG_SYSFUNC_FIRST);
1275 return(KIN_FIRST_SYSFUNC_ERR);
1276 }
1277
1278 fmax = KINScFNorm(kin_mem, kin_mem->kin_fval, kin_mem->kin_fscale);
1279 if (fmax <= (POINT01 * kin_mem->kin_fnormtol)) {
1280 kin_mem->kin_fnorm = N_VWL2Norm(kin_mem->kin_fval, kin_mem->kin_fscale);
1281 return(KIN_INITIAL_GUESS_OK);
1282 }
1283
1284 if (kin_mem->kin_printfl > 1)
1285 KINPrintInfo(kin_mem, PRNT_FMAX, "KINSOL", "KINSolInit", INFO_FMAX, fmax);
1286
1287 /* initialize the linear solver if linit != NULL */
1288
1289 if (kin_mem->kin_linit != NULL) {
1290 retval = kin_mem->kin_linit(kin_mem);
1291 if (retval != 0) {
1292 KINProcessError(kin_mem, KIN_LINIT_FAIL, "KINSOL", "KINSolInit", MSG_LINIT_FAIL);
1293 return(KIN_LINIT_FAIL);
1294 }
1295 }
1296
1297 /* initialize the L2 (Euclidean) norms of f for the linear iteration steps */
1298
1299 kin_mem->kin_fnorm = N_VWL2Norm(kin_mem->kin_fval, kin_mem->kin_fscale);
1300 kin_mem->kin_f1norm = HALF * kin_mem->kin_fnorm * kin_mem->kin_fnorm;
1301 kin_mem->kin_fnorm_sub = kin_mem->kin_fnorm;
1302
1303 if (kin_mem->kin_printfl > 0)
1304 KINPrintInfo(kin_mem, PRNT_NNI, "KINSOL", "KINSolInit",
1305 INFO_NNI, kin_mem->kin_nni, kin_mem->kin_nfe, kin_mem->kin_fnorm);
1306
1307 /* problem has now been successfully initialized */
1308
1309 return(KIN_SUCCESS);
1310 }
1311
1312 /*
1313 * -----------------------------------------------------------------
1314 * Step functions
1315 * -----------------------------------------------------------------
1316 */
1317
1318 /*
1319 * KINLinSolDrv
1320 *
1321 * This routine handles the process of solving for the approximate
1322 * solution of the Newton equations in the Newton iteration.
1323 * Subsequent routines handle the nonlinear aspects of its
1324 * application.
1325 */
1326
KINLinSolDrv(KINMem kin_mem)1327 static int KINLinSolDrv(KINMem kin_mem)
1328 {
1329 N_Vector x, b;
1330 int retval;
1331
1332 if ((kin_mem->kin_nni - kin_mem->kin_nnilset) >= kin_mem->kin_msbset) {
1333 kin_mem->kin_sthrsh = TWO;
1334 kin_mem->kin_update_fnorm_sub = SUNTRUE;
1335 }
1336
1337 for(;;){
1338
1339 kin_mem->kin_jacCurrent = SUNFALSE;
1340
1341 if ((kin_mem->kin_sthrsh > ONEPT5) && (kin_mem->kin_lsetup != NULL)) {
1342 retval = kin_mem->kin_lsetup(kin_mem);
1343 kin_mem->kin_jacCurrent = SUNTRUE;
1344 kin_mem->kin_nnilset = kin_mem->kin_nni;
1345 kin_mem->kin_nnilset_sub = kin_mem->kin_nni;
1346 if (retval != 0) return(KIN_LSETUP_FAIL);
1347 }
1348
1349 /* rename vectors for readability */
1350
1351 b = kin_mem->kin_unew;
1352 x = kin_mem->kin_pp;
1353
1354 /* load b with the current value of -fval */
1355
1356 N_VScale(-ONE, kin_mem->kin_fval, b);
1357
1358 /* call the generic 'lsolve' routine to solve the system Jx = b */
1359
1360 retval = kin_mem->kin_lsolve(kin_mem, x, b, &(kin_mem->kin_sJpnorm),
1361 &(kin_mem->kin_sFdotJp));
1362
1363 if (retval == 0) return(KIN_SUCCESS);
1364 else if (retval < 0) return(KIN_LSOLVE_FAIL);
1365 else if ((kin_mem->kin_lsetup == NULL) || (kin_mem->kin_jacCurrent)) return(KIN_LINSOLV_NO_RECOVERY);
1366
1367 /* loop back only if the linear solver setup is in use
1368 and Jacobian information is not current */
1369
1370 kin_mem->kin_sthrsh = TWO;
1371
1372 }
1373 }
1374
1375 /*
1376 * KINFullNewton
1377 *
1378 * This routine is the main driver for the Full Newton
1379 * algorithm. Its purpose is to compute unew = uu + pp in the
1380 * direction pp from uu, taking the full Newton step. The
1381 * step may be constrained if the constraint conditions are
1382 * violated, or if the norm of pp is greater than mxnewtstep.
1383 */
1384
KINFullNewton(KINMem kin_mem,realtype * fnormp,realtype * f1normp,booleantype * maxStepTaken)1385 static int KINFullNewton(KINMem kin_mem, realtype *fnormp, realtype *f1normp,
1386 booleantype *maxStepTaken)
1387 {
1388 realtype pnorm, ratio;
1389 booleantype fOK;
1390 int ircvr, retval;
1391
1392 *maxStepTaken = SUNFALSE;
1393 pnorm = N_VWL2Norm(kin_mem->kin_pp, kin_mem->kin_uscale);
1394 ratio = ONE;
1395 if (pnorm > kin_mem->kin_mxnewtstep) {
1396 ratio = kin_mem->kin_mxnewtstep / pnorm;
1397 N_VScale(ratio, kin_mem->kin_pp, kin_mem->kin_pp);
1398 pnorm = kin_mem->kin_mxnewtstep;
1399 }
1400
1401 if (kin_mem->kin_printfl > 0)
1402 KINPrintInfo(kin_mem, PRNT_PNORM, "KINSOL", "KINFullNewton", INFO_PNORM, pnorm);
1403
1404 /* If constraints are active, then constrain the step accordingly */
1405
1406 kin_mem->kin_stepl = pnorm;
1407 kin_mem->kin_stepmul = ONE;
1408 if (kin_mem->kin_constraintsSet) {
1409 retval = KINConstraint(kin_mem);
1410 if (retval == CONSTR_VIOLATED) {
1411 /* Apply stepmul set in KINConstraint */
1412 ratio *= kin_mem->kin_stepmul;
1413 N_VScale(kin_mem->kin_stepmul, kin_mem->kin_pp, kin_mem->kin_pp);
1414 pnorm *= kin_mem->kin_stepmul;
1415 kin_mem->kin_stepl = pnorm;
1416 if (kin_mem->kin_printfl > 0)
1417 KINPrintInfo(kin_mem, PRNT_PNORM, "KINSOL", "KINFullNewton", INFO_PNORM, pnorm);
1418 if (pnorm <= kin_mem->kin_scsteptol) {
1419 N_VLinearSum(ONE, kin_mem->kin_uu, ONE, kin_mem->kin_pp, kin_mem->kin_unew);
1420 return(STEP_TOO_SMALL);}
1421 }
1422 }
1423
1424 /* Attempt (at most MAX_RECVR times) to evaluate function at the new iterate */
1425
1426 fOK = SUNFALSE;
1427
1428 for (ircvr = 1; ircvr <= MAX_RECVR; ircvr++) {
1429
1430 /* compute the iterate unew = uu + pp */
1431 N_VLinearSum(ONE, kin_mem->kin_uu, ONE, kin_mem->kin_pp, kin_mem->kin_unew);
1432
1433 /* evaluate func(unew) and its norm, and return */
1434 retval = kin_mem->kin_func(kin_mem->kin_unew, kin_mem->kin_fval, kin_mem->kin_user_data); kin_mem->kin_nfe++;
1435
1436 /* if func was successful, accept pp */
1437 if (retval == 0) {fOK = SUNTRUE; break;}
1438
1439 /* if func failed unrecoverably, give up */
1440 else if (retval < 0) return(KIN_SYSFUNC_FAIL);
1441
1442 /* func failed recoverably; cut step in half and try again */
1443 ratio *= HALF;
1444 N_VScale(HALF, kin_mem->kin_pp, kin_mem->kin_pp);
1445 pnorm *= HALF;
1446 kin_mem->kin_stepl = pnorm;
1447 }
1448
1449 /* If func() failed recoverably MAX_RECVR times, give up */
1450
1451 if (!fOK) return(KIN_REPTD_SYSFUNC_ERR);
1452
1453 /* Evaluate function norms */
1454
1455 *fnormp = N_VWL2Norm(kin_mem->kin_fval, kin_mem->kin_fscale);
1456 *f1normp = HALF * (*fnormp) * (*fnormp);
1457
1458 /* scale sFdotJp and sJpnorm by ratio for later use in KINForcingTerm */
1459
1460 kin_mem->kin_sFdotJp *= ratio;
1461 kin_mem->kin_sJpnorm *= ratio;
1462
1463 if (kin_mem->kin_printfl > 1)
1464 KINPrintInfo(kin_mem, PRNT_FNORM, "KINSOL", "KINFullNewton", INFO_FNORM, *fnormp);
1465
1466 if (pnorm > (POINT99 * kin_mem->kin_mxnewtstep)) *maxStepTaken = SUNTRUE;
1467
1468 return(KIN_SUCCESS);
1469 }
1470
1471 /*
1472 * KINLineSearch
1473 *
1474 * The routine KINLineSearch implements the LineSearch algorithm.
1475 * Its purpose is to find unew = uu + rl * pp in the direction pp
1476 * from uu so that:
1477 * t
1478 * func(unew) <= func(uu) + alpha * g (unew - uu) (alpha = 1.e-4)
1479 *
1480 * and
1481 * t
1482 * func(unew) >= func(uu) + beta * g (unew - uu) (beta = 0.9)
1483 *
1484 * where 0 < rlmin <= rl <= rlmax.
1485 *
1486 * Note:
1487 * mxnewtstep
1488 * rlmax = ---------------- if uu+pp is feasible
1489 * ||uscale*pp||_L2
1490 *
1491 * rlmax = 1 otherwise
1492 *
1493 * and
1494 *
1495 * scsteptol
1496 * rlmin = --------------------------
1497 * || pp ||
1498 * || -------------------- ||_L-infinity
1499 * || (1/uscale + SUNRabs(uu)) ||
1500 *
1501 *
1502 * If the system function fails unrecoverably at any time, KINLineSearch
1503 * returns KIN_SYSFUNC_FAIL which will halt the solver.
1504 *
1505 * We attempt to corect recoverable system function failures only before
1506 * the alpha-condition loop; i.e. when the solution is updated with the
1507 * full Newton step (possibly reduced due to constraint violations).
1508 * Once we find a feasible pp, we assume that any update up to pp is
1509 * feasible.
1510 *
1511 * If the step size is limited due to constraint violations and/or
1512 * recoverable system function failures, we set rlmax=1 to ensure
1513 * that the update remains feasible during the attempts to enforce
1514 * the beta-condition (this is not an issue while enforcing the alpha
1515 * condition, as rl can only decrease from 1 at that stage)
1516 */
1517
KINLineSearch(KINMem kin_mem,realtype * fnormp,realtype * f1normp,booleantype * maxStepTaken)1518 static int KINLineSearch(KINMem kin_mem, realtype *fnormp, realtype *f1normp,
1519 booleantype *maxStepTaken)
1520 {
1521 realtype pnorm, ratio, slpi, rlmin, rlength, rl, rlmax, rldiff;
1522 realtype rltmp, rlprev, pt1trl, f1nprv, rllo, rlinc, alpha, beta;
1523 realtype alpha_cond, beta_cond, rl_a, tmp1, rl_b, tmp2, disc;
1524 int ircvr, nbktrk_l, retval;
1525 booleantype firstBacktrack, fOK;
1526
1527 /* Initializations */
1528
1529 nbktrk_l = 0; /* local backtracking counter */
1530 ratio = ONE; /* step change ratio */
1531 alpha = POINT0001;
1532 beta = POINT9;
1533
1534 firstBacktrack = SUNTRUE;
1535 *maxStepTaken = SUNFALSE;
1536
1537 rlprev = f1nprv = ZERO;
1538
1539 /* Compute length of Newton step */
1540
1541 pnorm = N_VWL2Norm(kin_mem->kin_pp, kin_mem->kin_uscale);
1542 rlmax = kin_mem->kin_mxnewtstep / pnorm;
1543 kin_mem->kin_stepl = pnorm;
1544
1545 /* If the full Newton step is too large, set it to the maximum allowable value */
1546
1547 if(pnorm > kin_mem->kin_mxnewtstep ) {
1548 ratio = kin_mem->kin_mxnewtstep / pnorm;
1549 N_VScale(ratio, kin_mem->kin_pp, kin_mem->kin_pp);
1550 pnorm = kin_mem->kin_mxnewtstep;
1551 rlmax = ONE;
1552 kin_mem->kin_stepl = pnorm;
1553 }
1554
1555 /* If constraint checking is activated, check and correct violations */
1556
1557 kin_mem->kin_stepmul = ONE;
1558
1559 if(kin_mem->kin_constraintsSet){
1560 retval = KINConstraint(kin_mem);
1561 if(retval == CONSTR_VIOLATED){
1562 /* Apply stepmul set in KINConstraint */
1563 N_VScale(kin_mem->kin_stepmul, kin_mem->kin_pp, kin_mem->kin_pp);
1564 ratio *= kin_mem->kin_stepmul;
1565 pnorm *= kin_mem->kin_stepmul;
1566 rlmax = ONE;
1567 kin_mem->kin_stepl = pnorm;
1568 if (kin_mem->kin_printfl > 0) KINPrintInfo(kin_mem, PRNT_PNORM1, "KINSOL", "KINLineSearch", INFO_PNORM1, pnorm);
1569 if (pnorm <= kin_mem->kin_scsteptol) {
1570 N_VLinearSum(ONE, kin_mem->kin_uu, ONE, kin_mem->kin_pp, kin_mem->kin_unew);
1571 return(STEP_TOO_SMALL);}
1572 }
1573 }
1574
1575 /* Attempt (at most MAX_RECVR times) to evaluate function at the new iterate */
1576
1577 fOK = SUNFALSE;
1578
1579 for (ircvr = 1; ircvr <= MAX_RECVR; ircvr++) {
1580
1581 /* compute the iterate unew = uu + pp */
1582 N_VLinearSum(ONE, kin_mem->kin_uu, ONE, kin_mem->kin_pp, kin_mem->kin_unew);
1583
1584 /* evaluate func(unew) and its norm, and return */
1585 retval = kin_mem->kin_func(kin_mem->kin_unew, kin_mem->kin_fval, kin_mem->kin_user_data); kin_mem->kin_nfe++;
1586
1587 /* if func was successful, accept pp */
1588 if (retval == 0) {fOK = SUNTRUE; break;}
1589
1590 /* if func failed unrecoverably, give up */
1591 else if (retval < 0) return(KIN_SYSFUNC_FAIL);
1592
1593 /* func failed recoverably; cut step in half and try again */
1594 N_VScale(HALF, kin_mem->kin_pp, kin_mem->kin_pp);
1595 ratio *= HALF;
1596 pnorm *= HALF;
1597 rlmax = ONE;
1598 kin_mem->kin_stepl = pnorm;
1599
1600 }
1601
1602 /* If func() failed recoverably MAX_RECVR times, give up */
1603
1604 if (!fOK) return(KIN_REPTD_SYSFUNC_ERR);
1605
1606 /* Evaluate function norms */
1607
1608 *fnormp = N_VWL2Norm(kin_mem->kin_fval, kin_mem->kin_fscale);
1609 *f1normp = HALF * (*fnormp) * (*fnormp) ;
1610
1611 /* Estimate the line search value rl (lambda) to satisfy both ALPHA and BETA conditions */
1612
1613 slpi = kin_mem->kin_sFdotJp * ratio;
1614 rlength = KINScSNorm(kin_mem, kin_mem->kin_pp, kin_mem->kin_uu);
1615 rlmin = (kin_mem->kin_scsteptol) / rlength;
1616 rl = ONE;
1617
1618 if (kin_mem->kin_printfl > 2)
1619 KINPrintInfo(kin_mem, PRNT_LAM, "KINSOL", "KINLineSearch", INFO_LAM, rlmin, kin_mem->kin_f1norm, pnorm);
1620
1621 /* Loop until the ALPHA condition is satisfied. Terminate if rl becomes too small */
1622
1623 for(;;) {
1624
1625 /* Evaluate test quantity */
1626
1627 alpha_cond = kin_mem->kin_f1norm + (alpha * slpi * rl);
1628
1629 if (kin_mem->kin_printfl > 2)
1630 KINPrintInfo(kin_mem, PRNT_ALPHA, "KINSOL", "KINLinesearch",
1631 INFO_ALPHA, *fnormp, *f1normp, alpha_cond, rl);
1632
1633 /* If ALPHA condition is satisfied, break out from loop */
1634
1635 if ((*f1normp) <= alpha_cond) break;
1636
1637 /* Backtracking. Use quadratic fit the first time and cubic fit afterwards. */
1638
1639 if (firstBacktrack) {
1640
1641 rltmp = -slpi / (TWO * ((*f1normp) - kin_mem->kin_f1norm - slpi));
1642 firstBacktrack = SUNFALSE;
1643
1644 } else {
1645
1646 tmp1 = (*f1normp) - kin_mem->kin_f1norm - (rl * slpi);
1647 tmp2 = f1nprv - kin_mem->kin_f1norm - (rlprev * slpi);
1648 rl_a = ((ONE / (rl * rl)) * tmp1) - ((ONE / (rlprev * rlprev)) * tmp2);
1649 rl_b = ((-rlprev / (rl * rl)) * tmp1) + ((rl / (rlprev * rlprev)) * tmp2);
1650 tmp1 = ONE / (rl - rlprev);
1651 rl_a *= tmp1;
1652 rl_b *= tmp1;
1653 disc = (rl_b * rl_b) - (THREE * rl_a * slpi);
1654
1655 if (SUNRabs(rl_a) < kin_mem->kin_uround) { /* cubic is actually just a quadratic (rl_a ~ 0) */
1656 rltmp = -slpi / (TWO * rl_b);
1657 } else { /* real cubic */
1658 rltmp = (-rl_b + SUNRsqrt(disc)) / (THREE * rl_a);
1659 }
1660 }
1661 if (rltmp > (HALF * rl)) rltmp = HALF * rl;
1662
1663 /* Set new rl (do not allow a reduction by a factor larger than 10) */
1664
1665 rlprev = rl;
1666 f1nprv = (*f1normp);
1667 pt1trl = POINT1 * rl;
1668 rl = SUNMAX(pt1trl, rltmp);
1669 nbktrk_l++;
1670
1671 /* Update unew and re-evaluate function */
1672
1673 N_VLinearSum(ONE, kin_mem->kin_uu, rl, kin_mem->kin_pp, kin_mem->kin_unew);
1674
1675 retval = kin_mem->kin_func(kin_mem->kin_unew, kin_mem->kin_fval, kin_mem->kin_user_data); kin_mem->kin_nfe++;
1676 if (retval != 0) return(KIN_SYSFUNC_FAIL);
1677
1678 *fnormp = N_VWL2Norm(kin_mem->kin_fval, kin_mem->kin_fscale);
1679 *f1normp = HALF * (*fnormp) * (*fnormp) ;
1680
1681 /* Check if rl (lambda) is too small */
1682
1683 if (rl < rlmin) {
1684 /* unew sufficiently distinct from uu cannot be found.
1685 copy uu into unew (step remains unchanged) and
1686 return STEP_TOO_SMALL */
1687 N_VScale(ONE, kin_mem->kin_uu, kin_mem->kin_unew);
1688 return(STEP_TOO_SMALL);
1689 }
1690
1691 } /* end ALPHA condition loop */
1692
1693
1694 /* ALPHA condition is satisfied. Now check the BETA condition */
1695
1696 beta_cond = kin_mem->kin_f1norm + (beta * slpi * rl);
1697
1698 if ((*f1normp) < beta_cond) {
1699
1700 /* BETA condition not satisfied */
1701
1702 if ((rl == ONE) && (pnorm < kin_mem->kin_mxnewtstep)) {
1703
1704 do {
1705
1706 rlprev = rl;
1707 f1nprv = *f1normp;
1708 rl = SUNMIN((TWO * rl), rlmax);
1709 nbktrk_l++;
1710
1711 N_VLinearSum(ONE, kin_mem->kin_uu, rl, kin_mem->kin_pp, kin_mem->kin_unew);
1712 retval = kin_mem->kin_func(kin_mem->kin_unew, kin_mem->kin_fval, kin_mem->kin_user_data); kin_mem->kin_nfe++;
1713 if (retval != 0) return(KIN_SYSFUNC_FAIL);
1714 *fnormp = N_VWL2Norm(kin_mem->kin_fval, kin_mem->kin_fscale);
1715 *f1normp = HALF * (*fnormp) * (*fnormp);
1716
1717 alpha_cond = kin_mem->kin_f1norm + (alpha * slpi * rl);
1718 beta_cond = kin_mem->kin_f1norm + (beta * slpi * rl);
1719
1720 if (kin_mem->kin_printfl > 2)
1721 KINPrintInfo(kin_mem, PRNT_BETA, "KINSOL", "KINLineSearch",
1722 INFO_BETA, *f1normp, beta_cond, rl);
1723
1724 } while (((*f1normp) <= alpha_cond) &&
1725 ((*f1normp) < beta_cond) && (rl < rlmax));
1726
1727 } /* end if (rl == ONE) block */
1728
1729 if ((rl < ONE) || ((rl > ONE) && (*f1normp > alpha_cond))) {
1730
1731 rllo = SUNMIN(rl, rlprev);
1732 rldiff = SUNRabs(rlprev - rl);
1733
1734 do {
1735
1736 rlinc = HALF * rldiff;
1737 rl = rllo + rlinc;
1738 nbktrk_l++;
1739
1740 N_VLinearSum(ONE, kin_mem->kin_uu, rl, kin_mem->kin_pp, kin_mem->kin_unew);
1741 retval = kin_mem->kin_func(kin_mem->kin_unew, kin_mem->kin_fval, kin_mem->kin_user_data); kin_mem->kin_nfe++;
1742 if (retval != 0) return(KIN_SYSFUNC_FAIL);
1743 *fnormp = N_VWL2Norm(kin_mem->kin_fval, kin_mem->kin_fscale);
1744 *f1normp = HALF * (*fnormp) * (*fnormp);
1745
1746 alpha_cond = kin_mem->kin_f1norm + (alpha * slpi * rl);
1747 beta_cond = kin_mem->kin_f1norm + (beta * slpi * rl);
1748
1749 if (kin_mem->kin_printfl > 2)
1750 KINPrintInfo(kin_mem, PRNT_ALPHABETA, "KINSOL", "KINLineSearch",
1751 INFO_ALPHABETA, *f1normp, alpha_cond, beta_cond, rl);
1752
1753 if ((*f1normp) > alpha_cond) rldiff = rlinc;
1754 else if (*f1normp < beta_cond) {
1755 rllo = rl;
1756 rldiff = rldiff - rlinc;
1757 }
1758
1759 } while ((*f1normp > alpha_cond) ||
1760 ((*f1normp < beta_cond) && (rldiff >= rlmin)));
1761
1762 if ( (*f1normp < beta_cond) || ((rldiff < rlmin) && (*f1normp > alpha_cond)) ) {
1763
1764 /* beta condition could not be satisfied or rldiff too small
1765 and alpha_cond not satisfied, so set unew to last u value
1766 that satisfied the alpha condition and continue */
1767
1768 N_VLinearSum(ONE, kin_mem->kin_uu, rllo, kin_mem->kin_pp, kin_mem->kin_unew);
1769 retval = kin_mem->kin_func(kin_mem->kin_unew, kin_mem->kin_fval, kin_mem->kin_user_data); kin_mem->kin_nfe++;
1770 if (retval != 0) return(KIN_SYSFUNC_FAIL);
1771 *fnormp = N_VWL2Norm(kin_mem->kin_fval, kin_mem->kin_fscale);
1772 *f1normp = HALF * (*fnormp) * (*fnormp);
1773
1774 /* increment beta-condition failures counter */
1775
1776 kin_mem->kin_nbcf++;
1777
1778 }
1779
1780 } /* end of if (rl < ONE) block */
1781
1782 } /* end of if (f1normp < beta_cond) block */
1783
1784 /* Update number of backtracking operations */
1785
1786 kin_mem->kin_nbktrk += nbktrk_l;
1787
1788 if (kin_mem->kin_printfl > 1)
1789 KINPrintInfo(kin_mem, PRNT_ADJ, "KINSOL", "KINLineSearch", INFO_ADJ, nbktrk_l);
1790
1791 /* scale sFdotJp and sJpnorm by rl * ratio for later use in KINForcingTerm */
1792
1793 kin_mem->kin_sFdotJp = kin_mem->kin_sFdotJp * rl * ratio;
1794 kin_mem->kin_sJpnorm = kin_mem->kin_sJpnorm * rl * ratio;
1795
1796 if ((rl * pnorm) > (POINT99 * kin_mem->kin_mxnewtstep)) *maxStepTaken = SUNTRUE;
1797
1798 return(KIN_SUCCESS);
1799 }
1800
1801 /*
1802 * Function : KINConstraint
1803 *
1804 * This routine checks if the proposed solution vector uu + pp
1805 * violates any constraints. If a constraint is violated, then the
1806 * scalar stepmul is determined such that uu + stepmul * pp does
1807 * not violate any constraints.
1808 *
1809 * Note: This routine is called by the functions
1810 * KINLineSearch and KINFullNewton.
1811 */
1812
KINConstraint(KINMem kin_mem)1813 static int KINConstraint(KINMem kin_mem)
1814 {
1815 N_VLinearSum(ONE, kin_mem->kin_uu, ONE, kin_mem->kin_pp, kin_mem->kin_vtemp1);
1816
1817 /* if vtemp1[i] violates constraint[i] then vtemp2[i] = 1
1818 else vtemp2[i] = 0 (vtemp2 is the mask vector) */
1819
1820 if(N_VConstrMask(kin_mem->kin_constraints, kin_mem->kin_vtemp1, kin_mem->kin_vtemp2)) return(KIN_SUCCESS);
1821
1822 /* vtemp1[i] = SUNRabs(pp[i]) */
1823
1824 N_VAbs(kin_mem->kin_pp, kin_mem->kin_vtemp1);
1825
1826 /* consider vtemp1[i] only if vtemp2[i] = 1 (constraint violated) */
1827
1828 N_VProd(kin_mem->kin_vtemp2, kin_mem->kin_vtemp1, kin_mem->kin_vtemp1);
1829
1830 N_VAbs(kin_mem->kin_uu, kin_mem->kin_vtemp2);
1831 kin_mem->kin_stepmul = POINT9 * N_VMinQuotient(kin_mem->kin_vtemp2, kin_mem->kin_vtemp1);
1832
1833 return(CONSTR_VIOLATED);
1834 }
1835
1836 /*
1837 * -----------------------------------------------------------------
1838 * Stopping tests
1839 * -----------------------------------------------------------------
1840 */
1841
1842 /*
1843 * KINStop
1844 *
1845 * This routine checks the current iterate unew to see if the
1846 * system func(unew) = 0 is satisfied by a variety of tests.
1847 *
1848 * strategy is one of KIN_NONE or KIN_LINESEARCH
1849 * sflag is one of KIN_SUCCESS, STEP_TOO_SMALL
1850 */
1851
KINStop(KINMem kin_mem,booleantype maxStepTaken,int sflag)1852 static int KINStop(KINMem kin_mem, booleantype maxStepTaken, int sflag)
1853 {
1854 realtype fmax, rlength, omexp;
1855 N_Vector delta;
1856
1857 /* Check for too small a step */
1858
1859 if (sflag == STEP_TOO_SMALL) {
1860
1861 if ((kin_mem->kin_lsetup != NULL) && !(kin_mem->kin_jacCurrent)) {
1862 /* If the Jacobian is out of date, update it and retry */
1863 kin_mem->kin_sthrsh = TWO;
1864 return(RETRY_ITERATION);
1865 } else {
1866 /* Give up */
1867 if (kin_mem->kin_globalstrategy == KIN_NONE) return(KIN_STEP_LT_STPTOL);
1868 else return(KIN_LINESEARCH_NONCONV);
1869 }
1870
1871 }
1872
1873 /* Check tolerance on scaled function norm at the current iterate */
1874
1875 fmax = KINScFNorm(kin_mem, kin_mem->kin_fval, kin_mem->kin_fscale);
1876
1877 if (kin_mem->kin_printfl > 1)
1878 KINPrintInfo(kin_mem, PRNT_FMAX, "KINSOL", "KINStop", INFO_FMAX, fmax);
1879
1880 if (fmax <= kin_mem->kin_fnormtol) return(KIN_SUCCESS);
1881
1882 /* Check if the scaled distance between the last two steps is too small */
1883 /* NOTE: pp used as work space to store this distance */
1884
1885 delta = kin_mem->kin_pp;
1886 N_VLinearSum(ONE, kin_mem->kin_unew, -ONE, kin_mem->kin_uu, delta);
1887 rlength = KINScSNorm(kin_mem, delta, kin_mem->kin_unew);
1888
1889 if (rlength <= kin_mem->kin_scsteptol) {
1890
1891 if ((kin_mem->kin_lsetup != NULL) && !(kin_mem->kin_jacCurrent)) {
1892 /* If the Jacobian is out of date, update it and retry */
1893 kin_mem->kin_sthrsh = TWO;
1894 return(CONTINUE_ITERATIONS);
1895 } else {
1896 /* give up */
1897 return(KIN_STEP_LT_STPTOL);
1898 }
1899
1900 }
1901
1902 /* Check if the maximum number of iterations is reached */
1903
1904 if (kin_mem->kin_nni >= kin_mem->kin_mxiter) return(KIN_MAXITER_REACHED);
1905
1906 /* Check for consecutive number of steps taken of size mxnewtstep
1907 and if not maxStepTaken, then set ncscmx to 0 */
1908
1909 if (maxStepTaken) kin_mem->kin_ncscmx++;
1910 else kin_mem->kin_ncscmx = 0;
1911
1912 if (kin_mem->kin_ncscmx == 5) return(KIN_MXNEWT_5X_EXCEEDED);
1913
1914 /* Proceed according to the type of linear solver used */
1915
1916 if (kin_mem->kin_inexact_ls) {
1917
1918 /* We're doing inexact Newton.
1919 Load threshold for reevaluating the Jacobian. */
1920
1921 kin_mem->kin_sthrsh = rlength;
1922
1923 } else if (!(kin_mem->kin_noResMon)) {
1924
1925 /* We're doing modified Newton and the user did not disable residual monitoring.
1926 Check if it is time to monitor residual. */
1927
1928 if ((kin_mem->kin_nni - kin_mem->kin_nnilset_sub) >= kin_mem->kin_msbset_sub) {
1929
1930 /* Residual monitoring needed */
1931
1932 kin_mem->kin_nnilset_sub = kin_mem->kin_nni;
1933
1934 /* If indicated, estimate new OMEGA value */
1935 if (kin_mem->kin_eval_omega) {
1936 omexp = SUNMAX(ZERO,((kin_mem->kin_fnorm)/(kin_mem->kin_fnormtol))-ONE);
1937 kin_mem->kin_omega = (omexp > TWELVE)? kin_mem->kin_omega_max : SUNMIN(kin_mem->kin_omega_min * SUNRexp(omexp), kin_mem->kin_omega_max);
1938 }
1939 /* Check if making satisfactory progress */
1940
1941 if (kin_mem->kin_fnorm > kin_mem->kin_omega * kin_mem->kin_fnorm_sub) {
1942 /* Insufficient progress */
1943 if ((kin_mem->kin_lsetup != NULL) && !(kin_mem->kin_jacCurrent)) {
1944 /* If the Jacobian is out of date, update it and retry */
1945 kin_mem->kin_sthrsh = TWO;
1946 return(CONTINUE_ITERATIONS);
1947 } else {
1948 /* Otherwise, we cannot do anything, so just return. */
1949 }
1950 } else {
1951 /* Sufficient progress */
1952 kin_mem->kin_fnorm_sub = kin_mem->kin_fnorm;
1953 kin_mem->kin_sthrsh = ONE;
1954 }
1955
1956 } else {
1957
1958 /* Residual monitoring not needed */
1959
1960 /* Reset sthrsh */
1961 if (kin_mem->kin_retry_nni || kin_mem->kin_update_fnorm_sub) kin_mem->kin_fnorm_sub = kin_mem->kin_fnorm;
1962 if (kin_mem->kin_update_fnorm_sub) kin_mem->kin_update_fnorm_sub = SUNFALSE;
1963 kin_mem->kin_sthrsh = ONE;
1964
1965 }
1966
1967 }
1968
1969 /* if made it to here, then the iteration process is not finished
1970 so return CONTINUE_ITERATIONS flag */
1971
1972 return(CONTINUE_ITERATIONS);
1973 }
1974
1975 /*
1976 * KINForcingTerm
1977 *
1978 * This routine computes eta, the scaling factor in the linear
1979 * convergence stopping tolerance eps when choice #1 or choice #2
1980 * forcing terms are used. Eta is computed here for all but the
1981 * first iterative step, which is set to the default in routine
1982 * KINSolInit.
1983 *
1984 * This routine was written by Homer Walker of Utah State
1985 * University with subsequent modifications by Allan Taylor @ LLNL.
1986 *
1987 * It is based on the concepts of the paper 'Choosing the forcing
1988 * terms in an inexact Newton method', SIAM J Sci Comput, 17
1989 * (1996), pp 16 - 32, or Utah State University Research Report
1990 * 6/94/75 of the same title.
1991 */
1992
KINForcingTerm(KINMem kin_mem,realtype fnormp)1993 static void KINForcingTerm(KINMem kin_mem, realtype fnormp)
1994 {
1995 realtype eta_max, eta_min, eta_safe, linmodel_norm;
1996
1997 eta_max = POINT9;
1998 eta_min = POINT0001;
1999 eta_safe = HALF;
2000
2001 /* choice #1 forcing term */
2002
2003 if (kin_mem->kin_etaflag == KIN_ETACHOICE1) {
2004
2005 /* compute the norm of f + Jp , scaled L2 norm */
2006
2007 linmodel_norm = SUNRsqrt((kin_mem->kin_fnorm * kin_mem->kin_fnorm) +
2008 (TWO * kin_mem->kin_sFdotJp) +
2009 (kin_mem->kin_sJpnorm * kin_mem->kin_sJpnorm));
2010
2011 /* form the safeguarded for choice #1 */
2012
2013 eta_safe = SUNRpowerR(kin_mem->kin_eta, kin_mem->kin_eta_alpha);
2014 kin_mem->kin_eta = SUNRabs(fnormp - linmodel_norm) / kin_mem->kin_fnorm;
2015 }
2016
2017 /* choice #2 forcing term */
2018
2019 if (kin_mem->kin_etaflag == KIN_ETACHOICE2) {
2020 eta_safe = kin_mem->kin_eta_gamma *
2021 SUNRpowerR(kin_mem->kin_eta, kin_mem->kin_eta_alpha);
2022
2023 kin_mem->kin_eta = kin_mem->kin_eta_gamma *
2024 SUNRpowerR((fnormp / kin_mem->kin_fnorm), kin_mem->kin_eta_alpha);
2025 }
2026
2027 /* apply safeguards */
2028
2029 if(eta_safe < POINT1) eta_safe = ZERO;
2030 kin_mem->kin_eta = SUNMAX(kin_mem->kin_eta, eta_safe);
2031 kin_mem->kin_eta = SUNMAX(kin_mem->kin_eta, eta_min);
2032 kin_mem->kin_eta = SUNMIN(kin_mem->kin_eta, eta_max);
2033
2034 return;
2035 }
2036
2037
2038 /*
2039 * -----------------------------------------------------------------
2040 * Norm functions
2041 * -----------------------------------------------------------------
2042 */
2043
2044 /*
2045 * Function : KINScFNorm
2046 *
2047 * This routine computes the max norm for scaled vectors. The
2048 * scaling vector is scale, and the vector of which the norm is to
2049 * be determined is vv. The returned value, fnormval, is the
2050 * resulting scaled vector norm. This routine uses N_Vector
2051 * functions from the vector module.
2052 */
2053
KINScFNorm(KINMem kin_mem,N_Vector v,N_Vector scale)2054 static realtype KINScFNorm(KINMem kin_mem, N_Vector v, N_Vector scale)
2055 {
2056 N_VProd(scale, v, kin_mem->kin_vtemp1);
2057 return(N_VMaxNorm(kin_mem->kin_vtemp1));
2058 }
2059
2060 /*
2061 * Function : KINScSNorm
2062 *
2063 * This routine computes the max norm of the scaled steplength, ss.
2064 * Here ucur is the current step and usc is the u scale factor.
2065 */
2066
KINScSNorm(KINMem kin_mem,N_Vector v,N_Vector u)2067 static realtype KINScSNorm(KINMem kin_mem, N_Vector v, N_Vector u)
2068 {
2069 realtype length;
2070
2071 N_VInv(kin_mem->kin_uscale, kin_mem->kin_vtemp1);
2072 N_VAbs(u, kin_mem->kin_vtemp2);
2073 N_VLinearSum(ONE, kin_mem->kin_vtemp1, ONE, kin_mem->kin_vtemp2, kin_mem->kin_vtemp1);
2074 N_VDiv(v, kin_mem->kin_vtemp1, kin_mem->kin_vtemp1);
2075
2076 length = N_VMaxNorm(kin_mem->kin_vtemp1);
2077
2078 return(length);
2079 }
2080
2081 /*
2082 * =================================================================
2083 * KINSOL Verbose output functions
2084 * =================================================================
2085 */
2086
2087 /*
2088 * KINPrintInfo
2089 *
2090 * KINPrintInfo is a high level error handling function
2091 * Based on the value info_code, it composes the info message and
2092 * passes it to the info handler function.
2093 */
2094
KINPrintInfo(KINMem kin_mem,int info_code,const char * module,const char * fname,const char * msgfmt,...)2095 void KINPrintInfo(KINMem kin_mem,
2096 int info_code, const char *module, const char *fname,
2097 const char *msgfmt, ...)
2098 {
2099 va_list ap;
2100 char msg[256], msg1[40];
2101 char retstr[30];
2102 int ret;
2103
2104 /* Initialize argument processing
2105 (msgfrmt is the last required argument) */
2106
2107 va_start(ap, msgfmt);
2108
2109 if (info_code == PRNT_RETVAL) {
2110
2111 /* If info_code = PRNT_RETVAL, decode the numeric value */
2112
2113 ret = va_arg(ap, int);
2114
2115 switch(ret) {
2116 case KIN_SUCCESS:
2117 sprintf(retstr, "KIN_SUCCESS");
2118 break;
2119 case KIN_SYSFUNC_FAIL:
2120 sprintf(retstr, "KIN_SYSFUNC_FAIL");
2121 break;
2122 case KIN_REPTD_SYSFUNC_ERR:
2123 sprintf(retstr, "KIN_REPTD_SYSFUNC_ERR");
2124 break;
2125 case KIN_STEP_LT_STPTOL:
2126 sprintf(retstr, "KIN_STEP_LT_STPTOL");
2127 break;
2128 case KIN_LINESEARCH_NONCONV:
2129 sprintf(retstr, "KIN_LINESEARCH_NONCONV");
2130 break;
2131 case KIN_LINESEARCH_BCFAIL:
2132 sprintf(retstr, "KIN_LINESEARCH_BCFAIL");
2133 break;
2134 case KIN_MAXITER_REACHED:
2135 sprintf(retstr, "KIN_MAXITER_REACHED");
2136 break;
2137 case KIN_MXNEWT_5X_EXCEEDED:
2138 sprintf(retstr, "KIN_MXNEWT_5X_EXCEEDED");
2139 break;
2140 case KIN_LINSOLV_NO_RECOVERY:
2141 sprintf(retstr, "KIN_LINSOLV_NO_RECOVERY");
2142 break;
2143 case KIN_LSETUP_FAIL:
2144 sprintf(retstr, "KIN_PRECONDSET_FAILURE");
2145 break;
2146 case KIN_LSOLVE_FAIL:
2147 sprintf(retstr, "KIN_PRECONDSOLVE_FAILURE");
2148 break;
2149 }
2150
2151 /* Compose the message */
2152
2153 sprintf(msg1, msgfmt, ret);
2154 sprintf(msg,"%s (%s)",msg1,retstr);
2155
2156
2157 } else {
2158
2159 /* Compose the message */
2160
2161 vsprintf(msg, msgfmt, ap);
2162
2163 }
2164
2165 /* call the info message handler */
2166
2167 kin_mem->kin_ihfun(module, fname, msg, kin_mem->kin_ih_data);
2168
2169 /* finalize argument processing */
2170
2171 va_end(ap);
2172
2173 return;
2174 }
2175
2176
2177 /*
2178 * KINInfoHandler
2179 *
2180 * This is the default KINSOL info handling function.
2181 * It sends the info message to the stream pointed to by kin_infofp
2182 */
2183
KINInfoHandler(const char * module,const char * function,char * msg,void * data)2184 void KINInfoHandler(const char *module, const char *function,
2185 char *msg, void *data)
2186 {
2187 KINMem kin_mem;
2188
2189 /* data points to kin_mem here */
2190
2191 kin_mem = (KINMem) data;
2192
2193 #ifndef NO_FPRINTF_OUTPUT
2194 if (kin_mem->kin_infofp != NULL) {
2195 STAN_SUNDIALS_FPRINTF(kin_mem->kin_infofp,"\n[%s] %s\n",module, function);
2196 STAN_SUNDIALS_FPRINTF(kin_mem->kin_infofp," %s\n",msg);
2197 }
2198 #endif
2199
2200 }
2201
2202 /*
2203 * =================================================================
2204 * KINSOL Error Handling functions
2205 * =================================================================
2206 */
2207
2208 /*
2209 * KINProcessError
2210 *
2211 * KINProcessError is a high level error handling function.
2212 * - If cv_mem==NULL it prints the error message to stderr.
2213 * - Otherwise, it sets up and calls the error handling function
2214 * pointed to by cv_ehfun.
2215 */
2216
KINProcessError(KINMem kin_mem,int error_code,const char * module,const char * fname,const char * msgfmt,...)2217 void KINProcessError(KINMem kin_mem,
2218 int error_code, const char *module, const char *fname,
2219 const char *msgfmt, ...)
2220 {
2221 va_list ap;
2222 char msg[256];
2223
2224 /* Initialize the argument pointer variable
2225 (msgfmt is the last required argument to KINProcessError) */
2226
2227 va_start(ap, msgfmt);
2228
2229 /* Compose the message */
2230
2231 vsprintf(msg, msgfmt, ap);
2232
2233 if (kin_mem == NULL) { /* We write to stderr */
2234 #ifndef NO_FPRINTF_OUTPUT
2235 STAN_SUNDIALS_FPRINTF(stderr, "\n[%s ERROR] %s\n ", module, fname);
2236 STAN_SUNDIALS_FPRINTF(stderr, "%s\n\n", msg);
2237 #endif
2238
2239 } else { /* We can call ehfun */
2240 kin_mem->kin_ehfun(error_code, module, fname, msg, kin_mem->kin_eh_data);
2241 }
2242
2243 /* Finalize argument processing */
2244 va_end(ap);
2245
2246 return;
2247 }
2248
2249 /*
2250 * KINErrHandler
2251 *
2252 * This is the default error handling function.
2253 * It sends the error message to the stream pointed to by kin_errfp
2254 */
2255
KINErrHandler(int error_code,const char * module,const char * function,char * msg,void * data)2256 void KINErrHandler(int error_code, const char *module,
2257 const char *function, char *msg, void *data)
2258 {
2259 KINMem kin_mem;
2260 char err_type[10];
2261
2262 /* data points to kin_mem here */
2263
2264 kin_mem = (KINMem) data;
2265
2266 if (error_code == KIN_WARNING)
2267 sprintf(err_type,"WARNING");
2268 else
2269 sprintf(err_type,"ERROR");
2270
2271 #ifndef NO_FPRINTF_OUTPUT
2272 if (kin_mem->kin_errfp != NULL) {
2273 STAN_SUNDIALS_FPRINTF(kin_mem->kin_errfp,"\n[%s %s] %s\n",module,err_type,function);
2274 STAN_SUNDIALS_FPRINTF(kin_mem->kin_errfp," %s\n\n",msg);
2275 }
2276 #endif
2277
2278 return;
2279 }
2280
2281
2282 /*
2283 * =======================================================================
2284 * Picard and fixed point solvers
2285 * =======================================================================
2286 */
2287
2288 /*
2289 * KINPicardAA
2290 *
2291 * This routine is the main driver for the Picard iteration with
2292 * acclerated fixed point.
2293 */
2294
KINPicardAA(KINMem kin_mem,long int * iterp,realtype * R,realtype * gamma,realtype * fmaxptr)2295 static int KINPicardAA(KINMem kin_mem, long int *iterp, realtype *R,
2296 realtype *gamma, realtype *fmaxptr)
2297 {
2298 int retval, ret;
2299 long int iter;
2300 realtype fmax, epsmin, fnormp;
2301 N_Vector delta, gval;
2302
2303 delta = kin_mem->kin_vtemp1;
2304 gval = kin_mem->kin_gval;
2305 ret = CONTINUE_ITERATIONS;
2306 fmax = kin_mem->kin_fnormtol + ONE;
2307 iter = 0;
2308 epsmin = ZERO;
2309 fnormp = -ONE;
2310
2311 N_VConst(ZERO, gval);
2312
2313 /* if eps is to be bounded from below, set the bound */
2314 if (kin_mem->kin_inexact_ls && !(kin_mem->kin_noMinEps)) epsmin = POINT01 * kin_mem->kin_fnormtol;
2315
2316 while (ret == CONTINUE_ITERATIONS) {
2317
2318 iter++;
2319
2320 /* Update the forcing term for the inexact linear solves */
2321 if (kin_mem->kin_inexact_ls) {
2322 kin_mem->kin_eps = (kin_mem->kin_eta + kin_mem->kin_uround) * kin_mem->kin_fnorm;
2323 if(!(kin_mem->kin_noMinEps)) kin_mem->kin_eps = SUNMAX(epsmin, kin_mem->kin_eps);
2324 }
2325
2326 /* evaluate g = uu - L^{-1}func(uu) and return if failed.
2327 For Picard, assume that the fval vector has been filled
2328 with an eval of the nonlinear residual prior to this call. */
2329 retval = KINPicardFcnEval(kin_mem, gval, kin_mem->kin_uu, kin_mem->kin_fval);
2330
2331 if (retval < 0) {
2332 ret = KIN_SYSFUNC_FAIL;
2333 break;
2334 }
2335
2336 if (kin_mem->kin_m_aa == 0) {
2337 N_VScale(ONE, gval, kin_mem->kin_unew);
2338 }
2339 else { /* use Anderson, if desired */
2340 N_VScale(ONE, kin_mem->kin_uu, kin_mem->kin_unew);
2341 AndersonAcc(kin_mem, gval, delta, kin_mem->kin_unew, kin_mem->kin_uu, iter-1, R, gamma);
2342 }
2343
2344 /* Fill the Newton residual based on the new solution iterate */
2345 retval = kin_mem->kin_func(kin_mem->kin_unew, kin_mem->kin_fval, kin_mem->kin_user_data); kin_mem->kin_nfe++;
2346
2347 if (retval < 0) {
2348 ret = KIN_SYSFUNC_FAIL;
2349 break;
2350 }
2351
2352 /* Evaluate function norms */
2353 fnormp = N_VWL2Norm(kin_mem->kin_fval, kin_mem->kin_fscale);
2354 fmax = KINScFNorm(kin_mem, kin_mem->kin_fval, kin_mem->kin_fscale); /* measure || F(x) ||_max */
2355 kin_mem->kin_fnorm = fmax;
2356 *fmaxptr = fmax;
2357
2358 if (kin_mem->kin_printfl > 1)
2359 KINPrintInfo(kin_mem, PRNT_FMAX, "KINSOL", "KINPicardAA", INFO_FMAX, fmax);
2360
2361 /* print the current iter, fnorm, and nfe values if printfl > 0 */
2362 if (kin_mem->kin_printfl > 0)
2363 KINPrintInfo(kin_mem, PRNT_NNI, "KINSOL", "KINPicardAA", INFO_NNI, iter, kin_mem->kin_nfe, kin_mem->kin_fnorm);
2364
2365 /* Check if the maximum number of iterations is reached */
2366 if (iter >= kin_mem->kin_mxiter) {
2367 ret = KIN_MAXITER_REACHED;
2368 }
2369 if (fmax <= kin_mem->kin_fnormtol) {
2370 ret = KIN_SUCCESS;
2371 }
2372
2373 /* Update with new iterate. */
2374 N_VScale(ONE, kin_mem->kin_unew, kin_mem->kin_uu);
2375
2376 if (ret == CONTINUE_ITERATIONS) {
2377 /* evaluate eta by calling the forcing term routine */
2378 if (kin_mem->kin_callForcingTerm) KINForcingTerm(kin_mem, fnormp);
2379 }
2380
2381 fflush(kin_mem->kin_errfp);
2382
2383 } /* end of loop; return */
2384
2385 *iterp = iter;
2386
2387 if (kin_mem->kin_printfl > 0)
2388 KINPrintInfo(kin_mem, PRNT_RETVAL, "KINSOL", "KINPicardAA", INFO_RETVAL, ret);
2389
2390 return(ret);
2391 }
2392
2393 /*
2394 * KINPicardFcnEval
2395 *
2396 * This routine evaluates the Picard fixed point function
2397 * using the linear solver, gval = u - L^{-1}F(u).
2398 * The function assumes the user has defined L either through
2399 * a user-supplied matvec if using a SPILS solver or through
2400 * a supplied matrix if using a dense solver. This assumption is
2401 * tested by a check on the strategy and the requisite functionality
2402 * within the linear solve routines.
2403 *
2404 * This routine fills gval = uu - L^{-1}F(uu) given uu and fval = F(uu).
2405 */
2406
KINPicardFcnEval(KINMem kin_mem,N_Vector gval,N_Vector uval,N_Vector fval1)2407 static int KINPicardFcnEval(KINMem kin_mem, N_Vector gval, N_Vector uval, N_Vector fval1)
2408 {
2409 int retval;
2410
2411 if ((kin_mem->kin_nni - kin_mem->kin_nnilset) >= kin_mem->kin_msbset) {
2412 kin_mem->kin_sthrsh = TWO;
2413 kin_mem->kin_update_fnorm_sub = SUNTRUE;
2414 }
2415
2416 for(;;){
2417
2418 kin_mem->kin_jacCurrent = SUNFALSE;
2419
2420 if ((kin_mem->kin_sthrsh > ONEPT5) && (kin_mem->kin_lsetup != NULL)) {
2421 retval = kin_mem->kin_lsetup(kin_mem);
2422 kin_mem->kin_jacCurrent = SUNTRUE;
2423 kin_mem->kin_nnilset = kin_mem->kin_nni;
2424 kin_mem->kin_nnilset_sub = kin_mem->kin_nni;
2425 if (retval != 0) return(KIN_LSETUP_FAIL);
2426 }
2427
2428 /* call the generic 'lsolve' routine to solve the system Lx = -fval
2429 Note that we are using gval to hold x. */
2430 N_VScale(-ONE, fval1, fval1);
2431 retval = kin_mem->kin_lsolve(kin_mem, gval, fval1, &(kin_mem->kin_sJpnorm),
2432 &(kin_mem->kin_sFdotJp));
2433
2434 if (retval == 0) {
2435 /* Update gval = uval + gval since gval = -L^{-1}F(uu) */
2436 N_VLinearSum(ONE, uval, ONE, gval, gval);
2437 return(KIN_SUCCESS);
2438 }
2439 else if (retval < 0) return(KIN_LSOLVE_FAIL);
2440 else if ((kin_mem->kin_lsetup == NULL) || (kin_mem->kin_jacCurrent)) return(KIN_LINSOLV_NO_RECOVERY);
2441
2442 /* loop back only if the linear solver setup is in use
2443 and matrix information is not current */
2444
2445 kin_mem->kin_sthrsh = TWO;
2446 }
2447
2448 }
2449
2450
2451 /*
2452 * KINFP
2453 *
2454 * This routine is the main driver for the fixed point iteration with
2455 * Anderson Acceleration.
2456 */
2457
KINFP(KINMem kin_mem)2458 static int KINFP(KINMem kin_mem)
2459 {
2460 int retval; /* return value from user func */
2461 int ret; /* iteration status */
2462 realtype fmax; /* max norm of residual func */
2463 N_Vector delta; /* temporary workspace vector */
2464
2465 delta = kin_mem->kin_vtemp1;
2466 ret = CONTINUE_ITERATIONS;
2467 fmax = kin_mem->kin_fnormtol + ONE;
2468
2469 /* initialize iteration count */
2470 kin_mem->kin_nni = 0;
2471
2472 while (ret == CONTINUE_ITERATIONS) {
2473
2474 /* update iteration count */
2475 kin_mem->kin_nni++;
2476
2477 /* evaluate func(uu) and return if failed */
2478 retval = kin_mem->kin_func(kin_mem->kin_uu, kin_mem->kin_fval,
2479 kin_mem->kin_user_data);
2480 kin_mem->kin_nfe++;
2481
2482 if (retval < 0) {
2483 ret = KIN_SYSFUNC_FAIL;
2484 break;
2485 }
2486
2487 /* compute new solution */
2488 if (kin_mem->kin_m_aa == 0) {
2489 /* standard fixed point */
2490 N_VScale(ONE, kin_mem->kin_fval, kin_mem->kin_unew);
2491 } else {
2492 /* apply Anderson acceleration */
2493 AndersonAcc(kin_mem, kin_mem->kin_fval, delta, kin_mem->kin_unew,
2494 kin_mem->kin_uu, kin_mem->kin_nni - 1, kin_mem->kin_R_aa,
2495 kin_mem->kin_gamma_aa);
2496 }
2497
2498 /* compute change between iterations */
2499 N_VLinearSum(ONE, kin_mem->kin_unew, -ONE, kin_mem->kin_uu, delta);
2500
2501 fmax = KINScFNorm(kin_mem, delta, kin_mem->kin_fscale); /* measure || g(x)-x || */
2502
2503 if (kin_mem->kin_printfl > 1)
2504 KINPrintInfo(kin_mem, PRNT_FMAX, "KINSOL", "KINFP", INFO_FMAX, fmax);
2505
2506 kin_mem->kin_fnorm = fmax;
2507
2508 /* print the current iter, fnorm, and nfe values if printfl > 0 */
2509 if (kin_mem->kin_printfl > 0)
2510 KINPrintInfo(kin_mem, PRNT_NNI, "KINSOL", "KINFP", INFO_NNI,
2511 kin_mem->kin_nni, kin_mem->kin_nfe, kin_mem->kin_fnorm);
2512
2513 /* Check if the maximum number of iterations is reached */
2514 if (kin_mem->kin_nni >= kin_mem->kin_mxiter) {
2515 ret = KIN_MAXITER_REACHED;
2516 }
2517 if (fmax <= kin_mem->kin_fnormtol) {
2518 ret = KIN_SUCCESS;
2519 }
2520
2521 if (ret == CONTINUE_ITERATIONS) {
2522 /* Only update solution if taking a next iteration. */
2523 /* CSW Should put in a conditional to send back the newest iterate or
2524 the one consistent with the fval */
2525 N_VScale(ONE, kin_mem->kin_unew, kin_mem->kin_uu);
2526 }
2527
2528 fflush(kin_mem->kin_errfp);
2529
2530 } /* end of loop; return */
2531
2532 if (kin_mem->kin_printfl > 0)
2533 KINPrintInfo(kin_mem, PRNT_RETVAL, "KINSOL", "KINFP", INFO_RETVAL, ret);
2534
2535 return(ret);
2536 }
2537
2538
2539 /* -----------------------------------------------------------------
2540 * Stopping tests
2541 * -----------------------------------------------------------------
2542 */
2543
2544
2545 /*
2546 * ========================================================================
2547 * Anderson Acceleration
2548 * ========================================================================
2549 */
2550
AndersonAcc(KINMem kin_mem,N_Vector gval,N_Vector fv,N_Vector x,N_Vector xold,long int iter,realtype * R,realtype * gamma)2551 static int AndersonAcc(KINMem kin_mem, N_Vector gval, N_Vector fv,
2552 N_Vector x, N_Vector xold,
2553 long int iter, realtype *R, realtype *gamma)
2554 {
2555 int retval;
2556 long int i_pt, i, j, lAA;
2557 long int *ipt_map;
2558 realtype alfa;
2559 realtype onembeta;
2560 realtype a, b, temp, c, s;
2561
2562 /* local shortcuts for fused vector operation */
2563 int nvec=0;
2564 realtype* cv=kin_mem->kin_cv;
2565 N_Vector* Xv=kin_mem->kin_Xv;
2566
2567 ipt_map = kin_mem->kin_ipt_map;
2568 i_pt = iter-1 - ((iter-1) / kin_mem->kin_m_aa) * kin_mem->kin_m_aa;
2569 N_VLinearSum(ONE, gval, -ONE, xold, fv);
2570 if (iter > 0) {
2571 /* compute dg_new = gval - gval_old */
2572 N_VLinearSum(ONE, gval, -ONE, kin_mem->kin_gold_aa, kin_mem->kin_dg_aa[i_pt]);
2573 /* compute df_new = fval - fval_old */
2574 N_VLinearSum(ONE, fv, -ONE, kin_mem->kin_fold_aa, kin_mem->kin_df_aa[i_pt]);
2575 }
2576
2577 N_VScale(ONE, gval, kin_mem->kin_gold_aa);
2578 N_VScale(ONE, fv, kin_mem->kin_fold_aa);
2579
2580 /* on first iteration, just do basic fixed-point update */
2581 if (iter == 0) {
2582 N_VScale(ONE, gval, x);
2583 return(0);
2584 }
2585
2586 /* update data structures based on current iteration index */
2587
2588 if (iter == 1) {
2589
2590 /* second iteration */
2591 R[0] = SUNRsqrt(N_VDotProd(kin_mem->kin_df_aa[i_pt], kin_mem->kin_df_aa[i_pt]));
2592 alfa = ONE/R[0];
2593 N_VScale(alfa, kin_mem->kin_df_aa[i_pt], kin_mem->kin_q_aa[i_pt]);
2594 ipt_map[0] = 0;
2595
2596 } else if (iter <= kin_mem->kin_m_aa) {
2597
2598 /* another iteration before we've reached maa */
2599 N_VScale(ONE, kin_mem->kin_df_aa[i_pt], kin_mem->kin_vtemp2);
2600 for (j=0; j < (iter-1); j++) {
2601 ipt_map[j] = j;
2602 R[(iter-1)*kin_mem->kin_m_aa+j] = N_VDotProd(kin_mem->kin_q_aa[j], kin_mem->kin_vtemp2);
2603 N_VLinearSum(ONE,kin_mem->kin_vtemp2, -R[(iter-1)*kin_mem->kin_m_aa+j], kin_mem->kin_q_aa[j], kin_mem->kin_vtemp2);
2604 }
2605 R[(iter-1)*kin_mem->kin_m_aa+iter-1] = SUNRsqrt(N_VDotProd(kin_mem->kin_vtemp2, kin_mem->kin_vtemp2));
2606 N_VScale((1/R[(iter-1)*kin_mem->kin_m_aa+iter-1]), kin_mem->kin_vtemp2, kin_mem->kin_q_aa[i_pt]);
2607 ipt_map[iter-1] = iter-1;
2608
2609 } else {
2610
2611 /* we've filled the acceleration subspace, so start recycling */
2612
2613 /* Delete left-most column vector from QR factorization */
2614 for (i=0; i < kin_mem->kin_m_aa-1; i++) {
2615 a = R[(i+1)*kin_mem->kin_m_aa + i];
2616 b = R[(i+1)*kin_mem->kin_m_aa + i+1];
2617 temp = SUNRsqrt(a*a + b*b);
2618 c = a / temp;
2619 s = b / temp;
2620 R[(i+1)*kin_mem->kin_m_aa + i] = temp;
2621 R[(i+1)*kin_mem->kin_m_aa + i+1] = ZERO;
2622 /* OK to re-use temp */
2623 if (i < kin_mem->kin_m_aa-1) {
2624 for (j = i+2; j < kin_mem->kin_m_aa; j++) {
2625 a = R[j*kin_mem->kin_m_aa + i];
2626 b = R[j*kin_mem->kin_m_aa + i+1];
2627 temp = c * a + s * b;
2628 R[j*kin_mem->kin_m_aa + i+1] = -s*a + c*b;
2629 R[j*kin_mem->kin_m_aa + i] = temp;
2630 }
2631 }
2632 N_VLinearSum(c, kin_mem->kin_q_aa[i], s, kin_mem->kin_q_aa[i+1], kin_mem->kin_vtemp2);
2633 N_VLinearSum(-s, kin_mem->kin_q_aa[i], c, kin_mem->kin_q_aa[i+1], kin_mem->kin_q_aa[i+1]);
2634 N_VScale(ONE, kin_mem->kin_vtemp2, kin_mem->kin_q_aa[i]);
2635 }
2636
2637 /* Shift R to the left by one. */
2638 for (i = 1; i < kin_mem->kin_m_aa; i++) {
2639 for (j = 0; j < kin_mem->kin_m_aa-1; j++) {
2640 R[(i-1)*kin_mem->kin_m_aa + j] = R[i*kin_mem->kin_m_aa + j];
2641 }
2642 }
2643
2644 /* Add the new df vector */
2645 N_VScale(ONE, kin_mem->kin_df_aa[i_pt], kin_mem->kin_vtemp2);
2646 for (j=0; j < (kin_mem->kin_m_aa-1); j++) {
2647 R[(kin_mem->kin_m_aa-1)*kin_mem->kin_m_aa+j] = N_VDotProd(kin_mem->kin_q_aa[j], kin_mem->kin_vtemp2);
2648 N_VLinearSum(ONE, kin_mem->kin_vtemp2, -R[(kin_mem->kin_m_aa-1)*kin_mem->kin_m_aa+j], kin_mem->kin_q_aa[j],kin_mem->kin_vtemp2);
2649 }
2650 R[(kin_mem->kin_m_aa-1)*kin_mem->kin_m_aa+kin_mem->kin_m_aa-1] = SUNRsqrt(N_VDotProd(kin_mem->kin_vtemp2, kin_mem->kin_vtemp2));
2651 N_VScale((1/R[(kin_mem->kin_m_aa-1)*kin_mem->kin_m_aa+kin_mem->kin_m_aa-1]), kin_mem->kin_vtemp2, kin_mem->kin_q_aa[kin_mem->kin_m_aa-1]);
2652
2653 /* Update the iteration map */
2654 j = 0;
2655 for (i=i_pt+1; i < kin_mem->kin_m_aa; i++)
2656 ipt_map[j++] = i;
2657 for (i=0; i < (i_pt+1); i++)
2658 ipt_map[j++] = i;
2659 }
2660
2661 /* Solve least squares problem and update solution */
2662 lAA = iter;
2663 if (kin_mem->kin_m_aa < iter) lAA = kin_mem->kin_m_aa;
2664
2665 retval = N_VDotProdMulti((int) lAA, fv, kin_mem->kin_q_aa, gamma);
2666 if (retval != KIN_SUCCESS) return(KIN_VECTOROP_ERR);
2667
2668 /* set arrays for fused vector operation */
2669 cv[0] = ONE;
2670 Xv[0] = gval;
2671 nvec = 1;
2672
2673 for (i=lAA-1; i > -1; i--) {
2674 for (j=i+1; j < lAA; j++) {
2675 gamma[i] = gamma[i]-R[j*kin_mem->kin_m_aa+i]*gamma[j];
2676 }
2677 gamma[i] = gamma[i]/R[i*kin_mem->kin_m_aa+i];
2678
2679 cv[nvec] = -gamma[i];
2680 Xv[nvec] = kin_mem->kin_dg_aa[ipt_map[i]];
2681 nvec += 1;
2682 }
2683
2684 /* if enabled, apply damping */
2685 if (kin_mem->kin_damping_aa) {
2686 onembeta = (ONE - kin_mem->kin_beta_aa);
2687 cv[nvec] = -onembeta;
2688 Xv[nvec] = fv;
2689 nvec += 1;
2690 for (i = lAA - 1; i > -1; i--) {
2691 cv[nvec] = onembeta * gamma[i];
2692 Xv[nvec] = kin_mem->kin_df_aa[ipt_map[i]];
2693 nvec += 1;
2694 }
2695 }
2696
2697 /* update solution */
2698 retval = N_VLinearCombination(nvec, cv, Xv, x);
2699 if (retval != KIN_SUCCESS) return(KIN_VECTOROP_ERR);
2700
2701 return 0;
2702 }
2703