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