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