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