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