xref: /petsc/src/snes/impls/vi/ss/viss.c (revision dcca6d9d80ebd869fe6029bd05a3aa9faafef49e)
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]) <= SNES_VI_NINF) && (PetscRealPart(u[i]) >= SNES_VI_INF)) { /* no constraints on variable */
80       phi_arr[i] = f_arr[i];
81     } else if (PetscRealPart(l[i]) <= SNES_VI_NINF) {                      /* upper bound on variable only */
82       phi_arr[i] = -Phi(u[i] - x_arr[i],-f_arr[i]);
83     } else if (PetscRealPart(u[i]) >= SNES_VI_INF) {                       /* 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]) <= SNES_VI_NINF) && (PetscRealPart(u[i]) >= SNES_VI_INF)) { /* no constraints on variable */
123       da[i] = 0;
124       db[i] = 1;
125     } else if (PetscRealPart(l[i]) <= SNES_VI_NINF) {                     /* 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]) >= SNES_VI_INF) {                      /* 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   MatStructure       flg = DIFFERENT_NONZERO_PATTERN;
243   PetscReal          gnorm,xnorm=0,ynorm;
244   Vec                Y,X,F;
245   KSPConvergedReason kspreason;
246   DM                 dm;
247   DMSNES             sdm;
248 
249   PetscFunctionBegin;
250   ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr);
251   ierr = DMGetDMSNES(dm,&sdm);CHKERRQ(ierr);
252 
253   vi->computeuserfunction   = sdm->ops->computefunction;
254   sdm->ops->computefunction = SNESVIComputeFunction;
255 
256   snes->numFailures            = 0;
257   snes->numLinearSolveFailures = 0;
258   snes->reason                 = SNES_CONVERGED_ITERATING;
259 
260   maxits = snes->max_its;               /* maximum number of iterations */
261   X      = snes->vec_sol;               /* solution vector */
262   F      = snes->vec_func;              /* residual vector */
263   Y      = snes->work[0];               /* work vectors */
264 
265   ierr       = PetscObjectSAWsTakeAccess((PetscObject)snes);CHKERRQ(ierr);
266   snes->iter = 0;
267   snes->norm = 0.0;
268   ierr       = PetscObjectSAWsGrantAccess((PetscObject)snes);CHKERRQ(ierr);
269 
270   ierr = SNESVIProjectOntoBounds(snes,X);CHKERRQ(ierr);
271   ierr = SNESComputeFunction(snes,X,vi->phi);CHKERRQ(ierr);
272   if (snes->domainerror) {
273     snes->reason              = SNES_DIVERGED_FUNCTION_DOMAIN;
274     sdm->ops->computefunction = vi->computeuserfunction;
275     PetscFunctionReturn(0);
276   }
277   /* Compute Merit function */
278   ierr = SNESVIComputeMeritFunction(vi->phi,&vi->merit,&vi->phinorm);CHKERRQ(ierr);
279 
280   ierr = VecNormBegin(X,NORM_2,&xnorm);CHKERRQ(ierr);        /* xnorm <- ||x||  */
281   ierr = VecNormEnd(X,NORM_2,&xnorm);CHKERRQ(ierr);
282   if (PetscIsInfOrNanReal(vi->merit)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number");
283 
284   ierr       = PetscObjectSAWsTakeAccess((PetscObject)snes);CHKERRQ(ierr);
285   snes->norm = vi->phinorm;
286   ierr       = PetscObjectSAWsGrantAccess((PetscObject)snes);CHKERRQ(ierr);
287   ierr       = SNESLogConvergenceHistory(snes,vi->phinorm,0);CHKERRQ(ierr);
288   ierr       = SNESMonitor(snes,0,vi->phinorm);CHKERRQ(ierr);
289 
290   /* test convergence */
291   ierr = (*snes->ops->converged)(snes,0,0.0,0.0,vi->phinorm,&snes->reason,snes->cnvP);CHKERRQ(ierr);
292   if (snes->reason) {
293     sdm->ops->computefunction = vi->computeuserfunction;
294     PetscFunctionReturn(0);
295   }
296 
297   for (i=0; i<maxits; i++) {
298 
299     /* Call general purpose update function */
300     if (snes->ops->update) {
301       ierr = (*snes->ops->update)(snes, snes->iter);CHKERRQ(ierr);
302     }
303 
304     /* Solve J Y = Phi, where J is the semismooth jacobian */
305 
306     /* Get the jacobian -- note that the function must be the original function for snes_fd and snes_fd_color to work for this*/
307     sdm->ops->computefunction = vi->computeuserfunction;
308     ierr                      = SNESComputeJacobian(snes,X,&snes->jacobian,&snes->jacobian_pre,&flg);CHKERRQ(ierr);
309     sdm->ops->computefunction = SNESVIComputeFunction;
310 
311     /* Get the diagonal shift and row scaling vectors */
312     ierr = SNESVIComputeBsubdifferentialVectors(snes,X,F,snes->jacobian,vi->Da,vi->Db);CHKERRQ(ierr);
313     /* Compute the semismooth jacobian */
314     ierr = SNESVIComputeJacobian(snes->jacobian,snes->jacobian_pre,vi->Da,vi->Db);CHKERRQ(ierr);
315     /* Compute the merit function gradient */
316     ierr = SNESVIComputeMeritFunctionGradient(snes->jacobian,vi->phi,vi->dpsi);CHKERRQ(ierr);
317     ierr = KSPSetOperators(snes->ksp,snes->jacobian,snes->jacobian_pre,flg);CHKERRQ(ierr);
318     ierr = KSPSolve(snes->ksp,vi->phi,Y);CHKERRQ(ierr);
319     ierr = KSPGetConvergedReason(snes->ksp,&kspreason);CHKERRQ(ierr);
320 
321     if (kspreason < 0) {
322       if (++snes->numLinearSolveFailures >= snes->maxLinearSolveFailures) {
323         ierr         = PetscInfo2(snes,"iter=%D, number linear solve failures %D greater than current SNES allowed, stopping solve\n",snes->iter,snes->numLinearSolveFailures);CHKERRQ(ierr);
324         snes->reason = SNES_DIVERGED_LINEAR_SOLVE;
325         break;
326       }
327     }
328     ierr              = KSPGetIterationNumber(snes->ksp,&lits);CHKERRQ(ierr);
329     snes->linear_its += lits;
330     ierr              = PetscInfo2(snes,"iter=%D, linear solve iterations=%D\n",snes->iter,lits);CHKERRQ(ierr);
331     /*
332     if (snes->ops->precheck) {
333       PetscBool changed_y = PETSC_FALSE;
334       ierr = (*snes->ops->precheck)(snes,X,Y,snes->precheck,&changed_y);CHKERRQ(ierr);
335     }
336 
337     if (PetscLogPrintInfo) {
338       ierr = SNESVICheckResidual_Private(snes,snes->jacobian,F,Y,G,W);CHKERRQ(ierr);
339     }
340     */
341     /* Compute a (scaled) negative update in the line search routine:
342          Y <- X - lambda*Y
343        and evaluate G = function(Y) (depends on the line search).
344     */
345     ierr  = VecCopy(Y,snes->vec_sol_update);CHKERRQ(ierr);
346     ynorm = 1; gnorm = vi->phinorm;
347     ierr = SNESLineSearchApply(snes->linesearch, X, vi->phi, &gnorm, Y);CHKERRQ(ierr);
348     ierr = SNESLineSearchGetSuccess(snes->linesearch, &lssucceed);CHKERRQ(ierr);
349     ierr = SNESLineSearchGetNorms(snes->linesearch, &xnorm, &gnorm, &ynorm);CHKERRQ(ierr);
350     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);
351     if (snes->reason == SNES_DIVERGED_FUNCTION_COUNT) break;
352     if (snes->domainerror) {
353       snes->reason              = SNES_DIVERGED_FUNCTION_DOMAIN;
354       sdm->ops->computefunction = vi->computeuserfunction;
355       PetscFunctionReturn(0);
356     }
357     if (!lssucceed) {
358       if (++snes->numFailures >= snes->maxFailures) {
359         PetscBool ismin;
360         snes->reason = SNES_DIVERGED_LINE_SEARCH;
361         ierr         = SNESVICheckLocalMin_Private(snes,snes->jacobian,vi->phi,X,gnorm,&ismin);CHKERRQ(ierr);
362         if (ismin) snes->reason = SNES_DIVERGED_LOCAL_MIN;
363         break;
364       }
365     }
366     /* Update function and solution vectors */
367     vi->phinorm = gnorm;
368     vi->merit   = 0.5*vi->phinorm*vi->phinorm;
369     /* Monitor convergence */
370     ierr       = PetscObjectSAWsTakeAccess((PetscObject)snes);CHKERRQ(ierr);
371     snes->iter = i+1;
372     snes->norm = vi->phinorm;
373     ierr       = PetscObjectSAWsGrantAccess((PetscObject)snes);CHKERRQ(ierr);
374     ierr       = SNESLogConvergenceHistory(snes,snes->norm,lits);CHKERRQ(ierr);
375     ierr       = SNESMonitor(snes,snes->iter,snes->norm);CHKERRQ(ierr);
376     /* Test for convergence, xnorm = || X || */
377     if (snes->ops->converged != SNESConvergedSkip) { ierr = VecNorm(X,NORM_2,&xnorm);CHKERRQ(ierr); }
378     ierr = (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,vi->phinorm,&snes->reason,snes->cnvP);CHKERRQ(ierr);
379     if (snes->reason) break;
380   }
381   if (i == maxits) {
382     ierr = PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",maxits);CHKERRQ(ierr);
383     if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
384   }
385   sdm->ops->computefunction = vi->computeuserfunction;
386   PetscFunctionReturn(0);
387 }
388 
389 /* -------------------------------------------------------------------------- */
390 /*
391    SNESSetUp_VINEWTONSSLS - Sets up the internal data structures for the later use
392    of the SNES nonlinear solver.
393 
394    Input Parameter:
395 .  snes - the SNES context
396 .  x - the solution vector
397 
398    Application Interface Routine: SNESSetUp()
399 
400    Notes:
401    For basic use of the SNES solvers, the user need not explicitly call
402    SNESSetUp(), since these actions will automatically occur during
403    the call to SNESSolve().
404  */
405 #undef __FUNCT__
406 #define __FUNCT__ "SNESSetUp_VINEWTONSSLS"
407 PetscErrorCode SNESSetUp_VINEWTONSSLS(SNES snes)
408 {
409   PetscErrorCode    ierr;
410   SNES_VINEWTONSSLS *vi = (SNES_VINEWTONSSLS*) snes->data;
411 
412   PetscFunctionBegin;
413   ierr = SNESSetUp_VI(snes);CHKERRQ(ierr);
414   ierr = VecDuplicate(snes->vec_sol, &vi->dpsi);CHKERRQ(ierr);
415   ierr = VecDuplicate(snes->vec_sol, &vi->phi);CHKERRQ(ierr);
416   ierr = VecDuplicate(snes->vec_sol, &vi->Da);CHKERRQ(ierr);
417   ierr = VecDuplicate(snes->vec_sol, &vi->Db);CHKERRQ(ierr);
418   ierr = VecDuplicate(snes->vec_sol, &vi->z);CHKERRQ(ierr);
419   ierr = VecDuplicate(snes->vec_sol, &vi->t);CHKERRQ(ierr);
420   PetscFunctionReturn(0);
421 }
422 /* -------------------------------------------------------------------------- */
423 #undef __FUNCT__
424 #define __FUNCT__ "SNESReset_VINEWTONSSLS"
425 PetscErrorCode SNESReset_VINEWTONSSLS(SNES snes)
426 {
427   SNES_VINEWTONSSLS *vi = (SNES_VINEWTONSSLS*) snes->data;
428   PetscErrorCode    ierr;
429 
430   PetscFunctionBegin;
431   ierr = SNESReset_VI(snes);CHKERRQ(ierr);
432   ierr = VecDestroy(&vi->dpsi);CHKERRQ(ierr);
433   ierr = VecDestroy(&vi->phi);CHKERRQ(ierr);
434   ierr = VecDestroy(&vi->Da);CHKERRQ(ierr);
435   ierr = VecDestroy(&vi->Db);CHKERRQ(ierr);
436   ierr = VecDestroy(&vi->z);CHKERRQ(ierr);
437   ierr = VecDestroy(&vi->t);CHKERRQ(ierr);
438   PetscFunctionReturn(0);
439 }
440 
441 /* -------------------------------------------------------------------------- */
442 /*
443    SNESSetFromOptions_VINEWTONSSLS - Sets various parameters for the SNESVI method.
444 
445    Input Parameter:
446 .  snes - the SNES context
447 
448    Application Interface Routine: SNESSetFromOptions()
449 */
450 #undef __FUNCT__
451 #define __FUNCT__ "SNESSetFromOptions_VINEWTONSSLS"
452 static PetscErrorCode SNESSetFromOptions_VINEWTONSSLS(SNES snes)
453 {
454   PetscErrorCode ierr;
455   SNESLineSearch linesearch;
456 
457   PetscFunctionBegin;
458   ierr = SNESSetFromOptions_VI(snes);CHKERRQ(ierr);
459   ierr = PetscOptionsHead("SNES semismooth method options");CHKERRQ(ierr);
460   ierr = PetscOptionsTail();CHKERRQ(ierr);
461   /* set up the default line search */
462   if (!snes->linesearch) {
463     ierr = SNESGetLineSearch(snes, &linesearch);CHKERRQ(ierr);
464     ierr = SNESLineSearchSetType(linesearch, SNESLINESEARCHBT);CHKERRQ(ierr);
465     ierr = SNESLineSearchBTSetAlpha(linesearch, 0.0);CHKERRQ(ierr);
466   }
467   PetscFunctionReturn(0);
468 }
469 
470 
471 /* -------------------------------------------------------------------------- */
472 /*MC
473       SNESVINEWTONSSLS - Semi-smooth solver for variational inequalities based on Newton's method
474 
475    Options Database:
476 +   -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
477 -   -snes_vi_monitor - prints the number of active constraints at each iteration.
478 
479    Level: beginner
480 
481    References:
482    - T. S. Munson, F. Facchinei, M. C. Ferris, A. Fischer, and C. Kanzow. The semismooth
483      algorithm for large scale complementarity problems. INFORMS Journal on Computing, 13 (2001).
484 
485 .seealso:  SNESVISetVariableBounds(), SNESVISetComputeVariableBounds(), SNESCreate(), SNES, SNESSetType(), SNESVINEWTONRSLS, SNESVINEWTONSSLS, SNESNEWTONTR, SNESLineSearchSet(),
486            SNESLineSearchSetPostCheck(), SNESLineSearchNo(), SNESLineSearchCubic(), SNESLineSearchQuadratic(),
487            SNESLineSearchSet(), SNESLineSearchNoNorms(), SNESLineSearchSetPreCheck(), SNESLineSearchSetParams(), SNESLineSearchGetParams()
488 
489 M*/
490 #undef __FUNCT__
491 #define __FUNCT__ "SNESCreate_VINEWTONSSLS"
492 PETSC_EXTERN PetscErrorCode SNESCreate_VINEWTONSSLS(SNES snes)
493 {
494   PetscErrorCode    ierr;
495   SNES_VINEWTONSSLS *vi;
496 
497   PetscFunctionBegin;
498   snes->ops->reset          = SNESReset_VINEWTONSSLS;
499   snes->ops->setup          = SNESSetUp_VINEWTONSSLS;
500   snes->ops->solve          = SNESSolve_VINEWTONSSLS;
501   snes->ops->destroy        = SNESDestroy_VI;
502   snes->ops->setfromoptions = SNESSetFromOptions_VINEWTONSSLS;
503   snes->ops->view           = NULL;
504 
505   snes->usesksp = PETSC_TRUE;
506   snes->usespc  = PETSC_FALSE;
507 
508   ierr       = PetscNewLog(snes,SNES_VINEWTONSSLS,&vi);CHKERRQ(ierr);
509   snes->data = (void*)vi;
510 
511   ierr = PetscObjectComposeFunction((PetscObject)snes,"SNESVISetVariableBounds_C",SNESVISetVariableBounds_VI);CHKERRQ(ierr);
512   ierr = PetscObjectComposeFunction((PetscObject)snes,"SNESVISetComputeVariableBounds_C",SNESVISetComputeVariableBounds_VI);CHKERRQ(ierr);
513   PetscFunctionReturn(0);
514 }
515 
516