xref: /petsc/src/snes/impls/vi/ss/viss.c (revision 4f02bc6a7df4e4b03c783003a8e0899ee35fcb05)
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(PetscOptionItems *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