1 /*
2 * -----------------------------------------------------------------
3 * $Revision: 4272 $
4 * $Date: 2014-12-02 11:19:41 -0800 (Tue, 02 Dec 2014) $
5 * -----------------------------------------------------------------
6 * Programmer(s): Radu Serban and Aaron Collier @ LLNL
7 * -----------------------------------------------------------------
8 * LLNS Copyright Start
9 * Copyright (c) 2014, Lawrence Livermore National Security
10 * This work was performed under the auspices of the U.S. Department
11 * of Energy by Lawrence Livermore National Laboratory in part under
12 * Contract W-7405-Eng-48 and in part under Contract DE-AC52-07NA27344.
13 * Produced at the Lawrence Livermore National Laboratory.
14 * All rights reserved.
15 * For details, see the LICENSE file.
16 * LLNS Copyright End
17 * -----------------------------------------------------------------
18 * This is the implementation file for the KINSPILS linear solvers.
19 * -----------------------------------------------------------------
20 */
21
22 #include <stdio.h>
23 #include <stdlib.h>
24 #include <stdarg.h>
25
26 #include "kinsol_impl.h"
27 #include "kinsol_spils_impl.h"
28
29 #include <sundials/sundials_math.h>
30
31 /*
32 * -----------------------------------------------------------------
33 * private constants
34 * -----------------------------------------------------------------
35 */
36
37 #define ZERO RCONST(0.0)
38 #define ONE RCONST(1.0)
39 #define TWO RCONST(2.0)
40
41
42 /*
43 * -----------------------------------------------------------------
44 * readability replacements
45 * -----------------------------------------------------------------
46 */
47
48 #define lrw1 (kin_mem->kin_lrw1)
49 #define liw1 (kin_mem->kin_liw1)
50 #define func (kin_mem->kin_func)
51 #define user_data (kin_mem->kin_user_data)
52 #define printfl (kin_mem->kin_printfl)
53 #define lmem (kin_mem->kin_lmem)
54 #define uu (kin_mem->kin_uu)
55 #define fval (kin_mem->kin_fval)
56 #define uscale (kin_mem->kin_uscale)
57 #define fscale (kin_mem->kin_fscale)
58 #define sqrt_relfunc (kin_mem->kin_sqrt_relfunc)
59 #define eps (kin_mem->kin_eps)
60 #define errfp (kin_mem->kin_errfp)
61 #define infofp (kin_mem->kin_infofp)
62 #define vtemp1 (kin_mem->kin_vtemp1)
63 #define vec_tmpl (kin_mem->kin_vtemp1)
64 #define vtemp2 (kin_mem->kin_vtemp2)
65
66 #define ils_type (kinspils_mem->s_type)
67 #define pretype (kinspils_mem->s_pretype)
68 #define gstype (kinspils_mem->s_gstype)
69 #define nli (kinspils_mem->s_nli)
70 #define npe (kinspils_mem->s_npe)
71 #define nps (kinspils_mem->s_nps)
72 #define ncfl (kinspils_mem->s_ncfl)
73 #define njtimes (kinspils_mem->s_njtimes)
74 #define nfes (kinspils_mem->s_nfes)
75 #define new_uu (kinspils_mem->s_new_uu)
76
77 #define jtimesDQ (kinspils_mem->s_jtimesDQ)
78 #define jtimes (kinspils_mem->s_jtimes)
79 #define J_data (kinspils_mem->s_J_data)
80
81 #define last_flag (kinspils_mem->s_last_flag)
82
83
84 /*
85 * -----------------------------------------------------------------
86 * Function : KINSpilsSetMaxRestarts
87 * -----------------------------------------------------------------
88 */
89
KINSpilsSetMaxRestarts(void * kinmem,int maxrs)90 int KINSpilsSetMaxRestarts(void *kinmem, int maxrs)
91 {
92 KINMem kin_mem;
93 KINSpilsMem kinspils_mem;
94
95 /* return immediately if kinmem is NULL */
96
97 if (kinmem == NULL) {
98 KINProcessError(NULL, KINSPILS_MEM_NULL, "KINSPILS", "KINSpilsSetMaxRestarts", MSGS_KINMEM_NULL);
99 return(KINSPILS_MEM_NULL);
100 }
101 kin_mem = (KINMem) kinmem;
102
103 if (lmem == NULL) {
104 KINProcessError(kin_mem, KINSPILS_LMEM_NULL, "KINSPILS", "KINSpilsSetMaxRestarts", MSGS_LMEM_NULL);
105 return(KINSPILS_LMEM_NULL);
106 }
107 kinspils_mem = (KINSpilsMem) lmem;
108
109 /* check for legal maxrs */
110
111 if (maxrs < 0) {
112 KINProcessError(kin_mem, KINSPILS_ILL_INPUT, "KINSPILS", "KINSpilsSetMaxRestarts", MSGS_NEG_MAXRS);
113 return(KINSPILS_ILL_INPUT);
114 }
115 kinspils_mem->s_maxlrst = maxrs;
116
117 return(KINSPILS_SUCCESS);
118 }
119
120 /*
121 * -----------------------------------------------------------------
122 * Function : KINSpilsSetPreconditioner
123 * -----------------------------------------------------------------
124 */
125
KINSpilsSetPreconditioner(void * kinmem,KINSpilsPrecSetupFn pset,KINSpilsPrecSolveFn psolve)126 int KINSpilsSetPreconditioner(void *kinmem,
127 KINSpilsPrecSetupFn pset, KINSpilsPrecSolveFn psolve)
128 {
129 KINMem kin_mem;
130 KINSpilsMem kinspils_mem;
131
132 /* return immediately if kinmem is NULL */
133
134 if (kinmem == NULL) {
135 KINProcessError(NULL, KINSPILS_MEM_NULL, "KINSPILS", "KINSpilsSetPreconditioner", MSGS_KINMEM_NULL);
136 return(KINSPILS_MEM_NULL);
137 }
138 kin_mem = (KINMem) kinmem;
139
140 if (lmem == NULL) {
141 KINProcessError(kin_mem, KINSPILS_LMEM_NULL, "KINSPILS", "KINSpilsSetPreconditioner", MSGS_LMEM_NULL);
142 return(KINSPILS_LMEM_NULL);
143 }
144 kinspils_mem = (KINSpilsMem) lmem;
145
146 kinspils_mem->s_pset = pset;
147 kinspils_mem->s_psolve = psolve;
148
149 return(KINSPILS_SUCCESS);
150 }
151
152 /*
153 * -----------------------------------------------------------------
154 * Function : KINSpilsSetJacTimesVecFn
155 * -----------------------------------------------------------------
156 */
157
KINSpilsSetJacTimesVecFn(void * kinmem,KINSpilsJacTimesVecFn jtv)158 int KINSpilsSetJacTimesVecFn(void *kinmem, KINSpilsJacTimesVecFn jtv)
159
160 {
161 KINMem kin_mem;
162 KINSpilsMem kinspils_mem;
163
164 /* return immediately if kinmem is NULL */
165
166 if (kinmem == NULL) {
167 KINProcessError(NULL, KINSPILS_MEM_NULL, "KINSPILS", "KINSpilsSetJacTimesVecFn", MSGS_KINMEM_NULL);
168 return(KINSPILS_MEM_NULL);
169 }
170 kin_mem = (KINMem) kinmem;
171
172 if (lmem == NULL) {
173 KINProcessError(kin_mem, KINSPILS_LMEM_NULL, "KINSPILS", "KINSpilsSetJacTimesVecFn", MSGS_LMEM_NULL);
174 return(KINSPILS_LMEM_NULL);
175 }
176 kinspils_mem = (KINSpilsMem) lmem;
177
178 if (jtv != NULL) {
179 jtimesDQ = FALSE;
180 jtimes = jtv;
181 } else {
182 jtimesDQ = TRUE;
183 }
184
185 return(KINSPILS_SUCCESS);
186 }
187
188 /*
189 * -----------------------------------------------------------------
190 * Function : KINSpilsGetWorkSpace
191 * -----------------------------------------------------------------
192 */
193
KINSpilsGetWorkSpace(void * kinmem,long int * lenrwSG,long int * leniwSG)194 int KINSpilsGetWorkSpace(void *kinmem, long int *lenrwSG, long int *leniwSG)
195 {
196 KINMem kin_mem;
197 KINSpilsMem kinspils_mem;
198 int maxl;
199
200 /* return immediately if kinmem is NULL */
201
202 if (kinmem == NULL) {
203 KINProcessError(NULL, KINSPILS_MEM_NULL, "KINSPILS", "KINSpilsGetWorkSpace", MSGS_KINMEM_NULL);
204 return(KINSPILS_MEM_NULL);
205 }
206 kin_mem = (KINMem) kinmem;
207
208 if (lmem == NULL) {
209 KINProcessError(kin_mem, KINSPILS_LMEM_NULL, "KINSPILS", "KINSpilsGetWorkSpace", MSGS_LMEM_NULL);
210 return(KINSPILS_LMEM_NULL);
211 }
212 kinspils_mem = (KINSpilsMem) lmem;
213
214 maxl = kinspils_mem->s_maxl;
215
216 switch(ils_type) {
217 case SPILS_SPGMR:
218 *lenrwSG = lrw1 * (maxl + 3) + (maxl * (maxl + 4)) + 1;
219 *leniwSG = liw1 * (maxl + 3);
220 break;
221 case SPILS_SPBCG:
222 *lenrwSG = lrw1 * 7;
223 *leniwSG = liw1 * 7;
224 break;
225 case SPILS_SPTFQMR:
226 *lenrwSG = lrw1 * 11;
227 *leniwSG = liw1 * 11;
228 break;
229 }
230
231
232 return(KINSPILS_SUCCESS);
233 }
234
235 /*
236 * -----------------------------------------------------------------
237 * Function : KINSpilsGetNumPrecEvals
238 * -----------------------------------------------------------------
239 */
240
KINSpilsGetNumPrecEvals(void * kinmem,long int * npevals)241 int KINSpilsGetNumPrecEvals(void *kinmem, long int *npevals)
242 {
243 KINMem kin_mem;
244 KINSpilsMem kinspils_mem;
245
246 /* return immediately if kinmem is NULL */
247
248 if (kinmem == NULL) {
249 KINProcessError(NULL, KINSPILS_MEM_NULL, "KINSPILS", "KINSpilsGetNumPrecEvals", MSGS_KINMEM_NULL);
250 return(KINSPILS_MEM_NULL);
251 }
252 kin_mem = (KINMem) kinmem;
253
254 if (lmem == NULL) {
255 KINProcessError(kin_mem, KINSPILS_LMEM_NULL, "KINSPILS", "KINSpilsGetNumPrecEvals", MSGS_LMEM_NULL);
256 return(KINSPILS_LMEM_NULL);
257 }
258 kinspils_mem = (KINSpilsMem) lmem;
259 *npevals = npe;
260
261 return(KINSPILS_SUCCESS);
262 }
263
264 /*
265 * -----------------------------------------------------------------
266 * Function : KINSpilsGetNumPrecSolves
267 * -----------------------------------------------------------------
268 */
269
KINSpilsGetNumPrecSolves(void * kinmem,long int * npsolves)270 int KINSpilsGetNumPrecSolves(void *kinmem, long int *npsolves)
271 {
272 KINMem kin_mem;
273 KINSpilsMem kinspils_mem;
274
275 /* return immediately if kinmem is NULL */
276
277 if (kinmem == NULL) {
278 KINProcessError(NULL, KINSPILS_MEM_NULL, "KINSPILS", "KINSpilsGetNumPrecSolves", MSGS_KINMEM_NULL);
279 return(KINSPILS_MEM_NULL);
280 }
281 kin_mem = (KINMem) kinmem;
282
283 if (lmem == NULL) {
284 KINProcessError(kin_mem, KINSPILS_LMEM_NULL, "KINSPILS", "KINSpilsGetNumPrecSolves", MSGS_LMEM_NULL);
285 return(KINSPILS_LMEM_NULL);
286 }
287 kinspils_mem = (KINSpilsMem) lmem;
288 *npsolves = nps;
289
290 return(KINSPILS_SUCCESS);
291 }
292
293 /*
294 * -----------------------------------------------------------------
295 * Function : KINSpilsGetNumLinIters
296 * -----------------------------------------------------------------
297 */
298
KINSpilsGetNumLinIters(void * kinmem,long int * nliters)299 int KINSpilsGetNumLinIters(void *kinmem, long int *nliters)
300 {
301 KINMem kin_mem;
302 KINSpilsMem kinspils_mem;
303
304 /* return immediately if kinmem is NULL */
305
306 if (kinmem == NULL) {
307 KINProcessError(NULL, KINSPILS_MEM_NULL, "KINSPILS", "KINSpilsGetNumLinIters", MSGS_KINMEM_NULL);
308 return(KINSPILS_MEM_NULL);
309 }
310 kin_mem = (KINMem) kinmem;
311
312 if (lmem == NULL) {
313 KINProcessError(kin_mem, KINSPILS_LMEM_NULL, "KINSPILS", "KINSpilsGetNumLinIters", MSGS_LMEM_NULL);
314 return(KINSPILS_LMEM_NULL);
315 }
316 kinspils_mem = (KINSpilsMem) lmem;
317 *nliters = nli;
318
319 return(KINSPILS_SUCCESS);
320 }
321
322 /*
323 * -----------------------------------------------------------------
324 * Function : KINSpilsGetNumConvFails
325 * -----------------------------------------------------------------
326 */
327
KINSpilsGetNumConvFails(void * kinmem,long int * nlcfails)328 int KINSpilsGetNumConvFails(void *kinmem, long int *nlcfails)
329 {
330 KINMem kin_mem;
331 KINSpilsMem kinspils_mem;
332
333 /* return immediately if kinmem is NULL */
334
335 if (kinmem == NULL) {
336 KINProcessError(NULL, KINSPILS_MEM_NULL, "KINSPILS", "KINSpilsGetNumConvFails", MSGS_KINMEM_NULL);
337 return(KINSPILS_MEM_NULL);
338 }
339 kin_mem = (KINMem) kinmem;
340
341 if (lmem == NULL) {
342 KINProcessError(kin_mem, KINSPILS_LMEM_NULL, "KINSPILS", "KINSpilsGetNumConvFails", MSGS_LMEM_NULL);
343 return(KINSPILS_LMEM_NULL);
344 }
345 kinspils_mem = (KINSpilsMem) lmem;
346 *nlcfails = ncfl;
347
348 return(KINSPILS_SUCCESS);
349 }
350
351 /*
352 * -----------------------------------------------------------------
353 * Function : KINSpilsGetNumJtimesEvals
354 * -----------------------------------------------------------------
355 */
356
KINSpilsGetNumJtimesEvals(void * kinmem,long int * njvevals)357 int KINSpilsGetNumJtimesEvals(void *kinmem, long int *njvevals)
358 {
359 KINMem kin_mem;
360 KINSpilsMem kinspils_mem;
361
362 /* return immediately if kinmem is NULL */
363
364 if (kinmem == NULL) {
365 KINProcessError(NULL, KINSPILS_MEM_NULL, "KINSPILS", "KINSpilsGetNumJtimesEvals", MSGS_KINMEM_NULL);
366 return(KINSPILS_MEM_NULL);
367 }
368 kin_mem = (KINMem) kinmem;
369
370 if (lmem == NULL) {
371 KINProcessError(kin_mem, KINSPILS_LMEM_NULL, "KINSPILS", "KINSpilsGetNumJtimesEvals", MSGS_LMEM_NULL);
372 return(KINSPILS_LMEM_NULL);
373 }
374 kinspils_mem = (KINSpilsMem) lmem;
375 *njvevals = njtimes;
376
377 return(KINSPILS_SUCCESS);
378 }
379
380 /*
381 * -----------------------------------------------------------------
382 * Function : KINSpilsGetNumFuncEvals
383 * -----------------------------------------------------------------
384 */
385
KINSpilsGetNumFuncEvals(void * kinmem,long int * nfevalsS)386 int KINSpilsGetNumFuncEvals(void *kinmem, long int *nfevalsS)
387 {
388 KINMem kin_mem;
389 KINSpilsMem kinspils_mem;
390
391 /* return immediately if kinmem is NULL */
392
393 if (kinmem == NULL) {
394 KINProcessError(NULL, KINSPILS_MEM_NULL, "KINSPILS", "KINSpilsGetNumFuncEvals", MSGS_KINMEM_NULL);
395 return(KINSPILS_MEM_NULL);
396 }
397 kin_mem = (KINMem) kinmem;
398
399 if (lmem == NULL) {
400 KINProcessError(kin_mem, KINSPILS_LMEM_NULL, "KINSPILS", "KINSpilsGetNumFuncEvals", MSGS_LMEM_NULL);
401 return(KINSPILS_LMEM_NULL);
402 }
403 kinspils_mem = (KINSpilsMem) lmem;
404 *nfevalsS = nfes;
405
406 return(KINSPILS_SUCCESS);
407 }
408
409 /*
410 * -----------------------------------------------------------------
411 * Function : KINSpilsGetLastFlag
412 * -----------------------------------------------------------------
413 */
414
KINSpilsGetLastFlag(void * kinmem,long int * flag)415 int KINSpilsGetLastFlag(void *kinmem, long int *flag)
416 {
417 KINMem kin_mem;
418 KINSpilsMem kinspils_mem;
419
420 /* return immediately if kinmem is NULL */
421
422 if (kinmem == NULL) {
423 KINProcessError(NULL, KINSPILS_MEM_NULL, "KINSPILS", "KINSpilsGetLastFlag", MSGS_KINMEM_NULL);
424 return(KINSPILS_MEM_NULL);
425 }
426 kin_mem = (KINMem) kinmem;
427
428 if (lmem == NULL) {
429 KINProcessError(kin_mem, KINSPILS_LMEM_NULL, "KINSPILS", "KINSpilsGetLastFlag", MSGS_LMEM_NULL);
430 return(KINSPILS_LMEM_NULL);
431 }
432 kinspils_mem = (KINSpilsMem) lmem;
433
434 *flag = last_flag;
435
436 return(KINSPILS_SUCCESS);
437 }
438
439 /*
440 * -----------------------------------------------------------------
441 * Function : KINSpilsGetReturnFlagName
442 * -----------------------------------------------------------------
443 */
444
KINSpilsGetReturnFlagName(long int flag)445 char *KINSpilsGetReturnFlagName(long int flag)
446 {
447 char *name;
448
449 name = (char *)malloc(30*sizeof(char));
450
451 switch(flag) {
452 case KINSPILS_SUCCESS:
453 sprintf(name, "KINSPILS_SUCCESS");
454 break;
455 case KINSPILS_MEM_NULL:
456 sprintf(name, "KINSPILS_MEM_NULL");
457 break;
458 case KINSPILS_LMEM_NULL:
459 sprintf(name, "KINSPILS_LMEM_NULL");
460 break;
461 case KINSPILS_ILL_INPUT:
462 sprintf(name, "KINSPILS_ILL_INPUT");
463 break;
464 case KINSPILS_MEM_FAIL:
465 sprintf(name, "KINSPILS_MEM_FAIL");
466 break;
467 case KINSPILS_PMEM_NULL:
468 sprintf(name, "KINSPILS_PMEM_NULL");
469 break;
470 default:
471 sprintf(name, "NONE");
472 }
473
474 return(name);
475 }
476
477 /*
478 * -----------------------------------------------------------------
479 * additional readability replacements
480 * -----------------------------------------------------------------
481 */
482
483 #define maxl (kinspils_mem->s_maxl)
484 #define maxlrst (kinspils_mem->s_maxlrst)
485 #define pset (kinspils_mem->s_pset)
486 #define psolve (kinspils_mem->s_psolve)
487 #define P_data (kinspils_mem->s_P_data)
488
489 /*
490 * -----------------------------------------------------------------
491 * Function : KINSpilsAtimes
492 * -----------------------------------------------------------------
493 * This routine coordinates the generation of the matrix-vector
494 * product z = J*v by calling either KINSpilsDQJtimes, which uses
495 * a difference quotient approximation for J*v, or by calling the
496 * user-supplied routine KINSpilsJacTimesVecFn if it is non-null.
497 * -----------------------------------------------------------------
498 */
499
KINSpilsAtimes(void * kinsol_mem,N_Vector v,N_Vector z)500 int KINSpilsAtimes(void *kinsol_mem, N_Vector v, N_Vector z)
501 {
502 KINMem kin_mem;
503 KINSpilsMem kinspils_mem;
504 int ret;
505
506 kin_mem = (KINMem) kinsol_mem;
507 kinspils_mem = (KINSpilsMem) lmem;
508
509 ret = jtimes(v, z, uu, &new_uu, J_data);
510 njtimes++;
511
512 return(ret);
513 }
514
515 /*
516 * -----------------------------------------------------------------
517 * Function : KINSpilsPSolve
518 * -----------------------------------------------------------------
519 * This routine interfaces between the generic Sp***Solve routine
520 * (within the SPGMR, SPBCG, or SPTFQMR solver) and the
521 * user's psolve routine. It passes to psolve all required state
522 * information from kinsol_mem. Its return value is the same as that
523 * returned by psolve. Note that the generic SP*** solver guarantees
524 * that KINSpilsPSolve will not be called in the case in which
525 * preconditioning is not done. This is the only case in which the
526 * user's psolve routine is allowed to be NULL.
527 * -----------------------------------------------------------------
528 */
529
KINSpilsPSolve(void * kinsol_mem,N_Vector r,N_Vector z,int lrdummy)530 int KINSpilsPSolve(void *kinsol_mem, N_Vector r, N_Vector z, int lrdummy)
531 {
532 KINMem kin_mem;
533 KINSpilsMem kinspils_mem;
534 int ret;
535
536 kin_mem = (KINMem) kinsol_mem;
537 kinspils_mem = (KINSpilsMem) lmem;
538
539 /* copy the rhs into z before the psolve call */
540 /* Note: z returns with the solution */
541
542 N_VScale(ONE, r, z);
543
544 /* this call is counted in nps within the KINSpilsSolve routine */
545
546 ret = psolve(uu, uscale, fval, fscale, z, P_data, vtemp1);
547
548 return(ret);
549 }
550
551 /*
552 * -----------------------------------------------------------------
553 * Function : KINSpilsDQJtimes
554 * -----------------------------------------------------------------
555 * This routine generates the matrix-vector product z = J*v using a
556 * difference quotient approximation. The approximation is
557 * J*v = [func(uu + sigma*v) - func(uu)]/sigma. Here sigma is based
558 * on the dot products (uscale*uu, uscale*v) and
559 * (uscale*v, uscale*v), the L1Norm(uscale*v), and on sqrt_relfunc
560 * (the square root of the relative error in the function). Note
561 * that v in the argument list has already been both preconditioned
562 * and unscaled.
563 *
564 * NOTE: Unlike the DQ Jacobian functions for direct linear solvers
565 * (which are called from within the lsetup function), this
566 * function is called from within the lsolve function and thus
567 * a recovery may still be possible even if the system function
568 * fails (recoverably).
569 * -----------------------------------------------------------------
570 */
571
KINSpilsDQJtimes(N_Vector v,N_Vector Jv,N_Vector u,booleantype * new_u,void * data)572 int KINSpilsDQJtimes(N_Vector v, N_Vector Jv,
573 N_Vector u, booleantype *new_u,
574 void *data)
575 {
576 realtype sigma, sigma_inv, sutsv, sq1norm, sign, vtv;
577 KINMem kin_mem;
578 KINSpilsMem kinspils_mem;
579 int retval;
580
581 /* data is kin_mem */
582
583 kin_mem = (KINMem) data;
584 kinspils_mem = (KINSpilsMem) lmem;
585
586 /* scale the vector v and put Du*v into vtemp1 */
587
588 N_VProd(v, uscale, vtemp1);
589
590 /* scale u and put into Jv (used as a temporary storage) */
591
592 N_VProd(u, uscale, Jv);
593
594 /* compute dot product (Du*u).(Du*v) */
595
596 sutsv = N_VDotProd(Jv, vtemp1);
597
598 /* compute dot product (Du*v).(Du*v) */
599
600 vtv = N_VDotProd(vtemp1, vtemp1);
601
602 sq1norm = N_VL1Norm(vtemp1);
603
604 sign = (sutsv >= ZERO) ? ONE : -ONE ;
605
606 /* this expression for sigma is from p. 469, Brown and Saad paper */
607
608 sigma = sign*sqrt_relfunc*SUNMAX(SUNRabs(sutsv),sq1norm)/vtv;
609
610 sigma_inv = ONE/sigma;
611
612 /* compute the u-prime at which to evaluate the function func */
613
614 N_VLinearSum(ONE, u, sigma, v, vtemp1);
615
616 /* call the system function to calculate func(u+sigma*v) */
617
618 retval = func(vtemp1, vtemp2, user_data);
619 nfes++;
620 if (retval != 0) return(retval);
621
622 /* finish the computation of the difference quotient */
623
624 N_VLinearSum(sigma_inv, vtemp2, -sigma_inv, fval, Jv);
625
626 return(0);
627 }
628
629