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