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