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