xref: /petsc/src/snes/impls/vi/ss/viss.c (revision cf1aed2ce99d23e50336629af3ca8cf096900abb)
1 
2 #include <../src/snes/impls/vi/ss/vissimpl.h> /*I "petscsnes.h" I*/
3 #include <../include/petsc-private/kspimpl.h>
4 #include <../include/petsc-private/matimpl.h>
5 #include <../include/petsc-private/dmimpl.h>
6 
7 
8 /*
9   SNESVIComputeMeritFunction - Evaluates the merit function for the mixed complementarity problem.
10 
11   Input Parameter:
12 . phi - the semismooth function
13 
14   Output Parameter:
15 . merit - the merit function
16 . phinorm - ||phi||
17 
18   Notes:
19   The merit function for the mixed complementarity problem is defined as
20      merit = 0.5*phi^T*phi
21 */
22 #undef __FUNCT__
23 #define __FUNCT__ "SNESVIComputeMeritFunction"
24 static PetscErrorCode SNESVIComputeMeritFunction(Vec phi, PetscReal *merit,PetscReal *phinorm)
25 {
26   PetscErrorCode ierr;
27 
28   PetscFunctionBegin;
29   ierr = VecNormBegin(phi,NORM_2,phinorm);CHKERRQ(ierr);
30   ierr = VecNormEnd(phi,NORM_2,phinorm);CHKERRQ(ierr);
31 
32   *merit = 0.5*(*phinorm)*(*phinorm);
33   PetscFunctionReturn(0);
34 }
35 
36 PETSC_STATIC_INLINE PetscScalar Phi(PetscScalar a,PetscScalar b)
37 {
38   return a + b - PetscSqrtScalar(a*a + b*b);
39 }
40 
41 PETSC_STATIC_INLINE PetscScalar DPhi(PetscScalar a,PetscScalar b)
42 {
43   if ((PetscAbsScalar(a) >= 1.e-6) || (PetscAbsScalar(b) >= 1.e-6)) return 1.0 - a/ PetscSqrtScalar(a*a + b*b);
44   else return .5;
45 }
46 
47 /*
48    SNESVIComputeFunction - Reformulates a system of nonlinear equations in mixed complementarity form to a system of nonlinear equations in semismooth form.
49 
50    Input Parameters:
51 .  snes - the SNES context
52 .  x - current iterate
53 .  functx - user defined function context
54 
55    Output Parameters:
56 .  phi - Semismooth function
57 
58 */
59 #undef __FUNCT__
60 #define __FUNCT__ "SNESVIComputeFunction"
61 static PetscErrorCode SNESVIComputeFunction(SNES snes,Vec X,Vec phi,void *functx)
62 {
63   PetscErrorCode    ierr;
64   SNES_VINEWTONSSLS *vi = (SNES_VINEWTONSSLS*)snes->data;
65   Vec               Xl  = snes->xl,Xu = snes->xu,F = snes->vec_func;
66   PetscScalar       *phi_arr,*x_arr,*f_arr,*l,*u;
67   PetscInt          i,nlocal;
68 
69   PetscFunctionBegin;
70   ierr = (*vi->computeuserfunction)(snes,X,F,functx);CHKERRQ(ierr);
71   ierr = VecGetLocalSize(X,&nlocal);CHKERRQ(ierr);
72   ierr = VecGetArray(X,&x_arr);CHKERRQ(ierr);
73   ierr = VecGetArray(F,&f_arr);CHKERRQ(ierr);
74   ierr = VecGetArray(Xl,&l);CHKERRQ(ierr);
75   ierr = VecGetArray(Xu,&u);CHKERRQ(ierr);
76   ierr = VecGetArray(phi,&phi_arr);CHKERRQ(ierr);
77 
78   for (i=0; i < nlocal; i++) {
79     if ((PetscRealPart(l[i]) <= PETSC_NINFINITY) && (PetscRealPart(u[i]) >= PETSC_INFINITY)) { /* no constraints on variable */
80       phi_arr[i] = f_arr[i];
81     } else if (PetscRealPart(l[i]) <= PETSC_NINFINITY) {                      /* upper bound on variable only */
82       phi_arr[i] = -Phi(u[i] - x_arr[i],-f_arr[i]);
83     } else if (PetscRealPart(u[i]) >= PETSC_INFINITY) {                       /* lower bound on variable only */
84       phi_arr[i] = Phi(x_arr[i] - l[i],f_arr[i]);
85     } else if (l[i] == u[i]) {
86       phi_arr[i] = l[i] - x_arr[i];
87     } else {                                                /* both bounds on variable */
88       phi_arr[i] = Phi(x_arr[i] - l[i],-Phi(u[i] - x_arr[i],-f_arr[i]));
89     }
90   }
91 
92   ierr = VecRestoreArray(X,&x_arr);CHKERRQ(ierr);
93   ierr = VecRestoreArray(F,&f_arr);CHKERRQ(ierr);
94   ierr = VecRestoreArray(Xl,&l);CHKERRQ(ierr);
95   ierr = VecRestoreArray(Xu,&u);CHKERRQ(ierr);
96   ierr = VecRestoreArray(phi,&phi_arr);CHKERRQ(ierr);
97   PetscFunctionReturn(0);
98 }
99 
100 /*
101    SNESVIComputeBsubdifferentialVectors - Computes the diagonal shift (Da) and row scaling (Db) vectors needed for the
102                                           the semismooth jacobian.
103 */
104 #undef __FUNCT__
105 #define __FUNCT__ "SNESVIComputeBsubdifferentialVectors"
106 PetscErrorCode SNESVIComputeBsubdifferentialVectors(SNES snes,Vec X,Vec F,Mat jac,Vec Da,Vec Db)
107 {
108   PetscErrorCode ierr;
109   PetscScalar    *l,*u,*x,*f,*da,*db,da1,da2,db1,db2;
110   PetscInt       i,nlocal;
111 
112   PetscFunctionBegin;
113   ierr = VecGetArray(X,&x);CHKERRQ(ierr);
114   ierr = VecGetArray(F,&f);CHKERRQ(ierr);
115   ierr = VecGetArray(snes->xl,&l);CHKERRQ(ierr);
116   ierr = VecGetArray(snes->xu,&u);CHKERRQ(ierr);
117   ierr = VecGetArray(Da,&da);CHKERRQ(ierr);
118   ierr = VecGetArray(Db,&db);CHKERRQ(ierr);
119   ierr = VecGetLocalSize(X,&nlocal);CHKERRQ(ierr);
120 
121   for (i=0; i< nlocal; i++) {
122     if ((PetscRealPart(l[i]) <= PETSC_NINFINITY) && (PetscRealPart(u[i]) >= PETSC_INFINITY)) { /* no constraints on variable */
123       da[i] = 0;
124       db[i] = 1;
125     } else if (PetscRealPart(l[i]) <= PETSC_NINFINITY) {                     /* upper bound on variable only */
126       da[i] = DPhi(u[i] - x[i], -f[i]);
127       db[i] = DPhi(-f[i],u[i] - x[i]);
128     } else if (PetscRealPart(u[i]) >= PETSC_INFINITY) {                      /* lower bound on variable only */
129       da[i] = DPhi(x[i] - l[i], f[i]);
130       db[i] = DPhi(f[i],x[i] - l[i]);
131     } else if (l[i] == u[i]) {                              /* fixed variable */
132       da[i] = 1;
133       db[i] = 0;
134     } else {                                                /* upper and lower bounds on variable */
135       da1   = DPhi(x[i] - l[i], -Phi(u[i] - x[i], -f[i]));
136       db1   = DPhi(-Phi(u[i] - x[i], -f[i]),x[i] - l[i]);
137       da2   = DPhi(u[i] - x[i], -f[i]);
138       db2   = DPhi(-f[i],u[i] - x[i]);
139       da[i] = da1 + db1*da2;
140       db[i] = db1*db2;
141     }
142   }
143 
144   ierr = VecRestoreArray(X,&x);CHKERRQ(ierr);
145   ierr = VecRestoreArray(F,&f);CHKERRQ(ierr);
146   ierr = VecRestoreArray(snes->xl,&l);CHKERRQ(ierr);
147   ierr = VecRestoreArray(snes->xu,&u);CHKERRQ(ierr);
148   ierr = VecRestoreArray(Da,&da);CHKERRQ(ierr);
149   ierr = VecRestoreArray(Db,&db);CHKERRQ(ierr);
150   PetscFunctionReturn(0);
151 }
152 
153 /*
154    SNESVIComputeJacobian - Computes the jacobian of the semismooth function.The Jacobian for the semismooth function is an element of the B-subdifferential of the Fischer-Burmeister function for complementarity problems.
155 
156    Input Parameters:
157 .  Da       - Diagonal shift vector for the semismooth jacobian.
158 .  Db       - Row scaling vector for the semismooth jacobian.
159 
160    Output Parameters:
161 .  jac      - semismooth jacobian
162 .  jac_pre  - optional preconditioning matrix
163 
164    Notes:
165    The semismooth jacobian matrix is given by
166    jac = Da + Db*jacfun
167    where Db is the row scaling matrix stored as a vector,
168          Da is the diagonal perturbation matrix stored as a vector
169    and   jacfun is the jacobian of the original nonlinear function.
170 */
171 #undef __FUNCT__
172 #define __FUNCT__ "SNESVIComputeJacobian"
173 PetscErrorCode SNESVIComputeJacobian(Mat jac, Mat jac_pre,Vec Da, Vec Db)
174 {
175   PetscErrorCode ierr;
176 
177   /* Do row scaling  and add diagonal perturbation */
178   ierr = MatDiagonalScale(jac,Db,NULL);CHKERRQ(ierr);
179   ierr = MatDiagonalSet(jac,Da,ADD_VALUES);CHKERRQ(ierr);
180   if (jac != jac_pre) { /* If jac and jac_pre are different */
181     ierr = MatDiagonalScale(jac_pre,Db,NULL);CHKERRQ(ierr);
182     ierr = MatDiagonalSet(jac_pre,Da,ADD_VALUES);CHKERRQ(ierr);
183   }
184   PetscFunctionReturn(0);
185 }
186 
187 /*
188    SNESVIComputeMeritFunctionGradient - Computes the gradient of the merit function psi.
189 
190    Input Parameters:
191    phi - semismooth function.
192    H   - semismooth jacobian
193 
194    Output Parameters:
195    dpsi - merit function gradient
196 
197    Notes:
198   The merit function gradient is computed as follows
199         dpsi = H^T*phi
200 */
201 #undef __FUNCT__
202 #define __FUNCT__ "SNESVIComputeMeritFunctionGradient"
203 PetscErrorCode SNESVIComputeMeritFunctionGradient(Mat H, Vec phi, Vec dpsi)
204 {
205   PetscErrorCode ierr;
206 
207   PetscFunctionBegin;
208   ierr = MatMultTranspose(H,phi,dpsi);CHKERRQ(ierr);
209   PetscFunctionReturn(0);
210 }
211 
212 
213 
214 /*
215    SNESSolve_VINEWTONSSLS - Solves the complementarity problem with a semismooth Newton
216    method using a line search.
217 
218    Input Parameters:
219 .  snes - the SNES context
220 
221    Output Parameter:
222 .  outits - number of iterations until termination
223 
224    Application Interface Routine: SNESSolve()
225 
226    Notes:
227    This implements essentially a semismooth Newton method with a
228    line search. The default line search does not do any line seach
229    but rather takes a full newton step.
230 
231    Developer Note: the code in this file should be slightly modified so that this routine need not exist and the SNESSolve_NEWTONLS() routine is called directly with the appropriate wrapped function and Jacobian evaluations
232 
233 */
234 #undef __FUNCT__
235 #define __FUNCT__ "SNESSolve_VINEWTONSSLS"
236 PetscErrorCode SNESSolve_VINEWTONSSLS(SNES snes)
237 {
238   SNES_VINEWTONSSLS  *vi = (SNES_VINEWTONSSLS*)snes->data;
239   PetscErrorCode     ierr;
240   PetscInt           maxits,i,lits;
241   PetscBool          lssucceed;
242   PetscReal          gnorm,xnorm=0,ynorm;
243   Vec                Y,X,F;
244   KSPConvergedReason kspreason;
245   DM                 dm;
246   DMSNES             sdm;
247 
248   PetscFunctionBegin;
249   ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr);
250   ierr = DMGetDMSNES(dm,&sdm);CHKERRQ(ierr);
251 
252   vi->computeuserfunction   = sdm->ops->computefunction;
253   sdm->ops->computefunction = SNESVIComputeFunction;
254 
255   snes->numFailures            = 0;
256   snes->numLinearSolveFailures = 0;
257   snes->reason                 = SNES_CONVERGED_ITERATING;
258 
259   maxits = snes->max_its;               /* maximum number of iterations */
260   X      = snes->vec_sol;               /* solution vector */
261   F      = snes->vec_func;              /* residual vector */
262   Y      = snes->work[0];               /* work vectors */
263 
264   ierr       = PetscObjectSAWsTakeAccess((PetscObject)snes);CHKERRQ(ierr);
265   snes->iter = 0;
266   snes->norm = 0.0;
267   ierr       = PetscObjectSAWsGrantAccess((PetscObject)snes);CHKERRQ(ierr);
268 
269   ierr = SNESVIProjectOntoBounds(snes,X);CHKERRQ(ierr);
270   ierr = SNESComputeFunction(snes,X,vi->phi);CHKERRQ(ierr);
271   if (snes->domainerror) {
272     snes->reason              = SNES_DIVERGED_FUNCTION_DOMAIN;
273     sdm->ops->computefunction = vi->computeuserfunction;
274     PetscFunctionReturn(0);
275   }
276   /* Compute Merit function */
277   ierr = SNESVIComputeMeritFunction(vi->phi,&vi->merit,&vi->phinorm);CHKERRQ(ierr);
278 
279   ierr = VecNormBegin(X,NORM_2,&xnorm);CHKERRQ(ierr);        /* xnorm <- ||x||  */
280   ierr = VecNormEnd(X,NORM_2,&xnorm);CHKERRQ(ierr);
281   if (PetscIsInfOrNanReal(vi->merit)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number");
282 
283   ierr       = PetscObjectSAWsTakeAccess((PetscObject)snes);CHKERRQ(ierr);
284   snes->norm = vi->phinorm;
285   ierr       = PetscObjectSAWsGrantAccess((PetscObject)snes);CHKERRQ(ierr);
286   ierr       = SNESLogConvergenceHistory(snes,vi->phinorm,0);CHKERRQ(ierr);
287   ierr       = SNESMonitor(snes,0,vi->phinorm);CHKERRQ(ierr);
288 
289   /* test convergence */
290   ierr = (*snes->ops->converged)(snes,0,0.0,0.0,vi->phinorm,&snes->reason,snes->cnvP);CHKERRQ(ierr);
291   if (snes->reason) {
292     sdm->ops->computefunction = vi->computeuserfunction;
293     PetscFunctionReturn(0);
294   }
295 
296   for (i=0; i<maxits; i++) {
297 
298     /* Call general purpose update function */
299     if (snes->ops->update) {
300       ierr = (*snes->ops->update)(snes, snes->iter);CHKERRQ(ierr);
301     }
302 
303     /* Solve J Y = Phi, where J is the semismooth jacobian */
304 
305     /* Get the jacobian -- note that the function must be the original function for snes_fd and snes_fd_color to work for this*/
306     sdm->ops->computefunction = vi->computeuserfunction;
307     ierr                      = SNESComputeJacobian(snes,X,snes->jacobian,snes->jacobian_pre);CHKERRQ(ierr);
308     sdm->ops->computefunction = SNESVIComputeFunction;
309 
310     /* Get the diagonal shift and row scaling vectors */
311     ierr = SNESVIComputeBsubdifferentialVectors(snes,X,F,snes->jacobian,vi->Da,vi->Db);CHKERRQ(ierr);
312     /* Compute the semismooth jacobian */
313     ierr = SNESVIComputeJacobian(snes->jacobian,snes->jacobian_pre,vi->Da,vi->Db);CHKERRQ(ierr);
314     /* Compute the merit function gradient */
315     ierr = SNESVIComputeMeritFunctionGradient(snes->jacobian,vi->phi,vi->dpsi);CHKERRQ(ierr);
316     ierr = KSPSetOperators(snes->ksp,snes->jacobian,snes->jacobian_pre);CHKERRQ(ierr);
317     ierr = KSPSolve(snes->ksp,vi->phi,Y);CHKERRQ(ierr);
318     ierr = KSPGetConvergedReason(snes->ksp,&kspreason);CHKERRQ(ierr);
319 
320     if (kspreason < 0) {
321       if (++snes->numLinearSolveFailures >= snes->maxLinearSolveFailures) {
322         ierr         = PetscInfo2(snes,"iter=%D, number linear solve failures %D greater than current SNES allowed, stopping solve\n",snes->iter,snes->numLinearSolveFailures);CHKERRQ(ierr);
323         snes->reason = SNES_DIVERGED_LINEAR_SOLVE;
324         break;
325       }
326     }
327     ierr              = KSPGetIterationNumber(snes->ksp,&lits);CHKERRQ(ierr);
328     snes->linear_its += lits;
329     ierr              = PetscInfo2(snes,"iter=%D, linear solve iterations=%D\n",snes->iter,lits);CHKERRQ(ierr);
330     /*
331     if (snes->ops->precheck) {
332       PetscBool changed_y = PETSC_FALSE;
333       ierr = (*snes->ops->precheck)(snes,X,Y,snes->precheck,&changed_y);CHKERRQ(ierr);
334     }
335 
336     if (PetscLogPrintInfo) {
337       ierr = SNESVICheckResidual_Private(snes,snes->jacobian,F,Y,G,W);CHKERRQ(ierr);
338     }
339     */
340     /* Compute a (scaled) negative update in the line search routine:
341          Y <- X - lambda*Y
342        and evaluate G = function(Y) (depends on the line search).
343     */
344     ierr  = VecCopy(Y,snes->vec_sol_update);CHKERRQ(ierr);
345     ynorm = 1; gnorm = vi->phinorm;
346     ierr = SNESLineSearchApply(snes->linesearch, X, vi->phi, &gnorm, Y);CHKERRQ(ierr);
347     ierr = SNESLineSearchGetSuccess(snes->linesearch, &lssucceed);CHKERRQ(ierr);
348     ierr = SNESLineSearchGetNorms(snes->linesearch, &xnorm, &gnorm, &ynorm);CHKERRQ(ierr);
349     ierr = PetscInfo4(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lssucceed=%d\n",(double)vi->phinorm,(double)gnorm,(double)ynorm,(int)lssucceed);CHKERRQ(ierr);
350     if (snes->reason == SNES_DIVERGED_FUNCTION_COUNT) break;
351     if (snes->domainerror) {
352       snes->reason              = SNES_DIVERGED_FUNCTION_DOMAIN;
353       sdm->ops->computefunction = vi->computeuserfunction;
354       PetscFunctionReturn(0);
355     }
356     if (!lssucceed) {
357       if (++snes->numFailures >= snes->maxFailures) {
358         PetscBool ismin;
359         snes->reason = SNES_DIVERGED_LINE_SEARCH;
360         ierr         = SNESVICheckLocalMin_Private(snes,snes->jacobian,vi->phi,X,gnorm,&ismin);CHKERRQ(ierr);
361         if (ismin) snes->reason = SNES_DIVERGED_LOCAL_MIN;
362         break;
363       }
364     }
365     /* Update function and solution vectors */
366     vi->phinorm = gnorm;
367     vi->merit   = 0.5*vi->phinorm*vi->phinorm;
368     /* Monitor convergence */
369     ierr       = PetscObjectSAWsTakeAccess((PetscObject)snes);CHKERRQ(ierr);
370     snes->iter = i+1;
371     snes->norm = vi->phinorm;
372     ierr       = PetscObjectSAWsGrantAccess((PetscObject)snes);CHKERRQ(ierr);
373     ierr       = SNESLogConvergenceHistory(snes,snes->norm,lits);CHKERRQ(ierr);
374     ierr       = SNESMonitor(snes,snes->iter,snes->norm);CHKERRQ(ierr);
375     /* Test for convergence, xnorm = || X || */
376     if (snes->ops->converged != SNESConvergedSkip) { ierr = VecNorm(X,NORM_2,&xnorm);CHKERRQ(ierr); }
377     ierr = (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,vi->phinorm,&snes->reason,snes->cnvP);CHKERRQ(ierr);
378     if (snes->reason) break;
379   }
380   if (i == maxits) {
381     ierr = PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",maxits);CHKERRQ(ierr);
382     if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
383   }
384   sdm->ops->computefunction = vi->computeuserfunction;
385   PetscFunctionReturn(0);
386 }
387 
388 /* -------------------------------------------------------------------------- */
389 /*
390    SNESSetUp_VINEWTONSSLS - Sets up the internal data structures for the later use
391    of the SNES nonlinear solver.
392 
393    Input Parameter:
394 .  snes - the SNES context
395 .  x - the solution vector
396 
397    Application Interface Routine: SNESSetUp()
398 
399    Notes:
400    For basic use of the SNES solvers, the user need not explicitly call
401    SNESSetUp(), since these actions will automatically occur during
402    the call to SNESSolve().
403  */
404 #undef __FUNCT__
405 #define __FUNCT__ "SNESSetUp_VINEWTONSSLS"
406 PetscErrorCode SNESSetUp_VINEWTONSSLS(SNES snes)
407 {
408   PetscErrorCode    ierr;
409   SNES_VINEWTONSSLS *vi = (SNES_VINEWTONSSLS*) snes->data;
410 
411   PetscFunctionBegin;
412   ierr = SNESSetUp_VI(snes);CHKERRQ(ierr);
413   ierr = VecDuplicate(snes->vec_sol, &vi->dpsi);CHKERRQ(ierr);
414   ierr = VecDuplicate(snes->vec_sol, &vi->phi);CHKERRQ(ierr);
415   ierr = VecDuplicate(snes->vec_sol, &vi->Da);CHKERRQ(ierr);
416   ierr = VecDuplicate(snes->vec_sol, &vi->Db);CHKERRQ(ierr);
417   ierr = VecDuplicate(snes->vec_sol, &vi->z);CHKERRQ(ierr);
418   ierr = VecDuplicate(snes->vec_sol, &vi->t);CHKERRQ(ierr);
419   PetscFunctionReturn(0);
420 }
421 /* -------------------------------------------------------------------------- */
422 #undef __FUNCT__
423 #define __FUNCT__ "SNESReset_VINEWTONSSLS"
424 PetscErrorCode SNESReset_VINEWTONSSLS(SNES snes)
425 {
426   SNES_VINEWTONSSLS *vi = (SNES_VINEWTONSSLS*) snes->data;
427   PetscErrorCode    ierr;
428 
429   PetscFunctionBegin;
430   ierr = SNESReset_VI(snes);CHKERRQ(ierr);
431   ierr = VecDestroy(&vi->dpsi);CHKERRQ(ierr);
432   ierr = VecDestroy(&vi->phi);CHKERRQ(ierr);
433   ierr = VecDestroy(&vi->Da);CHKERRQ(ierr);
434   ierr = VecDestroy(&vi->Db);CHKERRQ(ierr);
435   ierr = VecDestroy(&vi->z);CHKERRQ(ierr);
436   ierr = VecDestroy(&vi->t);CHKERRQ(ierr);
437   PetscFunctionReturn(0);
438 }
439 
440 /* -------------------------------------------------------------------------- */
441 /*
442    SNESSetFromOptions_VINEWTONSSLS - Sets various parameters for the SNESVI method.
443 
444    Input Parameter:
445 .  snes - the SNES context
446 
447    Application Interface Routine: SNESSetFromOptions()
448 */
449 #undef __FUNCT__
450 #define __FUNCT__ "SNESSetFromOptions_VINEWTONSSLS"
451 static PetscErrorCode SNESSetFromOptions_VINEWTONSSLS(SNES snes)
452 {
453   PetscErrorCode ierr;
454   SNESLineSearch linesearch;
455 
456   PetscFunctionBegin;
457   ierr = SNESSetFromOptions_VI(snes);CHKERRQ(ierr);
458   ierr = PetscOptionsHead("SNES semismooth method options");CHKERRQ(ierr);
459   ierr = PetscOptionsTail();CHKERRQ(ierr);
460   /* set up the default line search */
461   if (!snes->linesearch) {
462     ierr = SNESGetLineSearch(snes, &linesearch);CHKERRQ(ierr);
463     ierr = SNESLineSearchSetType(linesearch, SNESLINESEARCHBT);CHKERRQ(ierr);
464     ierr = SNESLineSearchBTSetAlpha(linesearch, 0.0);CHKERRQ(ierr);
465   }
466   PetscFunctionReturn(0);
467 }
468 
469 
470 /* -------------------------------------------------------------------------- */
471 /*MC
472       SNESVINEWTONSSLS - Semi-smooth solver for variational inequalities based on Newton's method
473 
474    Options Database:
475 +   -snes_vi_type <ss,rs,rsaug> a semi-smooth solver, a reduced space active set method, and a reduced space active set method that does not eliminate the active constraints from the Jacobian instead augments the Jacobian with additional variables that enforce the constraints
476 -   -snes_vi_monitor - prints the number of active constraints at each iteration.
477 
478    Level: beginner
479 
480    References:
481    - T. S. Munson, F. Facchinei, M. C. Ferris, A. Fischer, and C. Kanzow. The semismooth
482      algorithm for large scale complementarity problems. INFORMS Journal on Computing, 13 (2001).
483 
484 .seealso:  SNESVISetVariableBounds(), SNESVISetComputeVariableBounds(), SNESCreate(), SNES, SNESSetType(), SNESVINEWTONRSLS, SNESVINEWTONSSLS, SNESNEWTONTR, SNESLineSearchSet(),
485            SNESLineSearchSetPostCheck(), SNESLineSearchNo(), SNESLineSearchCubic(), SNESLineSearchQuadratic(),
486            SNESLineSearchSet(), SNESLineSearchNoNorms(), SNESLineSearchSetPreCheck(), SNESLineSearchSetParams(), SNESLineSearchGetParams()
487 
488 M*/
489 #undef __FUNCT__
490 #define __FUNCT__ "SNESCreate_VINEWTONSSLS"
491 PETSC_EXTERN PetscErrorCode SNESCreate_VINEWTONSSLS(SNES snes)
492 {
493   PetscErrorCode    ierr;
494   SNES_VINEWTONSSLS *vi;
495 
496   PetscFunctionBegin;
497   snes->ops->reset          = SNESReset_VINEWTONSSLS;
498   snes->ops->setup          = SNESSetUp_VINEWTONSSLS;
499   snes->ops->solve          = SNESSolve_VINEWTONSSLS;
500   snes->ops->destroy        = SNESDestroy_VI;
501   snes->ops->setfromoptions = SNESSetFromOptions_VINEWTONSSLS;
502   snes->ops->view           = NULL;
503 
504   snes->usesksp = PETSC_TRUE;
505   snes->usespc  = PETSC_FALSE;
506 
507   ierr       = PetscNewLog(snes,&vi);CHKERRQ(ierr);
508   snes->data = (void*)vi;
509 
510   ierr = PetscObjectComposeFunction((PetscObject)snes,"SNESVISetVariableBounds_C",SNESVISetVariableBounds_VI);CHKERRQ(ierr);
511   ierr = PetscObjectComposeFunction((PetscObject)snes,"SNESVISetComputeVariableBounds_C",SNESVISetComputeVariableBounds_VI);CHKERRQ(ierr);
512   PetscFunctionReturn(0);
513 }
514 
515