1 2 #include <../src/snes/impls/vi/viimpl.h> /*I "petscsnes.h" I*/ 3 #include <../include/private/kspimpl.h> 4 #include <../include/private/matimpl.h> 5 #include <../include/private/dmimpl.h> 6 7 typedef struct { 8 IS inactive; 9 Vec upper,lower,values,F; /* upper and lower bounds of all variables on this level, the values and the function values */ 10 PetscErrorCode (*getinterpolation)(DM,DM,Mat*,Vec*); /* DM's original routines */ 11 PetscErrorCode (*coarsen)(DM, MPI_Comm, DM*); 12 } DMSNESVI; 13 14 #undef __FUNCT__ 15 #define __FUNCT__ "SNESVIGetInActiveSetIS" 16 /* 17 SNESVIGetInActiveSetIS - Gets the global indices for the bogus inactive set variables 18 19 Input parameter 20 . snes - the SNES context 21 . X - the snes solution vector 22 23 Output parameter 24 . ISact - active set index set 25 26 */ 27 PetscErrorCode SNESVIGetInActiveSetIS(Vec upper,Vec lower,Vec X,Vec F,IS* inact) 28 { 29 PetscErrorCode ierr; 30 const PetscScalar *x,*xl,*xu,*f; 31 PetscInt *idx_act,i,nlocal,nloc_isact=0,ilow,ihigh,i1=0; 32 33 PetscFunctionBegin; 34 ierr = VecGetLocalSize(X,&nlocal);CHKERRQ(ierr); 35 ierr = VecGetOwnershipRange(X,&ilow,&ihigh);CHKERRQ(ierr); 36 ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr); 37 ierr = VecGetArrayRead(lower,&xl);CHKERRQ(ierr); 38 ierr = VecGetArrayRead(upper,&xu);CHKERRQ(ierr); 39 ierr = VecGetArrayRead(F,&f);CHKERRQ(ierr); 40 /* Compute inactive set size */ 41 for (i=0; i < nlocal;i++) { 42 if (((PetscRealPart(x[i]) > PetscRealPart(xl[i]) + 1.e-8 || (PetscRealPart(f[i]) < 0.0)) && ((PetscRealPart(x[i]) < PetscRealPart(xu[i]) - 1.e-8) || PetscRealPart(f[i]) > 0.0))) nloc_isact++; 43 } 44 45 ierr = PetscMalloc(nloc_isact*sizeof(PetscInt),&idx_act);CHKERRQ(ierr); 46 47 /* Set inactive set indices */ 48 for(i=0; i < nlocal; i++) { 49 if (((PetscRealPart(x[i]) > PetscRealPart(xl[i]) + 1.e-8 || (PetscRealPart(f[i]) < 0.0)) && ((PetscRealPart(x[i]) < PetscRealPart(xu[i]) - 1.e-8) || PetscRealPart(f[i]) > 0.0))) idx_act[i1++] = ilow+i; 50 } 51 52 /* Create inactive set IS */ 53 ierr = ISCreateGeneral(((PetscObject)upper)->comm,nloc_isact,idx_act,PETSC_OWN_POINTER,inact);CHKERRQ(ierr); 54 55 ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr); 56 ierr = VecRestoreArrayRead(lower,&xl);CHKERRQ(ierr); 57 ierr = VecRestoreArrayRead(upper,&xu);CHKERRQ(ierr); 58 ierr = VecRestoreArrayRead(F,&f);CHKERRQ(ierr); 59 PetscFunctionReturn(0); 60 } 61 62 #undef __FUNCT__ 63 #define __FUNCT__ "DMGetInterpolation_SNESVI" 64 /* 65 DMGetInterpolation_SNESVI - Modifieds the interpolation obtained from the DM by removing all rows and columns associated with active constraints. 66 67 */ 68 PetscErrorCode DMGetInterpolation_SNESVI(DM dm1,DM dm2,Mat *mat,Vec *vec) 69 { 70 PetscErrorCode ierr; 71 PetscContainer isnes; 72 DMSNESVI *dmsnesvi1,*dmsnesvi2; 73 Mat interp; 74 75 PetscFunctionBegin; 76 ierr = PetscObjectQuery((PetscObject)dm1,"VI",(PetscObject *)&isnes);CHKERRQ(ierr); 77 if (isnes) SETERRQ(((PetscObject)dm1)->comm,PETSC_ERR_PLIB,"Composed SNES is missing"); 78 ierr = PetscContainerGetPointer(isnes,(void**)&dmsnesvi1);CHKERRQ(ierr); 79 ierr = PetscObjectQuery((PetscObject)dm2,"VI",(PetscObject *)&isnes);CHKERRQ(ierr); 80 if (isnes) SETERRQ(((PetscObject)dm2)->comm,PETSC_ERR_PLIB,"Composed SNES is missing"); 81 ierr = PetscContainerGetPointer(isnes,(void**)&dmsnesvi2);CHKERRQ(ierr); 82 83 ierr = (*dmsnesvi1->getinterpolation)(dm1,dm2,&interp,PETSC_NULL);CHKERRQ(ierr); 84 ierr = MatGetSubMatrix(interp,dmsnesvi1->inactive,dmsnesvi2->inactive,MAT_INITIAL_MATRIX,mat);CHKERRQ(ierr); 85 ierr = MatDestroy(&interp);CHKERRQ(ierr); 86 PetscFunctionReturn(0); 87 } 88 89 extern PetscErrorCode DMSetVI(DM,Vec,Vec,Vec,Vec,IS); 90 91 #undef __FUNCT__ 92 #define __FUNCT__ "DMCoarsen_SNESVI" 93 /* 94 DMCoarsen_SNESVI - Computes the regular coarsened DM then computes additional information about its inactive set 95 96 */ 97 PetscErrorCode DMCoarsen_SNESVI(DM dm1,MPI_Comm comm,DM *dm2) 98 { 99 PetscErrorCode ierr; 100 PetscContainer isnes; 101 DMSNESVI *dmsnesvi1; 102 Vec upper,lower,values,F; 103 IS inactive; 104 VecScatter inject; 105 106 PetscFunctionBegin; 107 ierr = PetscObjectQuery((PetscObject)dm1,"VI",(PetscObject *)&isnes);CHKERRQ(ierr); 108 if (isnes) SETERRQ(((PetscObject)dm1)->comm,PETSC_ERR_PLIB,"Composed SNES is missing"); 109 ierr = PetscContainerGetPointer(isnes,(void**)&dmsnesvi1);CHKERRQ(ierr); 110 111 ierr = (*dmsnesvi1->coarsen)(dm1,comm,dm2);CHKERRQ(ierr); 112 ierr = DMGetInjection(*dm2,dm1,&inject);CHKERRQ(ierr); 113 ierr = DMCreateGlobalVector(*dm2,&upper);CHKERRQ(ierr); 114 ierr = VecScatterBegin(inject,dmsnesvi1->upper,upper,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 115 ierr = VecScatterEnd(inject,dmsnesvi1->upper,upper,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 116 ierr = DMCreateGlobalVector(*dm2,&lower);CHKERRQ(ierr); 117 ierr = VecScatterBegin(inject,dmsnesvi1->lower,lower,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 118 ierr = VecScatterEnd(inject,dmsnesvi1->lower,lower,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 119 ierr = DMCreateGlobalVector(*dm2,&values);CHKERRQ(ierr); 120 ierr = VecScatterBegin(inject,dmsnesvi1->values,values,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 121 ierr = VecScatterEnd(inject,dmsnesvi1->values,values,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 122 ierr = DMCreateGlobalVector(*dm2,&F);CHKERRQ(ierr); 123 ierr = VecScatterBegin(inject,dmsnesvi1->F,F,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 124 ierr = VecScatterEnd(inject,dmsnesvi1->F,F,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 125 ierr = VecScatterDestroy(&inject);CHKERRQ(ierr); 126 ierr = SNESVIGetInActiveSetIS(upper,lower,values,F,&inactive);CHKERRQ(ierr); 127 ierr = DMSetVI(*dm2,upper,lower,values,F,inactive);CHKERRQ(ierr); 128 PetscFunctionReturn(0); 129 } 130 131 #undef __FUNCT__ 132 #define __FUNCT__ "DMSNESVIDestroy" 133 PetscErrorCode DMSNESVIDestroy(DMSNESVI *dmsnesvi) 134 { 135 PetscErrorCode ierr; 136 137 PetscFunctionBegin; 138 ierr = VecDestroy(&dmsnesvi->upper);CHKERRQ(ierr); 139 ierr = VecDestroy(&dmsnesvi->lower);CHKERRQ(ierr); 140 ierr = VecDestroy(&dmsnesvi->values);CHKERRQ(ierr); 141 ierr = VecDestroy(&dmsnesvi->F);CHKERRQ(ierr); 142 ierr = ISDestroy(&dmsnesvi->inactive);CHKERRQ(ierr); 143 ierr = PetscFree(dmsnesvi);CHKERRQ(ierr); 144 PetscFunctionReturn(0); 145 } 146 147 #undef __FUNCT__ 148 #define __FUNCT__ "DMSetVI" 149 /* 150 DMSetVI - Marks a DM as associated with a VI problem. This causes the interpolation/restriction operators to 151 be restricted to only those variables NOT associated with active constraints. 152 153 */ 154 PetscErrorCode DMSetVI(DM dm,Vec upper,Vec lower,Vec values,Vec F,IS inactive) 155 { 156 PetscErrorCode ierr; 157 PetscContainer isnes; 158 DMSNESVI *dmsnesvi; 159 160 PetscFunctionBegin; 161 ierr = PetscObjectReference((PetscObject)upper);CHKERRQ(ierr); 162 ierr = PetscObjectReference((PetscObject)lower);CHKERRQ(ierr); 163 ierr = PetscObjectReference((PetscObject)values);CHKERRQ(ierr); 164 ierr = PetscObjectReference((PetscObject)F);CHKERRQ(ierr); 165 ierr = PetscObjectReference((PetscObject)inactive);CHKERRQ(ierr); 166 167 ierr = PetscObjectQuery((PetscObject)dm,"VI",(PetscObject *)&isnes);CHKERRQ(ierr); 168 if (!isnes) { 169 /* cannot just compose snes into dm because that will cause circular reference */ 170 ierr = PetscContainerCreate(((PetscObject)dm)->comm,&isnes);CHKERRQ(ierr); 171 ierr = PetscContainerSetUserDestroy(isnes,(PetscErrorCode (*)(void*))DMSNESVIDestroy);CHKERRQ(ierr); 172 ierr = PetscNew(DMSNESVI,&dmsnesvi);CHKERRQ(ierr); 173 ierr = PetscContainerSetPointer(isnes,(void*)dmsnesvi);CHKERRQ(ierr); 174 ierr = PetscObjectCompose((PetscObject)dm,"VI",(PetscObject)isnes);CHKERRQ(ierr); 175 dmsnesvi->getinterpolation = dm->ops->getinterpolation; 176 dm->ops->getinterpolation = DMGetInterpolation_SNESVI; 177 dmsnesvi->coarsen = dm->ops->coarsen; 178 dm->ops->coarsen = DMCoarsen_SNESVI; 179 } else { 180 ierr = PetscContainerGetPointer(isnes,(void**)&dmsnesvi);CHKERRQ(ierr); 181 ierr = VecDestroy(&dmsnesvi->upper);CHKERRQ(ierr); 182 ierr = VecDestroy(&dmsnesvi->lower);CHKERRQ(ierr); 183 ierr = VecDestroy(&dmsnesvi->values);CHKERRQ(ierr); 184 ierr = VecDestroy(&dmsnesvi->F);CHKERRQ(ierr); 185 ierr = ISDestroy(&dmsnesvi->inactive);CHKERRQ(ierr); 186 } 187 dmsnesvi->upper = upper; 188 dmsnesvi->lower = lower; 189 dmsnesvi->values = values; 190 dmsnesvi->F = F; 191 dmsnesvi->inactive = inactive; 192 PetscFunctionReturn(0); 193 } 194 195 196 197 #undef __FUNCT__ 198 #define __FUNCT__ "SNESMonitorVI" 199 PetscErrorCode SNESMonitorVI(SNES snes,PetscInt its,PetscReal fgnorm,void *dummy) 200 { 201 PetscErrorCode ierr; 202 SNES_VI *vi = (SNES_VI*)snes->data; 203 PetscViewerASCIIMonitor viewer = (PetscViewerASCIIMonitor) dummy; 204 const PetscScalar *x,*xl,*xu,*f; 205 PetscInt i,n,act = 0; 206 PetscReal rnorm,fnorm; 207 208 PetscFunctionBegin; 209 if (!dummy) { 210 ierr = PetscViewerASCIIMonitorCreate(((PetscObject)snes)->comm,"stdout",0,&viewer);CHKERRQ(ierr); 211 } 212 ierr = VecGetLocalSize(snes->vec_sol,&n);CHKERRQ(ierr); 213 ierr = VecGetArrayRead(vi->xl,&xl);CHKERRQ(ierr); 214 ierr = VecGetArrayRead(vi->xu,&xu);CHKERRQ(ierr); 215 ierr = VecGetArrayRead(snes->vec_sol,&x);CHKERRQ(ierr); 216 ierr = VecGetArrayRead(snes->vec_func,&f);CHKERRQ(ierr); 217 218 rnorm = 0.0; 219 for (i=0; i<n; i++) { 220 if (((PetscRealPart(x[i]) > PetscRealPart(xl[i]) + 1.e-8 || (PetscRealPart(f[i]) < 0.0)) && ((PetscRealPart(x[i]) < PetscRealPart(xu[i]) - 1.e-8) || PetscRealPart(f[i]) > 0.0))) rnorm += PetscRealPart(PetscConj(f[i])*f[i]); 221 else act++; 222 } 223 ierr = VecRestoreArrayRead(snes->vec_func,&f);CHKERRQ(ierr); 224 ierr = VecRestoreArrayRead(vi->xl,&xl);CHKERRQ(ierr); 225 ierr = VecRestoreArrayRead(vi->xu,&xu);CHKERRQ(ierr); 226 ierr = VecRestoreArrayRead(snes->vec_sol,&x);CHKERRQ(ierr); 227 ierr = MPI_Allreduce(&rnorm,&fnorm,1,MPIU_REAL,MPIU_SUM,((PetscObject)snes)->comm);CHKERRQ(ierr); 228 fnorm = sqrt(fnorm); 229 ierr = PetscViewerASCIIMonitorPrintf(viewer,"%3D SNES VI Function norm %14.12e Active constraints %D\n",its,fnorm,act);CHKERRQ(ierr); 230 if (!dummy) { 231 ierr = PetscViewerASCIIMonitorDestroy(&viewer);CHKERRQ(ierr); 232 } 233 PetscFunctionReturn(0); 234 } 235 236 /* 237 Checks if J^T F = 0 which implies we've found a local minimum of the norm of the function, 238 || F(u) ||_2 but not a zero, F(u) = 0. In the case when one cannot compute J^T F we use the fact that 239 0 = (J^T F)^T W = F^T J W iff W not in the null space of J. Thanks for Jorge More 240 for this trick. One assumes that the probability that W is in the null space of J is very, very small. 241 */ 242 #undef __FUNCT__ 243 #define __FUNCT__ "SNESVICheckLocalMin_Private" 244 PetscErrorCode SNESVICheckLocalMin_Private(SNES snes,Mat A,Vec F,Vec W,PetscReal fnorm,PetscBool *ismin) 245 { 246 PetscReal a1; 247 PetscErrorCode ierr; 248 PetscBool hastranspose; 249 250 PetscFunctionBegin; 251 *ismin = PETSC_FALSE; 252 ierr = MatHasOperation(A,MATOP_MULT_TRANSPOSE,&hastranspose);CHKERRQ(ierr); 253 if (hastranspose) { 254 /* Compute || J^T F|| */ 255 ierr = MatMultTranspose(A,F,W);CHKERRQ(ierr); 256 ierr = VecNorm(W,NORM_2,&a1);CHKERRQ(ierr); 257 ierr = PetscInfo1(snes,"|| J^T F|| %G near zero implies found a local minimum\n",a1/fnorm);CHKERRQ(ierr); 258 if (a1/fnorm < 1.e-4) *ismin = PETSC_TRUE; 259 } else { 260 Vec work; 261 PetscScalar result; 262 PetscReal wnorm; 263 264 ierr = VecSetRandom(W,PETSC_NULL);CHKERRQ(ierr); 265 ierr = VecNorm(W,NORM_2,&wnorm);CHKERRQ(ierr); 266 ierr = VecDuplicate(W,&work);CHKERRQ(ierr); 267 ierr = MatMult(A,W,work);CHKERRQ(ierr); 268 ierr = VecDot(F,work,&result);CHKERRQ(ierr); 269 ierr = VecDestroy(&work);CHKERRQ(ierr); 270 a1 = PetscAbsScalar(result)/(fnorm*wnorm); 271 ierr = PetscInfo1(snes,"(F^T J random)/(|| F ||*||random|| %G near zero implies found a local minimum\n",a1);CHKERRQ(ierr); 272 if (a1 < 1.e-4) *ismin = PETSC_TRUE; 273 } 274 PetscFunctionReturn(0); 275 } 276 277 /* 278 Checks if J^T(F - J*X) = 0 279 */ 280 #undef __FUNCT__ 281 #define __FUNCT__ "SNESVICheckResidual_Private" 282 PetscErrorCode SNESVICheckResidual_Private(SNES snes,Mat A,Vec F,Vec X,Vec W1,Vec W2) 283 { 284 PetscReal a1,a2; 285 PetscErrorCode ierr; 286 PetscBool hastranspose; 287 288 PetscFunctionBegin; 289 ierr = MatHasOperation(A,MATOP_MULT_TRANSPOSE,&hastranspose);CHKERRQ(ierr); 290 if (hastranspose) { 291 ierr = MatMult(A,X,W1);CHKERRQ(ierr); 292 ierr = VecAXPY(W1,-1.0,F);CHKERRQ(ierr); 293 294 /* Compute || J^T W|| */ 295 ierr = MatMultTranspose(A,W1,W2);CHKERRQ(ierr); 296 ierr = VecNorm(W1,NORM_2,&a1);CHKERRQ(ierr); 297 ierr = VecNorm(W2,NORM_2,&a2);CHKERRQ(ierr); 298 if (a1 != 0.0) { 299 ierr = PetscInfo1(snes,"||J^T(F-Ax)||/||F-AX|| %G near zero implies inconsistent rhs\n",a2/a1);CHKERRQ(ierr); 300 } 301 } 302 PetscFunctionReturn(0); 303 } 304 305 /* 306 SNESDefaultConverged_VI - Checks the convergence of the semismooth newton algorithm. 307 308 Notes: 309 The convergence criterion currently implemented is 310 merit < abstol 311 merit < rtol*merit_initial 312 */ 313 #undef __FUNCT__ 314 #define __FUNCT__ "SNESDefaultConverged_VI" 315 PetscErrorCode SNESDefaultConverged_VI(SNES snes,PetscInt it,PetscReal xnorm,PetscReal gradnorm,PetscReal fnorm,SNESConvergedReason *reason,void *dummy) 316 { 317 PetscErrorCode ierr; 318 319 PetscFunctionBegin; 320 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 321 PetscValidPointer(reason,6); 322 323 *reason = SNES_CONVERGED_ITERATING; 324 325 if (!it) { 326 /* set parameter for default relative tolerance convergence test */ 327 snes->ttol = fnorm*snes->rtol; 328 } 329 if (fnorm != fnorm) { 330 ierr = PetscInfo(snes,"Failed to converged, function norm is NaN\n");CHKERRQ(ierr); 331 *reason = SNES_DIVERGED_FNORM_NAN; 332 } else if (fnorm < snes->abstol) { 333 ierr = PetscInfo2(snes,"Converged due to function norm %G < %G\n",fnorm,snes->abstol);CHKERRQ(ierr); 334 *reason = SNES_CONVERGED_FNORM_ABS; 335 } else if (snes->nfuncs >= snes->max_funcs) { 336 ierr = PetscInfo2(snes,"Exceeded maximum number of function evaluations: %D > %D\n",snes->nfuncs,snes->max_funcs);CHKERRQ(ierr); 337 *reason = SNES_DIVERGED_FUNCTION_COUNT; 338 } 339 340 if (it && !*reason) { 341 if (fnorm < snes->ttol) { 342 ierr = PetscInfo2(snes,"Converged due to function norm %G < %G (relative tolerance)\n",fnorm,snes->ttol);CHKERRQ(ierr); 343 *reason = SNES_CONVERGED_FNORM_RELATIVE; 344 } 345 } 346 PetscFunctionReturn(0); 347 } 348 349 /* 350 SNESVIComputeMeritFunction - Evaluates the merit function for the mixed complementarity problem. 351 352 Input Parameter: 353 . phi - the semismooth function 354 355 Output Parameter: 356 . merit - the merit function 357 . phinorm - ||phi|| 358 359 Notes: 360 The merit function for the mixed complementarity problem is defined as 361 merit = 0.5*phi^T*phi 362 */ 363 #undef __FUNCT__ 364 #define __FUNCT__ "SNESVIComputeMeritFunction" 365 static PetscErrorCode SNESVIComputeMeritFunction(Vec phi, PetscReal* merit,PetscReal* phinorm) 366 { 367 PetscErrorCode ierr; 368 369 PetscFunctionBegin; 370 ierr = VecNormBegin(phi,NORM_2,phinorm);CHKERRQ(ierr); 371 ierr = VecNormEnd(phi,NORM_2,phinorm);CHKERRQ(ierr); 372 373 *merit = 0.5*(*phinorm)*(*phinorm); 374 PetscFunctionReturn(0); 375 } 376 377 PETSC_STATIC_INLINE PetscScalar Phi(PetscScalar a,PetscScalar b) 378 { 379 return a + b - sqrt(a*a + b*b); 380 } 381 382 PETSC_STATIC_INLINE PetscScalar DPhi(PetscScalar a,PetscScalar b) 383 { 384 if ((PetscAbsScalar(a) >= 1.e-6) || (PetscAbsScalar(b) >= 1.e-6)) return 1.0 - a/ sqrt(a*a + b*b); 385 else return .5; 386 } 387 388 /* 389 SNESVIComputeFunction - Reformulates a system of nonlinear equations in mixed complementarity form to a system of nonlinear equations in semismooth form. 390 391 Input Parameters: 392 . snes - the SNES context 393 . x - current iterate 394 . functx - user defined function context 395 396 Output Parameters: 397 . phi - Semismooth function 398 399 */ 400 #undef __FUNCT__ 401 #define __FUNCT__ "SNESVIComputeFunction" 402 static PetscErrorCode SNESVIComputeFunction(SNES snes,Vec X,Vec phi,void* functx) 403 { 404 PetscErrorCode ierr; 405 SNES_VI *vi = (SNES_VI*)snes->data; 406 Vec Xl = vi->xl,Xu = vi->xu,F = snes->vec_func; 407 PetscScalar *phi_arr,*x_arr,*f_arr,*l,*u; 408 PetscInt i,nlocal; 409 410 PetscFunctionBegin; 411 ierr = (*vi->computeuserfunction)(snes,X,F,functx);CHKERRQ(ierr); 412 ierr = VecGetLocalSize(X,&nlocal);CHKERRQ(ierr); 413 ierr = VecGetArray(X,&x_arr);CHKERRQ(ierr); 414 ierr = VecGetArray(F,&f_arr);CHKERRQ(ierr); 415 ierr = VecGetArray(Xl,&l);CHKERRQ(ierr); 416 ierr = VecGetArray(Xu,&u);CHKERRQ(ierr); 417 ierr = VecGetArray(phi,&phi_arr);CHKERRQ(ierr); 418 419 for (i=0;i < nlocal;i++) { 420 if ((PetscRealPart(l[i]) <= SNES_VI_NINF) && (PetscRealPart(u[i]) >= SNES_VI_INF)) { /* no constraints on variable */ 421 phi_arr[i] = f_arr[i]; 422 } else if (PetscRealPart(l[i]) <= SNES_VI_NINF) { /* upper bound on variable only */ 423 phi_arr[i] = -Phi(u[i] - x_arr[i],-f_arr[i]); 424 } else if (PetscRealPart(u[i]) >= SNES_VI_INF) { /* lower bound on variable only */ 425 phi_arr[i] = Phi(x_arr[i] - l[i],f_arr[i]); 426 } else if (l[i] == u[i]) { 427 phi_arr[i] = l[i] - x_arr[i]; 428 } else { /* both bounds on variable */ 429 phi_arr[i] = Phi(x_arr[i] - l[i],-Phi(u[i] - x_arr[i],-f_arr[i])); 430 } 431 } 432 433 ierr = VecRestoreArray(X,&x_arr);CHKERRQ(ierr); 434 ierr = VecRestoreArray(F,&f_arr);CHKERRQ(ierr); 435 ierr = VecRestoreArray(Xl,&l);CHKERRQ(ierr); 436 ierr = VecRestoreArray(Xu,&u);CHKERRQ(ierr); 437 ierr = VecRestoreArray(phi,&phi_arr);CHKERRQ(ierr); 438 PetscFunctionReturn(0); 439 } 440 441 /* 442 SNESVIComputeBsubdifferentialVectors - Computes the diagonal shift (Da) and row scaling (Db) vectors needed for the 443 the semismooth jacobian. 444 */ 445 #undef __FUNCT__ 446 #define __FUNCT__ "SNESVIComputeBsubdifferentialVectors" 447 PetscErrorCode SNESVIComputeBsubdifferentialVectors(SNES snes,Vec X,Vec F,Mat jac,Vec Da,Vec Db) 448 { 449 PetscErrorCode ierr; 450 SNES_VI *vi = (SNES_VI*)snes->data; 451 PetscScalar *l,*u,*x,*f,*da,*db,da1,da2,db1,db2; 452 PetscInt i,nlocal; 453 454 PetscFunctionBegin; 455 456 ierr = VecGetArray(X,&x);CHKERRQ(ierr); 457 ierr = VecGetArray(F,&f);CHKERRQ(ierr); 458 ierr = VecGetArray(vi->xl,&l);CHKERRQ(ierr); 459 ierr = VecGetArray(vi->xu,&u);CHKERRQ(ierr); 460 ierr = VecGetArray(Da,&da);CHKERRQ(ierr); 461 ierr = VecGetArray(Db,&db);CHKERRQ(ierr); 462 ierr = VecGetLocalSize(X,&nlocal);CHKERRQ(ierr); 463 464 for (i=0;i< nlocal;i++) { 465 if ((PetscRealPart(l[i]) <= SNES_VI_NINF) && (PetscRealPart(u[i]) >= SNES_VI_INF)) {/* no constraints on variable */ 466 da[i] = 0; 467 db[i] = 1; 468 } else if (PetscRealPart(l[i]) <= SNES_VI_NINF) { /* upper bound on variable only */ 469 da[i] = DPhi(u[i] - x[i], -f[i]); 470 db[i] = DPhi(-f[i],u[i] - x[i]); 471 } else if (PetscRealPart(u[i]) >= SNES_VI_INF) { /* lower bound on variable only */ 472 da[i] = DPhi(x[i] - l[i], f[i]); 473 db[i] = DPhi(f[i],x[i] - l[i]); 474 } else if (l[i] == u[i]) { /* fixed variable */ 475 da[i] = 1; 476 db[i] = 0; 477 } else { /* upper and lower bounds on variable */ 478 da1 = DPhi(x[i] - l[i], -Phi(u[i] - x[i], -f[i])); 479 db1 = DPhi(-Phi(u[i] - x[i], -f[i]),x[i] - l[i]); 480 da2 = DPhi(u[i] - x[i], -f[i]); 481 db2 = DPhi(-f[i],u[i] - x[i]); 482 da[i] = da1 + db1*da2; 483 db[i] = db1*db2; 484 } 485 } 486 487 ierr = VecRestoreArray(X,&x);CHKERRQ(ierr); 488 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr); 489 ierr = VecRestoreArray(vi->xl,&l);CHKERRQ(ierr); 490 ierr = VecRestoreArray(vi->xu,&u);CHKERRQ(ierr); 491 ierr = VecRestoreArray(Da,&da);CHKERRQ(ierr); 492 ierr = VecRestoreArray(Db,&db);CHKERRQ(ierr); 493 PetscFunctionReturn(0); 494 } 495 496 /* 497 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. 498 499 Input Parameters: 500 . Da - Diagonal shift vector for the semismooth jacobian. 501 . Db - Row scaling vector for the semismooth jacobian. 502 503 Output Parameters: 504 . jac - semismooth jacobian 505 . jac_pre - optional preconditioning matrix 506 507 Notes: 508 The semismooth jacobian matrix is given by 509 jac = Da + Db*jacfun 510 where Db is the row scaling matrix stored as a vector, 511 Da is the diagonal perturbation matrix stored as a vector 512 and jacfun is the jacobian of the original nonlinear function. 513 */ 514 #undef __FUNCT__ 515 #define __FUNCT__ "SNESVIComputeJacobian" 516 PetscErrorCode SNESVIComputeJacobian(Mat jac, Mat jac_pre,Vec Da, Vec Db) 517 { 518 PetscErrorCode ierr; 519 520 /* Do row scaling and add diagonal perturbation */ 521 ierr = MatDiagonalScale(jac,Db,PETSC_NULL);CHKERRQ(ierr); 522 ierr = MatDiagonalSet(jac,Da,ADD_VALUES);CHKERRQ(ierr); 523 if (jac != jac_pre) { /* If jac and jac_pre are different */ 524 ierr = MatDiagonalScale(jac_pre,Db,PETSC_NULL); 525 ierr = MatDiagonalSet(jac_pre,Da,ADD_VALUES);CHKERRQ(ierr); 526 } 527 PetscFunctionReturn(0); 528 } 529 530 /* 531 SNESVIComputeMeritFunctionGradient - Computes the gradient of the merit function psi. 532 533 Input Parameters: 534 phi - semismooth function. 535 H - semismooth jacobian 536 537 Output Parameters: 538 dpsi - merit function gradient 539 540 Notes: 541 The merit function gradient is computed as follows 542 dpsi = H^T*phi 543 */ 544 #undef __FUNCT__ 545 #define __FUNCT__ "SNESVIComputeMeritFunctionGradient" 546 PetscErrorCode SNESVIComputeMeritFunctionGradient(Mat H, Vec phi, Vec dpsi) 547 { 548 PetscErrorCode ierr; 549 550 PetscFunctionBegin; 551 ierr = MatMultTranspose(H,phi,dpsi);CHKERRQ(ierr); 552 PetscFunctionReturn(0); 553 } 554 555 /* -------------------------------------------------------------------------- */ 556 /* 557 SNESVIProjectOntoBounds - Projects X onto the feasible region so that Xl[i] <= X[i] <= Xu[i] for i = 1...n. 558 559 Input Parameters: 560 . SNES - nonlinear solver context 561 562 Output Parameters: 563 . X - Bound projected X 564 565 */ 566 567 #undef __FUNCT__ 568 #define __FUNCT__ "SNESVIProjectOntoBounds" 569 PetscErrorCode SNESVIProjectOntoBounds(SNES snes,Vec X) 570 { 571 PetscErrorCode ierr; 572 SNES_VI *vi = (SNES_VI*)snes->data; 573 const PetscScalar *xl,*xu; 574 PetscScalar *x; 575 PetscInt i,n; 576 577 PetscFunctionBegin; 578 ierr = VecGetLocalSize(X,&n);CHKERRQ(ierr); 579 ierr = VecGetArray(X,&x);CHKERRQ(ierr); 580 ierr = VecGetArrayRead(vi->xl,&xl);CHKERRQ(ierr); 581 ierr = VecGetArrayRead(vi->xu,&xu);CHKERRQ(ierr); 582 583 for(i = 0;i<n;i++) { 584 if (PetscRealPart(x[i]) < PetscRealPart(xl[i])) x[i] = xl[i]; 585 else if (PetscRealPart(x[i]) > PetscRealPart(xu[i])) x[i] = xu[i]; 586 } 587 ierr = VecRestoreArray(X,&x);CHKERRQ(ierr); 588 ierr = VecRestoreArrayRead(vi->xl,&xl);CHKERRQ(ierr); 589 ierr = VecRestoreArrayRead(vi->xu,&xu);CHKERRQ(ierr); 590 PetscFunctionReturn(0); 591 } 592 593 /* -------------------------------------------------------------------- 594 595 This file implements a semismooth truncated Newton method with a line search, 596 for solving a system of nonlinear equations in complementarity form, using the KSP, Vec, 597 and Mat interfaces for linear solvers, vectors, and matrices, 598 respectively. 599 600 The following basic routines are required for each nonlinear solver: 601 SNESCreate_XXX() - Creates a nonlinear solver context 602 SNESSetFromOptions_XXX() - Sets runtime options 603 SNESSolve_XXX() - Solves the nonlinear system 604 SNESDestroy_XXX() - Destroys the nonlinear solver context 605 The suffix "_XXX" denotes a particular implementation, in this case 606 we use _VI (e.g., SNESCreate_VI, SNESSolve_VI) for solving 607 systems of nonlinear equations with a line search (LS) method. 608 These routines are actually called via the common user interface 609 routines SNESCreate(), SNESSetFromOptions(), SNESSolve(), and 610 SNESDestroy(), so the application code interface remains identical 611 for all nonlinear solvers. 612 613 Another key routine is: 614 SNESSetUp_XXX() - Prepares for the use of a nonlinear solver 615 by setting data structures and options. The interface routine SNESSetUp() 616 is not usually called directly by the user, but instead is called by 617 SNESSolve() if necessary. 618 619 Additional basic routines are: 620 SNESView_XXX() - Prints details of runtime options that 621 have actually been used. 622 These are called by application codes via the interface routines 623 SNESView(). 624 625 The various types of solvers (preconditioners, Krylov subspace methods, 626 nonlinear solvers, timesteppers) are all organized similarly, so the 627 above description applies to these categories also. 628 629 -------------------------------------------------------------------- */ 630 /* 631 SNESSolveVI_SS - Solves the complementarity problem with a semismooth Newton 632 method using a line search. 633 634 Input Parameters: 635 . snes - the SNES context 636 637 Output Parameter: 638 . outits - number of iterations until termination 639 640 Application Interface Routine: SNESSolve() 641 642 Notes: 643 This implements essentially a semismooth Newton method with a 644 line search. The default line search does not do any line seach 645 but rather takes a full newton step. 646 */ 647 #undef __FUNCT__ 648 #define __FUNCT__ "SNESSolveVI_SS" 649 PetscErrorCode SNESSolveVI_SS(SNES snes) 650 { 651 SNES_VI *vi = (SNES_VI*)snes->data; 652 PetscErrorCode ierr; 653 PetscInt maxits,i,lits; 654 PetscBool lssucceed; 655 MatStructure flg = DIFFERENT_NONZERO_PATTERN; 656 PetscReal gnorm,xnorm=0,ynorm; 657 Vec Y,X,F,G,W; 658 KSPConvergedReason kspreason; 659 660 PetscFunctionBegin; 661 vi->computeuserfunction = snes->ops->computefunction; 662 snes->ops->computefunction = SNESVIComputeFunction; 663 664 snes->numFailures = 0; 665 snes->numLinearSolveFailures = 0; 666 snes->reason = SNES_CONVERGED_ITERATING; 667 668 maxits = snes->max_its; /* maximum number of iterations */ 669 X = snes->vec_sol; /* solution vector */ 670 F = snes->vec_func; /* residual vector */ 671 Y = snes->work[0]; /* work vectors */ 672 G = snes->work[1]; 673 W = snes->work[2]; 674 675 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 676 snes->iter = 0; 677 snes->norm = 0.0; 678 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 679 680 ierr = SNESVIProjectOntoBounds(snes,X);CHKERRQ(ierr); 681 ierr = SNESComputeFunction(snes,X,vi->phi);CHKERRQ(ierr); 682 if (snes->domainerror) { 683 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 684 snes->ops->computefunction = vi->computeuserfunction; 685 PetscFunctionReturn(0); 686 } 687 /* Compute Merit function */ 688 ierr = SNESVIComputeMeritFunction(vi->phi,&vi->merit,&vi->phinorm);CHKERRQ(ierr); 689 690 ierr = VecNormBegin(X,NORM_2,&xnorm);CHKERRQ(ierr); /* xnorm <- ||x|| */ 691 ierr = VecNormEnd(X,NORM_2,&xnorm);CHKERRQ(ierr); 692 if (PetscIsInfOrNanReal(vi->merit)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 693 694 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 695 snes->norm = vi->phinorm; 696 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 697 SNESLogConvHistory(snes,vi->phinorm,0); 698 ierr = SNESMonitor(snes,0,vi->phinorm);CHKERRQ(ierr); 699 700 /* set parameter for default relative tolerance convergence test */ 701 snes->ttol = vi->phinorm*snes->rtol; 702 /* test convergence */ 703 ierr = (*snes->ops->converged)(snes,0,0.0,0.0,vi->phinorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 704 if (snes->reason) { 705 snes->ops->computefunction = vi->computeuserfunction; 706 PetscFunctionReturn(0); 707 } 708 709 for (i=0; i<maxits; i++) { 710 711 /* Call general purpose update function */ 712 if (snes->ops->update) { 713 ierr = (*snes->ops->update)(snes, snes->iter);CHKERRQ(ierr); 714 } 715 716 /* Solve J Y = Phi, where J is the semismooth jacobian */ 717 /* Get the nonlinear function jacobian */ 718 ierr = SNESComputeJacobian(snes,X,&snes->jacobian,&snes->jacobian_pre,&flg);CHKERRQ(ierr); 719 /* Get the diagonal shift and row scaling vectors */ 720 ierr = SNESVIComputeBsubdifferentialVectors(snes,X,F,snes->jacobian,vi->Da,vi->Db);CHKERRQ(ierr); 721 /* Compute the semismooth jacobian */ 722 ierr = SNESVIComputeJacobian(snes->jacobian,snes->jacobian_pre,vi->Da,vi->Db);CHKERRQ(ierr); 723 /* Compute the merit function gradient */ 724 ierr = SNESVIComputeMeritFunctionGradient(snes->jacobian,vi->phi,vi->dpsi);CHKERRQ(ierr); 725 ierr = KSPSetOperators(snes->ksp,snes->jacobian,snes->jacobian_pre,flg);CHKERRQ(ierr); 726 ierr = SNES_KSPSolve(snes,snes->ksp,vi->phi,Y);CHKERRQ(ierr); 727 ierr = KSPGetConvergedReason(snes->ksp,&kspreason);CHKERRQ(ierr); 728 729 if (kspreason < 0) { 730 if (++snes->numLinearSolveFailures >= snes->maxLinearSolveFailures) { 731 ierr = PetscInfo2(snes,"iter=%D, number linear solve failures %D greater than current SNES allowed, stopping solve\n",snes->iter,snes->numLinearSolveFailures);CHKERRQ(ierr); 732 snes->reason = SNES_DIVERGED_LINEAR_SOLVE; 733 break; 734 } 735 } 736 ierr = KSPGetIterationNumber(snes->ksp,&lits);CHKERRQ(ierr); 737 snes->linear_its += lits; 738 ierr = PetscInfo2(snes,"iter=%D, linear solve iterations=%D\n",snes->iter,lits);CHKERRQ(ierr); 739 /* 740 if (vi->precheckstep) { 741 PetscBool changed_y = PETSC_FALSE; 742 ierr = (*vi->precheckstep)(snes,X,Y,vi->precheck,&changed_y);CHKERRQ(ierr); 743 } 744 745 if (PetscLogPrintInfo){ 746 ierr = SNESVICheckResidual_Private(snes,snes->jacobian,F,Y,G,W);CHKERRQ(ierr); 747 } 748 */ 749 /* Compute a (scaled) negative update in the line search routine: 750 Y <- X - lambda*Y 751 and evaluate G = function(Y) (depends on the line search). 752 */ 753 ierr = VecCopy(Y,snes->vec_sol_update);CHKERRQ(ierr); 754 ynorm = 1; gnorm = vi->phinorm; 755 ierr = (*vi->LineSearch)(snes,vi->lsP,X,vi->phi,G,Y,W,vi->phinorm,xnorm,&ynorm,&gnorm,&lssucceed);CHKERRQ(ierr); 756 ierr = PetscInfo4(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lssucceed=%d\n",vi->phinorm,gnorm,ynorm,(int)lssucceed);CHKERRQ(ierr); 757 if (snes->reason == SNES_DIVERGED_FUNCTION_COUNT) break; 758 if (snes->domainerror) { 759 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 760 snes->ops->computefunction = vi->computeuserfunction; 761 PetscFunctionReturn(0); 762 } 763 if (!lssucceed) { 764 if (++snes->numFailures >= snes->maxFailures) { 765 PetscBool ismin; 766 snes->reason = SNES_DIVERGED_LINE_SEARCH; 767 ierr = SNESVICheckLocalMin_Private(snes,snes->jacobian,G,W,gnorm,&ismin);CHKERRQ(ierr); 768 if (ismin) snes->reason = SNES_DIVERGED_LOCAL_MIN; 769 break; 770 } 771 } 772 /* Update function and solution vectors */ 773 vi->phinorm = gnorm; 774 vi->merit = 0.5*vi->phinorm*vi->phinorm; 775 ierr = VecCopy(G,vi->phi);CHKERRQ(ierr); 776 ierr = VecCopy(W,X);CHKERRQ(ierr); 777 /* Monitor convergence */ 778 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 779 snes->iter = i+1; 780 snes->norm = vi->phinorm; 781 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 782 SNESLogConvHistory(snes,snes->norm,lits); 783 ierr = SNESMonitor(snes,snes->iter,snes->norm);CHKERRQ(ierr); 784 /* Test for convergence, xnorm = || X || */ 785 if (snes->ops->converged != SNESSkipConverged) { ierr = VecNorm(X,NORM_2,&xnorm);CHKERRQ(ierr); } 786 ierr = (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,vi->phinorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 787 if (snes->reason) break; 788 } 789 if (i == maxits) { 790 ierr = PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",maxits);CHKERRQ(ierr); 791 if(!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT; 792 } 793 snes->ops->computefunction = vi->computeuserfunction; 794 PetscFunctionReturn(0); 795 } 796 797 #undef __FUNCT__ 798 #define __FUNCT__ "SNESVIGetActiveSetIS" 799 /* 800 SNESVIGetActiveSetIndices - Gets the global indices for the active set variables 801 802 Input parameter 803 . snes - the SNES context 804 . X - the snes solution vector 805 . F - the nonlinear function vector 806 807 Output parameter 808 . ISact - active set index set 809 */ 810 PetscErrorCode SNESVIGetActiveSetIS(SNES snes,Vec X,Vec F,IS* ISact) 811 { 812 PetscErrorCode ierr; 813 SNES_VI *vi = (SNES_VI*)snes->data; 814 Vec Xl=vi->xl,Xu=vi->xu; 815 const PetscScalar *x,*f,*xl,*xu; 816 PetscInt *idx_act,i,nlocal,nloc_isact=0,ilow,ihigh,i1=0; 817 818 PetscFunctionBegin; 819 ierr = VecGetLocalSize(X,&nlocal);CHKERRQ(ierr); 820 ierr = VecGetOwnershipRange(X,&ilow,&ihigh);CHKERRQ(ierr); 821 ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr); 822 ierr = VecGetArrayRead(Xl,&xl);CHKERRQ(ierr); 823 ierr = VecGetArrayRead(Xu,&xu);CHKERRQ(ierr); 824 ierr = VecGetArrayRead(F,&f);CHKERRQ(ierr); 825 /* Compute active set size */ 826 for (i=0; i < nlocal;i++) { 827 if (!((PetscRealPart(x[i]) > PetscRealPart(xl[i]) + 1.e-8 || (PetscRealPart(f[i]) < 0.0)) && ((PetscRealPart(x[i]) < PetscRealPart(xu[i]) - 1.e-8) || PetscRealPart(f[i]) > 0.0))) nloc_isact++; 828 } 829 830 ierr = PetscMalloc(nloc_isact*sizeof(PetscInt),&idx_act);CHKERRQ(ierr); 831 832 /* Set active set indices */ 833 for(i=0; i < nlocal; i++) { 834 if (!((PetscRealPart(x[i]) > PetscRealPart(xl[i]) + 1.e-8 || (PetscRealPart(f[i]) < 0.0)) && ((PetscRealPart(x[i]) < PetscRealPart(xu[i]) - 1.e-8) || PetscRealPart(f[i]) > 0.0))) idx_act[i1++] = ilow+i; 835 } 836 837 /* Create active set IS */ 838 ierr = ISCreateGeneral(((PetscObject)snes)->comm,nloc_isact,idx_act,PETSC_OWN_POINTER,ISact);CHKERRQ(ierr); 839 840 ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr); 841 ierr = VecRestoreArrayRead(Xl,&xl);CHKERRQ(ierr); 842 ierr = VecRestoreArrayRead(Xu,&xu);CHKERRQ(ierr); 843 ierr = VecRestoreArrayRead(F,&f);CHKERRQ(ierr); 844 PetscFunctionReturn(0); 845 } 846 847 #undef __FUNCT__ 848 #define __FUNCT__ "SNESVICreateIndexSets_RS" 849 PetscErrorCode SNESVICreateIndexSets_RS(SNES snes,Vec X,Vec F,IS* ISact,IS* ISinact) 850 { 851 PetscErrorCode ierr; 852 853 PetscFunctionBegin; 854 ierr = SNESVIGetActiveSetIS(snes,X,F,ISact);CHKERRQ(ierr); 855 ierr = ISComplement(*ISact,X->map->rstart,X->map->rend,ISinact);CHKERRQ(ierr); 856 PetscFunctionReturn(0); 857 } 858 859 /* Create active and inactive set vectors. The local size of this vector is set and petsc computes the global size */ 860 #undef __FUNCT__ 861 #define __FUNCT__ "SNESVICreateSubVectors" 862 PetscErrorCode SNESVICreateSubVectors(SNES snes,PetscInt n,Vec* newv) 863 { 864 PetscErrorCode ierr; 865 Vec v; 866 867 PetscFunctionBegin; 868 ierr = VecCreate(((PetscObject)snes)->comm,&v);CHKERRQ(ierr); 869 ierr = VecSetSizes(v,n,PETSC_DECIDE);CHKERRQ(ierr); 870 ierr = VecSetFromOptions(v);CHKERRQ(ierr); 871 *newv = v; 872 873 PetscFunctionReturn(0); 874 } 875 876 /* Resets the snes PC and KSP when the active set sizes change */ 877 #undef __FUNCT__ 878 #define __FUNCT__ "SNESVIResetPCandKSP" 879 PetscErrorCode SNESVIResetPCandKSP(SNES snes,Mat Amat,Mat Pmat) 880 { 881 PetscErrorCode ierr; 882 KSP snesksp; 883 884 PetscFunctionBegin; 885 ierr = SNESGetKSP(snes,&snesksp);CHKERRQ(ierr); 886 ierr = KSPReset(snesksp);CHKERRQ(ierr); 887 888 /* 889 KSP kspnew; 890 PC pcnew; 891 const MatSolverPackage stype; 892 893 894 ierr = KSPCreate(((PetscObject)snes)->comm,&kspnew);CHKERRQ(ierr); 895 kspnew->pc_side = snesksp->pc_side; 896 kspnew->rtol = snesksp->rtol; 897 kspnew->abstol = snesksp->abstol; 898 kspnew->max_it = snesksp->max_it; 899 ierr = KSPSetType(kspnew,((PetscObject)snesksp)->type_name);CHKERRQ(ierr); 900 ierr = KSPGetPC(kspnew,&pcnew);CHKERRQ(ierr); 901 ierr = PCSetType(kspnew->pc,((PetscObject)snesksp->pc)->type_name);CHKERRQ(ierr); 902 ierr = PCSetOperators(kspnew->pc,Amat,Pmat,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); 903 ierr = PCFactorGetMatSolverPackage(snesksp->pc,&stype);CHKERRQ(ierr); 904 ierr = PCFactorSetMatSolverPackage(kspnew->pc,stype);CHKERRQ(ierr); 905 ierr = KSPDestroy(&snesksp);CHKERRQ(ierr); 906 snes->ksp = kspnew; 907 ierr = PetscLogObjectParent(snes,kspnew);CHKERRQ(ierr); 908 ierr = KSPSetFromOptions(kspnew);CHKERRQ(ierr);*/ 909 PetscFunctionReturn(0); 910 } 911 912 913 #undef __FUNCT__ 914 #define __FUNCT__ "SNESVIComputeInactiveSetFnorm" 915 PetscErrorCode SNESVIComputeInactiveSetFnorm(SNES snes,Vec F,Vec X,PetscReal *fnorm) 916 { 917 PetscErrorCode ierr; 918 SNES_VI *vi = (SNES_VI*)snes->data; 919 const PetscScalar *x,*xl,*xu,*f; 920 PetscInt i,n; 921 PetscReal rnorm; 922 923 PetscFunctionBegin; 924 ierr = VecGetLocalSize(X,&n);CHKERRQ(ierr); 925 ierr = VecGetArrayRead(vi->xl,&xl);CHKERRQ(ierr); 926 ierr = VecGetArrayRead(vi->xu,&xu);CHKERRQ(ierr); 927 ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr); 928 ierr = VecGetArrayRead(F,&f);CHKERRQ(ierr); 929 rnorm = 0.0; 930 for (i=0; i<n; i++) { 931 if (((PetscRealPart(x[i]) > PetscRealPart(xl[i]) + 1.e-8 || (PetscRealPart(f[i]) < 0.0)) && ((PetscRealPart(x[i]) < PetscRealPart(xu[i]) - 1.e-8) || PetscRealPart(f[i]) > 0.0))) rnorm += PetscRealPart(PetscConj(f[i])*f[i]); 932 } 933 ierr = VecRestoreArrayRead(F,&f);CHKERRQ(ierr); 934 ierr = VecRestoreArrayRead(vi->xl,&xl);CHKERRQ(ierr); 935 ierr = VecRestoreArrayRead(vi->xu,&xu);CHKERRQ(ierr); 936 ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr); 937 ierr = MPI_Allreduce(&rnorm,fnorm,1,MPIU_REAL,MPIU_SUM,((PetscObject)snes)->comm);CHKERRQ(ierr); 938 *fnorm = sqrt(*fnorm); 939 PetscFunctionReturn(0); 940 } 941 942 /* Variational Inequality solver using reduce space method. No semismooth algorithm is 943 implemented in this algorithm. It basically identifies the active constraints and does 944 a linear solve on the other variables (those not associated with the active constraints). */ 945 #undef __FUNCT__ 946 #define __FUNCT__ "SNESSolveVI_RS" 947 PetscErrorCode SNESSolveVI_RS(SNES snes) 948 { 949 SNES_VI *vi = (SNES_VI*)snes->data; 950 PetscErrorCode ierr; 951 PetscInt maxits,i,lits; 952 PetscBool lssucceed; 953 MatStructure flg = DIFFERENT_NONZERO_PATTERN; 954 PetscReal fnorm,gnorm,xnorm=0,ynorm; 955 Vec Y,X,F,G,W; 956 KSPConvergedReason kspreason; 957 958 PetscFunctionBegin; 959 snes->numFailures = 0; 960 snes->numLinearSolveFailures = 0; 961 snes->reason = SNES_CONVERGED_ITERATING; 962 963 maxits = snes->max_its; /* maximum number of iterations */ 964 X = snes->vec_sol; /* solution vector */ 965 F = snes->vec_func; /* residual vector */ 966 Y = snes->work[0]; /* work vectors */ 967 G = snes->work[1]; 968 W = snes->work[2]; 969 970 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 971 snes->iter = 0; 972 snes->norm = 0.0; 973 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 974 975 ierr = SNESVIProjectOntoBounds(snes,X);CHKERRQ(ierr); 976 ierr = SNESComputeFunction(snes,X,F);CHKERRQ(ierr); 977 if (snes->domainerror) { 978 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 979 PetscFunctionReturn(0); 980 } 981 ierr = SNESVIComputeInactiveSetFnorm(snes,F,X,&fnorm);CHKERRQ(ierr); 982 ierr = VecNormBegin(X,NORM_2,&xnorm);CHKERRQ(ierr); /* xnorm <- ||x|| */ 983 ierr = VecNormEnd(X,NORM_2,&xnorm);CHKERRQ(ierr); 984 if (PetscIsInfOrNanReal(fnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 985 986 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 987 snes->norm = fnorm; 988 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 989 SNESLogConvHistory(snes,fnorm,0); 990 ierr = SNESMonitor(snes,0,fnorm);CHKERRQ(ierr); 991 992 /* set parameter for default relative tolerance convergence test */ 993 snes->ttol = fnorm*snes->rtol; 994 /* test convergence */ 995 ierr = (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 996 if (snes->reason) PetscFunctionReturn(0); 997 998 for (i=0; i<maxits; i++) { 999 1000 IS IS_act,IS_inact; /* _act -> active set _inact -> inactive set */ 1001 VecScatter scat_act,scat_inact; 1002 PetscInt nis_act,nis_inact; 1003 Vec Y_act,Y_inact,F_inact; 1004 Mat jac_inact_inact,prejac_inact_inact; 1005 IS keptrows; 1006 PetscBool isequal; 1007 1008 /* Call general purpose update function */ 1009 if (snes->ops->update) { 1010 ierr = (*snes->ops->update)(snes, snes->iter);CHKERRQ(ierr); 1011 } 1012 ierr = SNESComputeJacobian(snes,X,&snes->jacobian,&snes->jacobian_pre,&flg);CHKERRQ(ierr); 1013 1014 /* Create active and inactive index sets */ 1015 ierr = SNESVICreateIndexSets_RS(snes,X,F,&IS_act,&IS_inact);CHKERRQ(ierr); 1016 1017 /* Create inactive set submatrix */ 1018 ierr = MatGetSubMatrix(snes->jacobian,IS_inact,IS_inact,MAT_INITIAL_MATRIX,&jac_inact_inact);CHKERRQ(ierr); 1019 ierr = MatFindNonzeroRows(jac_inact_inact,&keptrows);CHKERRQ(ierr); 1020 if (keptrows) { 1021 PetscInt cnt,*nrows,k; 1022 const PetscInt *krows,*inact; 1023 PetscInt rstart=jac_inact_inact->rmap->rstart; 1024 1025 ierr = MatDestroy(&jac_inact_inact);CHKERRQ(ierr); 1026 ierr = ISDestroy(&IS_act);CHKERRQ(ierr); 1027 1028 ierr = ISGetLocalSize(keptrows,&cnt);CHKERRQ(ierr); 1029 ierr = ISGetIndices(keptrows,&krows);CHKERRQ(ierr); 1030 ierr = ISGetIndices(IS_inact,&inact);CHKERRQ(ierr); 1031 ierr = PetscMalloc(cnt*sizeof(PetscInt),&nrows);CHKERRQ(ierr); 1032 for (k=0; k<cnt; k++) { 1033 nrows[k] = inact[krows[k]-rstart]; 1034 } 1035 ierr = ISRestoreIndices(keptrows,&krows);CHKERRQ(ierr); 1036 ierr = ISRestoreIndices(IS_inact,&inact);CHKERRQ(ierr); 1037 ierr = ISDestroy(&keptrows);CHKERRQ(ierr); 1038 ierr = ISDestroy(&IS_inact);CHKERRQ(ierr); 1039 1040 ierr = ISCreateGeneral(PETSC_COMM_WORLD,cnt,nrows,PETSC_OWN_POINTER,&IS_inact);CHKERRQ(ierr); 1041 ierr = ISComplement(IS_inact,F->map->rstart,F->map->rend,&IS_act);CHKERRQ(ierr); 1042 ierr = MatGetSubMatrix(snes->jacobian,IS_inact,IS_inact,MAT_INITIAL_MATRIX,&jac_inact_inact);CHKERRQ(ierr); 1043 } 1044 1045 /* Get sizes of active and inactive sets */ 1046 ierr = ISGetLocalSize(IS_act,&nis_act);CHKERRQ(ierr); 1047 ierr = ISGetLocalSize(IS_inact,&nis_inact);CHKERRQ(ierr); 1048 1049 /* Create active and inactive set vectors */ 1050 ierr = SNESVICreateSubVectors(snes,nis_inact,&F_inact);CHKERRQ(ierr); 1051 ierr = SNESVICreateSubVectors(snes,nis_act,&Y_act);CHKERRQ(ierr); 1052 ierr = SNESVICreateSubVectors(snes,nis_inact,&Y_inact);CHKERRQ(ierr); 1053 1054 /* Create scatter contexts */ 1055 ierr = VecScatterCreate(Y,IS_act,Y_act,PETSC_NULL,&scat_act);CHKERRQ(ierr); 1056 ierr = VecScatterCreate(Y,IS_inact,Y_inact,PETSC_NULL,&scat_inact);CHKERRQ(ierr); 1057 1058 /* Do a vec scatter to active and inactive set vectors */ 1059 ierr = VecScatterBegin(scat_inact,F,F_inact,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1060 ierr = VecScatterEnd(scat_inact,F,F_inact,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1061 1062 ierr = VecScatterBegin(scat_act,Y,Y_act,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1063 ierr = VecScatterEnd(scat_act,Y,Y_act,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1064 1065 ierr = VecScatterBegin(scat_inact,Y,Y_inact,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1066 ierr = VecScatterEnd(scat_inact,Y,Y_inact,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1067 1068 /* Active set direction = 0 */ 1069 ierr = VecSet(Y_act,0);CHKERRQ(ierr); 1070 if (snes->jacobian != snes->jacobian_pre) { 1071 ierr = MatGetSubMatrix(snes->jacobian_pre,IS_inact,IS_inact,MAT_INITIAL_MATRIX,&prejac_inact_inact);CHKERRQ(ierr); 1072 } else prejac_inact_inact = jac_inact_inact; 1073 1074 ierr = ISEqual(vi->IS_inact_prev,IS_inact,&isequal);CHKERRQ(ierr); 1075 if (!isequal) { 1076 ierr = SNESVIResetPCandKSP(snes,jac_inact_inact,prejac_inact_inact);CHKERRQ(ierr); 1077 } 1078 1079 /* ierr = ISView(IS_inact,0);CHKERRQ(ierr); */ 1080 /* ierr = ISView(IS_act,0);CHKERRQ(ierr);*/ 1081 /* ierr = MatView(snes->jacobian_pre,0); */ 1082 1083 ierr = KSPSetOperators(snes->ksp,jac_inact_inact,prejac_inact_inact,flg);CHKERRQ(ierr); 1084 ierr = KSPSetUp(snes->ksp);CHKERRQ(ierr); 1085 { 1086 PC pc; 1087 PetscBool flg; 1088 ierr = KSPGetPC(snes->ksp,&pc);CHKERRQ(ierr); 1089 ierr = PetscTypeCompare((PetscObject)pc,PCFIELDSPLIT,&flg);CHKERRQ(ierr); 1090 if (flg) { 1091 KSP *subksps; 1092 ierr = PCFieldSplitGetSubKSP(pc,PETSC_NULL,&subksps);CHKERRQ(ierr); 1093 ierr = KSPGetPC(subksps[0],&pc);CHKERRQ(ierr); 1094 ierr = PetscFree(subksps);CHKERRQ(ierr); 1095 ierr = PetscTypeCompare((PetscObject)pc,PCBJACOBI,&flg);CHKERRQ(ierr); 1096 if (flg) { 1097 PetscInt n,N = 101*101,j,cnts[3] = {0,0,0}; 1098 const PetscInt *ii; 1099 1100 ierr = ISGetSize(IS_inact,&n);CHKERRQ(ierr); 1101 ierr = ISGetIndices(IS_inact,&ii);CHKERRQ(ierr); 1102 for (j=0; j<n; j++) { 1103 if (ii[j] < N) cnts[0]++; 1104 else if (ii[j] < 2*N) cnts[1]++; 1105 else if (ii[j] < 3*N) cnts[2]++; 1106 } 1107 ierr = ISRestoreIndices(IS_inact,&ii);CHKERRQ(ierr); 1108 1109 ierr = PCBJacobiSetTotalBlocks(pc,3,cnts);CHKERRQ(ierr); 1110 } 1111 } 1112 } 1113 1114 ierr = SNES_KSPSolve(snes,snes->ksp,F_inact,Y_inact);CHKERRQ(ierr); 1115 ierr = KSPGetConvergedReason(snes->ksp,&kspreason);CHKERRQ(ierr); 1116 if (kspreason < 0) { 1117 if (++snes->numLinearSolveFailures >= snes->maxLinearSolveFailures) { 1118 ierr = PetscInfo2(snes,"iter=%D, number linear solve failures %D greater than current SNES allowed, stopping solve\n",snes->iter,snes->numLinearSolveFailures);CHKERRQ(ierr); 1119 snes->reason = SNES_DIVERGED_LINEAR_SOLVE; 1120 break; 1121 } 1122 } 1123 1124 ierr = VecScatterBegin(scat_act,Y_act,Y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1125 ierr = VecScatterEnd(scat_act,Y_act,Y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1126 ierr = VecScatterBegin(scat_inact,Y_inact,Y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1127 ierr = VecScatterEnd(scat_inact,Y_inact,Y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1128 1129 ierr = VecDestroy(&F_inact);CHKERRQ(ierr); 1130 ierr = VecDestroy(&Y_act);CHKERRQ(ierr); 1131 ierr = VecDestroy(&Y_inact);CHKERRQ(ierr); 1132 ierr = VecScatterDestroy(&scat_act);CHKERRQ(ierr); 1133 ierr = VecScatterDestroy(&scat_inact);CHKERRQ(ierr); 1134 ierr = ISDestroy(&IS_act);CHKERRQ(ierr); 1135 if (!isequal) { 1136 ierr = ISDestroy(&vi->IS_inact_prev);CHKERRQ(ierr); 1137 ierr = ISDuplicate(IS_inact,&vi->IS_inact_prev);CHKERRQ(ierr); 1138 } 1139 ierr = ISDestroy(&IS_inact);CHKERRQ(ierr); 1140 ierr = MatDestroy(&jac_inact_inact);CHKERRQ(ierr); 1141 if (snes->jacobian != snes->jacobian_pre) { 1142 ierr = MatDestroy(&prejac_inact_inact);CHKERRQ(ierr); 1143 } 1144 ierr = KSPGetIterationNumber(snes->ksp,&lits);CHKERRQ(ierr); 1145 snes->linear_its += lits; 1146 ierr = PetscInfo2(snes,"iter=%D, linear solve iterations=%D\n",snes->iter,lits);CHKERRQ(ierr); 1147 /* 1148 if (vi->precheckstep) { 1149 PetscBool changed_y = PETSC_FALSE; 1150 ierr = (*vi->precheckstep)(snes,X,Y,vi->precheck,&changed_y);CHKERRQ(ierr); 1151 } 1152 1153 if (PetscLogPrintInfo){ 1154 ierr = SNESVICheckResidual_Private(snes,snes->jacobian,F,Y,G,W);CHKERRQ(ierr); 1155 } 1156 */ 1157 /* Compute a (scaled) negative update in the line search routine: 1158 Y <- X - lambda*Y 1159 and evaluate G = function(Y) (depends on the line search). 1160 */ 1161 ierr = VecCopy(Y,snes->vec_sol_update);CHKERRQ(ierr); 1162 ynorm = 1; gnorm = fnorm; 1163 ierr = (*vi->LineSearch)(snes,vi->lsP,X,F,G,Y,W,fnorm,xnorm,&ynorm,&gnorm,&lssucceed);CHKERRQ(ierr); 1164 ierr = PetscInfo4(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lssucceed=%d\n",fnorm,gnorm,ynorm,(int)lssucceed);CHKERRQ(ierr); 1165 if (snes->reason == SNES_DIVERGED_FUNCTION_COUNT) break; 1166 if (snes->domainerror) { 1167 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 1168 PetscFunctionReturn(0); 1169 } 1170 if (!lssucceed) { 1171 if (++snes->numFailures >= snes->maxFailures) { 1172 PetscBool ismin; 1173 snes->reason = SNES_DIVERGED_LINE_SEARCH; 1174 ierr = SNESVICheckLocalMin_Private(snes,snes->jacobian,G,W,gnorm,&ismin);CHKERRQ(ierr); 1175 if (ismin) snes->reason = SNES_DIVERGED_LOCAL_MIN; 1176 break; 1177 } 1178 } 1179 /* Update function and solution vectors */ 1180 fnorm = gnorm; 1181 ierr = VecCopy(G,F);CHKERRQ(ierr); 1182 ierr = VecCopy(W,X);CHKERRQ(ierr); 1183 /* Monitor convergence */ 1184 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 1185 snes->iter = i+1; 1186 snes->norm = fnorm; 1187 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 1188 SNESLogConvHistory(snes,snes->norm,lits); 1189 ierr = SNESMonitor(snes,snes->iter,snes->norm);CHKERRQ(ierr); 1190 /* Test for convergence, xnorm = || X || */ 1191 if (snes->ops->converged != SNESSkipConverged) { ierr = VecNorm(X,NORM_2,&xnorm);CHKERRQ(ierr); } 1192 ierr = (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 1193 if (snes->reason) break; 1194 } 1195 if (i == maxits) { 1196 ierr = PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",maxits);CHKERRQ(ierr); 1197 if(!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT; 1198 } 1199 PetscFunctionReturn(0); 1200 } 1201 1202 #undef __FUNCT__ 1203 #define __FUNCT__ "SNESVISetRedundancyCheck" 1204 PetscErrorCode SNESVISetRedundancyCheck(SNES snes,PetscErrorCode (*func)(SNES,IS,IS*,void*),void *ctx) 1205 { 1206 SNES_VI *vi = (SNES_VI*)snes->data; 1207 1208 PetscFunctionBegin; 1209 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1210 vi->checkredundancy = func; 1211 vi->ctxP = ctx; 1212 PetscFunctionReturn(0); 1213 } 1214 1215 #if defined(PETSC_HAVE_MATLAB_ENGINE) 1216 #include <engine.h> 1217 #include <mex.h> 1218 typedef struct {char *funcname; mxArray *ctx;} SNESMatlabContext; 1219 1220 #undef __FUNCT__ 1221 #define __FUNCT__ "SNESVIRedundancyCheck_Matlab" 1222 PetscErrorCode SNESVIRedundancyCheck_Matlab(SNES snes,IS is_act,IS* is_redact,void* ctx) 1223 { 1224 PetscErrorCode ierr; 1225 SNESMatlabContext *sctx = (SNESMatlabContext*)ctx; 1226 int nlhs = 1, nrhs = 5; 1227 mxArray *plhs[1], *prhs[5]; 1228 long long int l1 = 0, l2 = 0, ls = 0; 1229 PetscInt *indices=PETSC_NULL; 1230 1231 PetscFunctionBegin; 1232 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1233 PetscValidHeaderSpecific(is_act,IS_CLASSID,2); 1234 PetscValidPointer(is_redact,3); 1235 PetscCheckSameComm(snes,1,is_act,2); 1236 1237 /* Create IS for reduced active set of size 0, its size and indices will 1238 bet set by the Matlab function */ 1239 ierr = ISCreateGeneral(((PetscObject)snes)->comm,0,indices,PETSC_OWN_POINTER,is_redact);CHKERRQ(ierr); 1240 /* call Matlab function in ctx */ 1241 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 1242 ierr = PetscMemcpy(&l1,&is_act,sizeof(is_act));CHKERRQ(ierr); 1243 ierr = PetscMemcpy(&l2,is_redact,sizeof(is_act));CHKERRQ(ierr); 1244 prhs[0] = mxCreateDoubleScalar((double)ls); 1245 prhs[1] = mxCreateDoubleScalar((double)l1); 1246 prhs[2] = mxCreateDoubleScalar((double)l2); 1247 prhs[3] = mxCreateString(sctx->funcname); 1248 prhs[4] = sctx->ctx; 1249 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESVIRedundancyCheckInternal");CHKERRQ(ierr); 1250 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 1251 mxDestroyArray(prhs[0]); 1252 mxDestroyArray(prhs[1]); 1253 mxDestroyArray(prhs[2]); 1254 mxDestroyArray(prhs[3]); 1255 mxDestroyArray(plhs[0]); 1256 PetscFunctionReturn(0); 1257 } 1258 1259 #undef __FUNCT__ 1260 #define __FUNCT__ "SNESVISetRedundancyCheckMatlab" 1261 PetscErrorCode SNESVISetRedundancyCheckMatlab(SNES snes,const char* func,mxArray* ctx) 1262 { 1263 PetscErrorCode ierr; 1264 SNESMatlabContext *sctx; 1265 1266 PetscFunctionBegin; 1267 /* currently sctx is memory bleed */ 1268 ierr = PetscMalloc(sizeof(SNESMatlabContext),&sctx);CHKERRQ(ierr); 1269 ierr = PetscStrallocpy(func,&sctx->funcname);CHKERRQ(ierr); 1270 sctx->ctx = mxDuplicateArray(ctx); 1271 ierr = SNESVISetRedundancyCheck(snes,SNESVIRedundancyCheck_Matlab,sctx);CHKERRQ(ierr); 1272 PetscFunctionReturn(0); 1273 } 1274 1275 #endif 1276 1277 1278 /* Variational Inequality solver using augmented space method. It does the opposite of the 1279 reduced space method i.e. it identifies the active set variables and instead of discarding 1280 them it augments the original system by introducing additional equality 1281 constraint equations for active set variables. The user can optionally provide an IS for 1282 a subset of 'redundant' active set variables which will treated as free variables. 1283 Specific implementation for Allen-Cahn problem 1284 */ 1285 #undef __FUNCT__ 1286 #define __FUNCT__ "SNESSolveVI_RSAUG" 1287 PetscErrorCode SNESSolveVI_RSAUG(SNES snes) 1288 { 1289 SNES_VI *vi = (SNES_VI*)snes->data; 1290 PetscErrorCode ierr; 1291 PetscInt maxits,i,lits; 1292 PetscBool lssucceed; 1293 MatStructure flg = DIFFERENT_NONZERO_PATTERN; 1294 PetscReal fnorm,gnorm,xnorm=0,ynorm; 1295 Vec Y,X,F,G,W; 1296 KSPConvergedReason kspreason; 1297 1298 PetscFunctionBegin; 1299 snes->numFailures = 0; 1300 snes->numLinearSolveFailures = 0; 1301 snes->reason = SNES_CONVERGED_ITERATING; 1302 1303 maxits = snes->max_its; /* maximum number of iterations */ 1304 X = snes->vec_sol; /* solution vector */ 1305 F = snes->vec_func; /* residual vector */ 1306 Y = snes->work[0]; /* work vectors */ 1307 G = snes->work[1]; 1308 W = snes->work[2]; 1309 1310 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 1311 snes->iter = 0; 1312 snes->norm = 0.0; 1313 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 1314 1315 ierr = SNESVIProjectOntoBounds(snes,X);CHKERRQ(ierr); 1316 ierr = SNESComputeFunction(snes,X,F);CHKERRQ(ierr); 1317 if (snes->domainerror) { 1318 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 1319 PetscFunctionReturn(0); 1320 } 1321 ierr = SNESVIComputeInactiveSetFnorm(snes,F,X,&fnorm);CHKERRQ(ierr); 1322 ierr = VecNormBegin(X,NORM_2,&xnorm);CHKERRQ(ierr); /* xnorm <- ||x|| */ 1323 ierr = VecNormEnd(X,NORM_2,&xnorm);CHKERRQ(ierr); 1324 if (PetscIsInfOrNanReal(fnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1325 1326 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 1327 snes->norm = fnorm; 1328 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 1329 SNESLogConvHistory(snes,fnorm,0); 1330 ierr = SNESMonitor(snes,0,fnorm);CHKERRQ(ierr); 1331 1332 /* set parameter for default relative tolerance convergence test */ 1333 snes->ttol = fnorm*snes->rtol; 1334 /* test convergence */ 1335 ierr = (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 1336 if (snes->reason) PetscFunctionReturn(0); 1337 1338 for (i=0; i<maxits; i++) { 1339 IS IS_act,IS_inact; /* _act -> active set _inact -> inactive set */ 1340 IS IS_redact; /* redundant active set */ 1341 Mat J_aug,Jpre_aug; 1342 Vec F_aug,Y_aug; 1343 PetscInt nis_redact,nis_act; 1344 const PetscInt *idx_redact,*idx_act; 1345 PetscInt k; 1346 PetscInt *idx_actkept=PETSC_NULL,nkept=0; /* list of kept active set */ 1347 PetscScalar *f,*f2; 1348 PetscBool isequal; 1349 PetscInt i1=0,j1=0; 1350 1351 /* Call general purpose update function */ 1352 if (snes->ops->update) { 1353 ierr = (*snes->ops->update)(snes, snes->iter);CHKERRQ(ierr); 1354 } 1355 ierr = SNESComputeJacobian(snes,X,&snes->jacobian,&snes->jacobian_pre,&flg);CHKERRQ(ierr); 1356 1357 /* Create active and inactive index sets */ 1358 ierr = SNESVICreateIndexSets_RS(snes,X,F,&IS_act,&IS_inact);CHKERRQ(ierr); 1359 1360 /* Get local active set size */ 1361 ierr = ISGetLocalSize(IS_act,&nis_act);CHKERRQ(ierr); 1362 if (nis_act) { 1363 ierr = ISGetIndices(IS_act,&idx_act);CHKERRQ(ierr); 1364 IS_redact = PETSC_NULL; 1365 if(vi->checkredundancy) { 1366 (*vi->checkredundancy)(snes,IS_act,&IS_redact,vi->ctxP); 1367 } 1368 1369 if(!IS_redact) { 1370 /* User called checkredundancy function but didn't create IS_redact because 1371 there were no redundant active set variables */ 1372 /* Copy over all active set indices to the list */ 1373 ierr = PetscMalloc(nis_act*sizeof(PetscInt),&idx_actkept);CHKERRQ(ierr); 1374 for(k=0;k < nis_act;k++) idx_actkept[k] = idx_act[k]; 1375 nkept = nis_act; 1376 } else { 1377 ierr = ISGetLocalSize(IS_redact,&nis_redact);CHKERRQ(ierr); 1378 ierr = PetscMalloc((nis_act-nis_redact)*sizeof(PetscInt),&idx_actkept);CHKERRQ(ierr); 1379 1380 /* Create reduced active set list */ 1381 ierr = ISGetIndices(IS_act,&idx_act);CHKERRQ(ierr); 1382 ierr = ISGetIndices(IS_redact,&idx_redact);CHKERRQ(ierr); 1383 j1 = 0;nkept = 0; 1384 for(k=0;k<nis_act;k++) { 1385 if(j1 < nis_redact && idx_act[k] == idx_redact[j1]) j1++; 1386 else idx_actkept[nkept++] = idx_act[k]; 1387 } 1388 ierr = ISRestoreIndices(IS_act,&idx_act);CHKERRQ(ierr); 1389 ierr = ISRestoreIndices(IS_redact,&idx_redact);CHKERRQ(ierr); 1390 1391 ierr = ISDestroy(&IS_redact);CHKERRQ(ierr); 1392 } 1393 1394 /* Create augmented F and Y */ 1395 ierr = VecCreate(((PetscObject)snes)->comm,&F_aug);CHKERRQ(ierr); 1396 ierr = VecSetSizes(F_aug,F->map->n+nkept,PETSC_DECIDE);CHKERRQ(ierr); 1397 ierr = VecSetFromOptions(F_aug);CHKERRQ(ierr); 1398 ierr = VecDuplicate(F_aug,&Y_aug);CHKERRQ(ierr); 1399 1400 /* Copy over F to F_aug in the first n locations */ 1401 ierr = VecGetArray(F,&f);CHKERRQ(ierr); 1402 ierr = VecGetArray(F_aug,&f2);CHKERRQ(ierr); 1403 ierr = PetscMemcpy(f2,f,F->map->n*sizeof(PetscScalar));CHKERRQ(ierr); 1404 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr); 1405 ierr = VecRestoreArray(F_aug,&f2);CHKERRQ(ierr); 1406 1407 /* Create the augmented jacobian matrix */ 1408 ierr = MatCreate(((PetscObject)snes)->comm,&J_aug);CHKERRQ(ierr); 1409 ierr = MatSetSizes(J_aug,X->map->n+nkept,X->map->n+nkept,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr); 1410 ierr = MatSetFromOptions(J_aug);CHKERRQ(ierr); 1411 1412 1413 { /* local vars */ 1414 /* Preallocate augmented matrix and set values in it...Doing only sequential case first*/ 1415 PetscInt ncols; 1416 const PetscInt *cols; 1417 const PetscScalar *vals; 1418 PetscScalar value[2]; 1419 PetscInt row,col[2]; 1420 PetscInt *d_nnz; 1421 value[0] = 1.0; value[1] = 0.0; 1422 ierr = PetscMalloc((X->map->n+nkept)*sizeof(PetscInt),&d_nnz);CHKERRQ(ierr); 1423 ierr = PetscMemzero(d_nnz,(X->map->n+nkept)*sizeof(PetscInt));CHKERRQ(ierr); 1424 for(row=0;row<snes->jacobian->rmap->n;row++) { 1425 ierr = MatGetRow(snes->jacobian,row,&ncols,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 1426 d_nnz[row] += ncols; 1427 ierr = MatRestoreRow(snes->jacobian,row,&ncols,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 1428 } 1429 1430 for(k=0;k<nkept;k++) { 1431 d_nnz[idx_actkept[k]] += 1; 1432 d_nnz[snes->jacobian->rmap->n+k] += 2; 1433 } 1434 ierr = MatSeqAIJSetPreallocation(J_aug,PETSC_NULL,d_nnz);CHKERRQ(ierr); 1435 1436 ierr = PetscFree(d_nnz);CHKERRQ(ierr); 1437 1438 /* Set the original jacobian matrix in J_aug */ 1439 for(row=0;row<snes->jacobian->rmap->n;row++) { 1440 ierr = MatGetRow(snes->jacobian,row,&ncols,&cols,&vals);CHKERRQ(ierr); 1441 ierr = MatSetValues(J_aug,1,&row,ncols,cols,vals,INSERT_VALUES);CHKERRQ(ierr); 1442 ierr = MatRestoreRow(snes->jacobian,row,&ncols,&cols,&vals);CHKERRQ(ierr); 1443 } 1444 /* Add the augmented part */ 1445 for(k=0;k<nkept;k++) { 1446 row = snes->jacobian->rmap->n+k; 1447 col[0] = idx_actkept[k]; col[1] = snes->jacobian->rmap->n+k; 1448 ierr = MatSetValues(J_aug,1,&row,1,col,value,INSERT_VALUES);CHKERRQ(ierr); 1449 ierr = MatSetValues(J_aug,1,&col[0],1,&row,&value[0],INSERT_VALUES);CHKERRQ(ierr); 1450 } 1451 ierr = MatAssemblyBegin(J_aug,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1452 ierr = MatAssemblyEnd(J_aug,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1453 /* Only considering prejac = jac for now */ 1454 Jpre_aug = J_aug; 1455 } /* local vars*/ 1456 ierr = ISRestoreIndices(IS_act,&idx_act);CHKERRQ(ierr); 1457 ierr = ISDestroy(&IS_act);CHKERRQ(ierr); 1458 } else { 1459 F_aug = F; J_aug = snes->jacobian; Y_aug = Y; Jpre_aug = snes->jacobian_pre; 1460 } 1461 1462 ierr = ISEqual(vi->IS_inact_prev,IS_inact,&isequal);CHKERRQ(ierr); 1463 if (!isequal) { 1464 ierr = SNESVIResetPCandKSP(snes,J_aug,Jpre_aug);CHKERRQ(ierr); 1465 } 1466 ierr = KSPSetOperators(snes->ksp,J_aug,Jpre_aug,flg);CHKERRQ(ierr); 1467 ierr = KSPSetUp(snes->ksp);CHKERRQ(ierr); 1468 /* { 1469 PC pc; 1470 PetscBool flg; 1471 ierr = KSPGetPC(snes->ksp,&pc);CHKERRQ(ierr); 1472 ierr = PetscTypeCompare((PetscObject)pc,PCFIELDSPLIT,&flg);CHKERRQ(ierr); 1473 if (flg) { 1474 KSP *subksps; 1475 ierr = PCFieldSplitGetSubKSP(pc,PETSC_NULL,&subksps);CHKERRQ(ierr); 1476 ierr = KSPGetPC(subksps[0],&pc);CHKERRQ(ierr); 1477 ierr = PetscFree(subksps);CHKERRQ(ierr); 1478 ierr = PetscTypeCompare((PetscObject)pc,PCBJACOBI,&flg);CHKERRQ(ierr); 1479 if (flg) { 1480 PetscInt n,N = 101*101,j,cnts[3] = {0,0,0}; 1481 const PetscInt *ii; 1482 1483 ierr = ISGetSize(IS_inact,&n);CHKERRQ(ierr); 1484 ierr = ISGetIndices(IS_inact,&ii);CHKERRQ(ierr); 1485 for (j=0; j<n; j++) { 1486 if (ii[j] < N) cnts[0]++; 1487 else if (ii[j] < 2*N) cnts[1]++; 1488 else if (ii[j] < 3*N) cnts[2]++; 1489 } 1490 ierr = ISRestoreIndices(IS_inact,&ii);CHKERRQ(ierr); 1491 1492 ierr = PCBJacobiSetTotalBlocks(pc,3,cnts);CHKERRQ(ierr); 1493 } 1494 } 1495 } 1496 */ 1497 ierr = SNES_KSPSolve(snes,snes->ksp,F_aug,Y_aug);CHKERRQ(ierr); 1498 ierr = KSPGetConvergedReason(snes->ksp,&kspreason);CHKERRQ(ierr); 1499 if (kspreason < 0) { 1500 if (++snes->numLinearSolveFailures >= snes->maxLinearSolveFailures) { 1501 ierr = PetscInfo2(snes,"iter=%D, number linear solve failures %D greater than current SNES allowed, stopping solve\n",snes->iter,snes->numLinearSolveFailures);CHKERRQ(ierr); 1502 snes->reason = SNES_DIVERGED_LINEAR_SOLVE; 1503 break; 1504 } 1505 } 1506 1507 if(nis_act) { 1508 PetscScalar *y1,*y2; 1509 ierr = VecGetArray(Y,&y1);CHKERRQ(ierr); 1510 ierr = VecGetArray(Y_aug,&y2);CHKERRQ(ierr); 1511 /* Copy over inactive Y_aug to Y */ 1512 j1 = 0; 1513 for(i1=Y->map->rstart;i1 < Y->map->rend;i1++) { 1514 if(j1 < nkept && idx_actkept[j1] == i1) j1++; 1515 else y1[i1-Y->map->rstart] = y2[i1-Y->map->rstart]; 1516 } 1517 ierr = VecRestoreArray(Y,&y1);CHKERRQ(ierr); 1518 ierr = VecRestoreArray(Y_aug,&y2);CHKERRQ(ierr); 1519 1520 ierr = VecDestroy(&F_aug);CHKERRQ(ierr); 1521 ierr = VecDestroy(&Y_aug);CHKERRQ(ierr); 1522 ierr = MatDestroy(&J_aug);CHKERRQ(ierr); 1523 ierr = PetscFree(idx_actkept);CHKERRQ(ierr); 1524 } 1525 1526 if (!isequal) { 1527 ierr = ISDestroy(&vi->IS_inact_prev);CHKERRQ(ierr); 1528 ierr = ISDuplicate(IS_inact,&vi->IS_inact_prev);CHKERRQ(ierr); 1529 } 1530 ierr = ISDestroy(&IS_inact);CHKERRQ(ierr); 1531 1532 ierr = KSPGetIterationNumber(snes->ksp,&lits);CHKERRQ(ierr); 1533 snes->linear_its += lits; 1534 ierr = PetscInfo2(snes,"iter=%D, linear solve iterations=%D\n",snes->iter,lits);CHKERRQ(ierr); 1535 /* 1536 if (vi->precheckstep) { 1537 PetscBool changed_y = PETSC_FALSE; 1538 ierr = (*vi->precheckstep)(snes,X,Y,vi->precheck,&changed_y);CHKERRQ(ierr); 1539 } 1540 1541 if (PetscLogPrintInfo){ 1542 ierr = SNESVICheckResidual_Private(snes,snes->jacobian,F,Y,G,W);CHKERRQ(ierr); 1543 } 1544 */ 1545 /* Compute a (scaled) negative update in the line search routine: 1546 Y <- X - lambda*Y 1547 and evaluate G = function(Y) (depends on the line search). 1548 */ 1549 ierr = VecCopy(Y,snes->vec_sol_update);CHKERRQ(ierr); 1550 ynorm = 1; gnorm = fnorm; 1551 ierr = (*vi->LineSearch)(snes,vi->lsP,X,F,G,Y,W,fnorm,xnorm,&ynorm,&gnorm,&lssucceed);CHKERRQ(ierr); 1552 ierr = PetscInfo4(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lssucceed=%d\n",fnorm,gnorm,ynorm,(int)lssucceed);CHKERRQ(ierr); 1553 if (snes->reason == SNES_DIVERGED_FUNCTION_COUNT) break; 1554 if (snes->domainerror) { 1555 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 1556 PetscFunctionReturn(0); 1557 } 1558 if (!lssucceed) { 1559 if (++snes->numFailures >= snes->maxFailures) { 1560 PetscBool ismin; 1561 snes->reason = SNES_DIVERGED_LINE_SEARCH; 1562 ierr = SNESVICheckLocalMin_Private(snes,snes->jacobian,G,W,gnorm,&ismin);CHKERRQ(ierr); 1563 if (ismin) snes->reason = SNES_DIVERGED_LOCAL_MIN; 1564 break; 1565 } 1566 } 1567 /* Update function and solution vectors */ 1568 fnorm = gnorm; 1569 ierr = VecCopy(G,F);CHKERRQ(ierr); 1570 ierr = VecCopy(W,X);CHKERRQ(ierr); 1571 /* Monitor convergence */ 1572 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 1573 snes->iter = i+1; 1574 snes->norm = fnorm; 1575 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 1576 SNESLogConvHistory(snes,snes->norm,lits); 1577 ierr = SNESMonitor(snes,snes->iter,snes->norm);CHKERRQ(ierr); 1578 /* Test for convergence, xnorm = || X || */ 1579 if (snes->ops->converged != SNESSkipConverged) { ierr = VecNorm(X,NORM_2,&xnorm);CHKERRQ(ierr); } 1580 ierr = (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 1581 if (snes->reason) break; 1582 } 1583 if (i == maxits) { 1584 ierr = PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",maxits);CHKERRQ(ierr); 1585 if(!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT; 1586 } 1587 PetscFunctionReturn(0); 1588 } 1589 1590 /* -------------------------------------------------------------------------- */ 1591 /* 1592 SNESSetUp_VI - Sets up the internal data structures for the later use 1593 of the SNESVI nonlinear solver. 1594 1595 Input Parameter: 1596 . snes - the SNES context 1597 . x - the solution vector 1598 1599 Application Interface Routine: SNESSetUp() 1600 1601 Notes: 1602 For basic use of the SNES solvers, the user need not explicitly call 1603 SNESSetUp(), since these actions will automatically occur during 1604 the call to SNESSolve(). 1605 */ 1606 #undef __FUNCT__ 1607 #define __FUNCT__ "SNESSetUp_VI" 1608 PetscErrorCode SNESSetUp_VI(SNES snes) 1609 { 1610 PetscErrorCode ierr; 1611 SNES_VI *vi = (SNES_VI*) snes->data; 1612 PetscInt i_start[3],i_end[3]; 1613 1614 PetscFunctionBegin; 1615 1616 ierr = SNESDefaultGetWork(snes,3);CHKERRQ(ierr); 1617 1618 /* If the lower and upper bound on variables are not set, set it to 1619 -Inf and Inf */ 1620 if (!vi->xl && !vi->xu) { 1621 vi->usersetxbounds = PETSC_FALSE; 1622 ierr = VecDuplicate(snes->vec_sol, &vi->xl); CHKERRQ(ierr); 1623 ierr = VecSet(vi->xl,SNES_VI_NINF);CHKERRQ(ierr); 1624 ierr = VecDuplicate(snes->vec_sol, &vi->xu); CHKERRQ(ierr); 1625 ierr = VecSet(vi->xu,SNES_VI_INF);CHKERRQ(ierr); 1626 } else { 1627 /* Check if lower bound, upper bound and solution vector distribution across the processors is identical */ 1628 ierr = VecGetOwnershipRange(snes->vec_sol,i_start,i_end);CHKERRQ(ierr); 1629 ierr = VecGetOwnershipRange(vi->xl,i_start+1,i_end+1);CHKERRQ(ierr); 1630 ierr = VecGetOwnershipRange(vi->xu,i_start+2,i_end+2);CHKERRQ(ierr); 1631 if ((i_start[0] != i_start[1]) || (i_start[0] != i_start[2]) || (i_end[0] != i_end[1]) || (i_end[0] != i_end[2])) 1632 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Distribution of lower bound, upper bound and the solution vector should be identical across all the processors."); 1633 } 1634 if (snes->ops->solve != SNESSolveVI_SS) { 1635 /* Set up previous active index set for the first snes solve 1636 vi->IS_inact_prev = 0,1,2,....N */ 1637 PetscInt *indices; 1638 PetscInt i,n,rstart,rend; 1639 1640 ierr = VecGetOwnershipRange(snes->vec_sol,&rstart,&rend);CHKERRQ(ierr); 1641 ierr = VecGetLocalSize(snes->vec_sol,&n);CHKERRQ(ierr); 1642 ierr = PetscMalloc(n*sizeof(PetscInt),&indices);CHKERRQ(ierr); 1643 for(i=0;i < n; i++) indices[i] = rstart + i; 1644 ierr = ISCreateGeneral(((PetscObject)snes)->comm,n,indices,PETSC_OWN_POINTER,&vi->IS_inact_prev); 1645 } 1646 1647 if (snes->ops->solve == SNESSolveVI_SS) { 1648 ierr = VecDuplicate(snes->vec_sol, &vi->dpsi);CHKERRQ(ierr); 1649 ierr = VecDuplicate(snes->vec_sol, &vi->phi); CHKERRQ(ierr); 1650 ierr = VecDuplicate(snes->vec_sol, &vi->Da); CHKERRQ(ierr); 1651 ierr = VecDuplicate(snes->vec_sol, &vi->Db); CHKERRQ(ierr); 1652 ierr = VecDuplicate(snes->vec_sol, &vi->z);CHKERRQ(ierr); 1653 ierr = VecDuplicate(snes->vec_sol, &vi->t); CHKERRQ(ierr); 1654 } 1655 PetscFunctionReturn(0); 1656 } 1657 /* -------------------------------------------------------------------------- */ 1658 /* 1659 SNESDestroy_VI - Destroys the private SNES_VI context that was created 1660 with SNESCreate_VI(). 1661 1662 Input Parameter: 1663 . snes - the SNES context 1664 1665 Application Interface Routine: SNESDestroy() 1666 */ 1667 #undef __FUNCT__ 1668 #define __FUNCT__ "SNESDestroy_VI" 1669 PetscErrorCode SNESDestroy_VI(SNES snes) 1670 { 1671 SNES_VI *vi = (SNES_VI*) snes->data; 1672 PetscErrorCode ierr; 1673 1674 PetscFunctionBegin; 1675 if (snes->ops->solve != SNESSolveVI_SS) { 1676 ierr = ISDestroy(&vi->IS_inact_prev); 1677 } 1678 1679 if (snes->ops->solve == SNESSolveVI_SS) { 1680 /* clear vectors */ 1681 ierr = VecDestroy(&vi->dpsi);CHKERRQ(ierr); 1682 ierr = VecDestroy(&vi->phi); CHKERRQ(ierr); 1683 ierr = VecDestroy(&vi->Da); CHKERRQ(ierr); 1684 ierr = VecDestroy(&vi->Db); CHKERRQ(ierr); 1685 ierr = VecDestroy(&vi->z); CHKERRQ(ierr); 1686 ierr = VecDestroy(&vi->t); CHKERRQ(ierr); 1687 } 1688 1689 if (!vi->usersetxbounds) { 1690 ierr = VecDestroy(&vi->xl); CHKERRQ(ierr); 1691 ierr = VecDestroy(&vi->xu); CHKERRQ(ierr); 1692 } 1693 1694 ierr = PetscViewerASCIIMonitorDestroy(&vi->lsmonitor);CHKERRQ(ierr); 1695 ierr = PetscFree(snes->data);CHKERRQ(ierr); 1696 1697 /* clear composed functions */ 1698 ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESLineSearchSet_C","",PETSC_NULL);CHKERRQ(ierr); 1699 ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESLineSearchSetMonitor_C","",PETSC_NULL);CHKERRQ(ierr); 1700 PetscFunctionReturn(0); 1701 } 1702 1703 /* -------------------------------------------------------------------------- */ 1704 #undef __FUNCT__ 1705 #define __FUNCT__ "SNESLineSearchNo_VI" 1706 1707 /* 1708 This routine does not actually do a line search but it takes a full newton 1709 step while ensuring that the new iterates remain within the constraints. 1710 1711 */ 1712 PetscErrorCode SNESLineSearchNo_VI(SNES snes,void *lsctx,Vec x,Vec f,Vec g,Vec y,Vec w,PetscReal fnorm,PetscReal xnorm,PetscReal *ynorm,PetscReal *gnorm,PetscBool *flag) 1713 { 1714 PetscErrorCode ierr; 1715 SNES_VI *vi = (SNES_VI*)snes->data; 1716 PetscBool changed_w = PETSC_FALSE,changed_y = PETSC_FALSE; 1717 1718 PetscFunctionBegin; 1719 *flag = PETSC_TRUE; 1720 ierr = PetscLogEventBegin(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1721 ierr = VecNorm(y,NORM_2,ynorm);CHKERRQ(ierr); /* ynorm = || y || */ 1722 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); /* w <- x - y */ 1723 ierr = SNESVIProjectOntoBounds(snes,w);CHKERRQ(ierr); 1724 if (vi->postcheckstep) { 1725 ierr = (*vi->postcheckstep)(snes,x,y,w,vi->postcheck,&changed_y,&changed_w);CHKERRQ(ierr); 1726 } 1727 if (changed_y) { 1728 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); /* w <- x - y */ 1729 ierr = SNESVIProjectOntoBounds(snes,w);CHKERRQ(ierr); 1730 } 1731 ierr = SNESVIProjectOntoBounds(snes,w);CHKERRQ(ierr); 1732 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1733 if (!snes->domainerror) { 1734 if (snes->ops->solve != SNESSolveVI_SS) { 1735 ierr = SNESVIComputeInactiveSetFnorm(snes,g,w,gnorm);CHKERRQ(ierr); 1736 } else { 1737 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); /* gnorm = || g || */ 1738 } 1739 if (PetscIsInfOrNanReal(*gnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1740 } 1741 if (vi->lsmonitor) { 1742 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: Using full step: fnorm %G gnorm %G\n",fnorm,*gnorm);CHKERRQ(ierr); 1743 } 1744 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1745 PetscFunctionReturn(0); 1746 } 1747 1748 /* -------------------------------------------------------------------------- */ 1749 #undef __FUNCT__ 1750 #define __FUNCT__ "SNESLineSearchNoNorms_VI" 1751 1752 /* 1753 This routine is a copy of SNESLineSearchNoNorms in snes/impls/ls/ls.c 1754 */ 1755 PetscErrorCode SNESLineSearchNoNorms_VI(SNES snes,void *lsctx,Vec x,Vec f,Vec g,Vec y,Vec w,PetscReal fnorm,PetscReal xnorm,PetscReal *ynorm,PetscReal *gnorm,PetscBool *flag) 1756 { 1757 PetscErrorCode ierr; 1758 SNES_VI *vi = (SNES_VI*)snes->data; 1759 PetscBool changed_w = PETSC_FALSE,changed_y = PETSC_FALSE; 1760 1761 PetscFunctionBegin; 1762 *flag = PETSC_TRUE; 1763 ierr = PetscLogEventBegin(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1764 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); /* w <- x - y */ 1765 ierr = SNESVIProjectOntoBounds(snes,w);CHKERRQ(ierr); 1766 if (vi->postcheckstep) { 1767 ierr = (*vi->postcheckstep)(snes,x,y,w,vi->postcheck,&changed_y,&changed_w);CHKERRQ(ierr); 1768 } 1769 if (changed_y) { 1770 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); /* w <- x - y */ 1771 ierr = SNESVIProjectOntoBounds(snes,w);CHKERRQ(ierr); 1772 } 1773 1774 /* don't evaluate function the last time through */ 1775 if (snes->iter < snes->max_its-1) { 1776 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1777 } 1778 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1779 PetscFunctionReturn(0); 1780 } 1781 1782 /* -------------------------------------------------------------------------- */ 1783 #undef __FUNCT__ 1784 #define __FUNCT__ "SNESLineSearchCubic_VI" 1785 /* 1786 This routine implements a cubic line search while doing a projection on the variable bounds 1787 */ 1788 PetscErrorCode SNESLineSearchCubic_VI(SNES snes,void *lsctx,Vec x,Vec f,Vec g,Vec y,Vec w,PetscReal fnorm,PetscReal xnorm,PetscReal *ynorm,PetscReal *gnorm,PetscBool *flag) 1789 { 1790 PetscReal initslope,lambdaprev,gnormprev,a,b,d,t1,t2,rellength; 1791 PetscReal minlambda,lambda,lambdatemp; 1792 #if defined(PETSC_USE_COMPLEX) 1793 PetscScalar cinitslope; 1794 #endif 1795 PetscErrorCode ierr; 1796 PetscInt count; 1797 SNES_VI *vi = (SNES_VI*)snes->data; 1798 PetscBool changed_w = PETSC_FALSE,changed_y = PETSC_FALSE; 1799 MPI_Comm comm; 1800 1801 PetscFunctionBegin; 1802 ierr = PetscObjectGetComm((PetscObject)snes,&comm);CHKERRQ(ierr); 1803 ierr = PetscLogEventBegin(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1804 *flag = PETSC_TRUE; 1805 1806 ierr = VecNorm(y,NORM_2,ynorm);CHKERRQ(ierr); 1807 if (*ynorm == 0.0) { 1808 if (vi->lsmonitor) { 1809 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: Initial direction and size is 0\n");CHKERRQ(ierr); 1810 } 1811 *gnorm = fnorm; 1812 ierr = VecCopy(x,w);CHKERRQ(ierr); 1813 ierr = VecCopy(f,g);CHKERRQ(ierr); 1814 *flag = PETSC_FALSE; 1815 goto theend1; 1816 } 1817 if (*ynorm > vi->maxstep) { /* Step too big, so scale back */ 1818 if (vi->lsmonitor) { 1819 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: Scaling step by %G old ynorm %G\n",vi->maxstep/(*ynorm),*ynorm);CHKERRQ(ierr); 1820 } 1821 ierr = VecScale(y,vi->maxstep/(*ynorm));CHKERRQ(ierr); 1822 *ynorm = vi->maxstep; 1823 } 1824 ierr = VecMaxPointwiseDivide(y,x,&rellength);CHKERRQ(ierr); 1825 minlambda = vi->minlambda/rellength; 1826 ierr = MatMult(snes->jacobian,y,w);CHKERRQ(ierr); 1827 #if defined(PETSC_USE_COMPLEX) 1828 ierr = VecDot(f,w,&cinitslope);CHKERRQ(ierr); 1829 initslope = PetscRealPart(cinitslope); 1830 #else 1831 ierr = VecDot(f,w,&initslope);CHKERRQ(ierr); 1832 #endif 1833 if (initslope > 0.0) initslope = -initslope; 1834 if (initslope == 0.0) initslope = -1.0; 1835 1836 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); 1837 ierr = SNESVIProjectOntoBounds(snes,w);CHKERRQ(ierr); 1838 if (snes->nfuncs >= snes->max_funcs) { 1839 ierr = PetscInfo(snes,"Exceeded maximum function evaluations, while checking full step length!\n");CHKERRQ(ierr); 1840 *flag = PETSC_FALSE; 1841 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 1842 goto theend1; 1843 } 1844 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1845 if (snes->ops->solve != SNESSolveVI_SS) { 1846 ierr = SNESVIComputeInactiveSetFnorm(snes,g,w,gnorm);CHKERRQ(ierr); 1847 } else { 1848 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); 1849 } 1850 if (snes->domainerror) { 1851 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1852 PetscFunctionReturn(0); 1853 } 1854 if (PetscIsInfOrNanReal(*gnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1855 ierr = PetscInfo2(snes,"Initial fnorm %G gnorm %G\n",fnorm,*gnorm);CHKERRQ(ierr); 1856 if (.5*(*gnorm)*(*gnorm) <= .5*fnorm*fnorm + vi->alpha*initslope) { /* Sufficient reduction */ 1857 if (vi->lsmonitor) { 1858 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: Using full step: fnorm %G gnorm %G\n",fnorm,*gnorm);CHKERRQ(ierr); 1859 } 1860 goto theend1; 1861 } 1862 1863 /* Fit points with quadratic */ 1864 lambda = 1.0; 1865 lambdatemp = -initslope/((*gnorm)*(*gnorm) - fnorm*fnorm - 2.0*initslope); 1866 lambdaprev = lambda; 1867 gnormprev = *gnorm; 1868 if (lambdatemp > .5*lambda) lambdatemp = .5*lambda; 1869 if (lambdatemp <= .1*lambda) lambda = .1*lambda; 1870 else lambda = lambdatemp; 1871 1872 ierr = VecWAXPY(w,-lambda,y,x);CHKERRQ(ierr); 1873 ierr = SNESVIProjectOntoBounds(snes,w);CHKERRQ(ierr); 1874 if (snes->nfuncs >= snes->max_funcs) { 1875 ierr = PetscInfo1(snes,"Exceeded maximum function evaluations, while attempting quadratic backtracking! %D \n",snes->nfuncs);CHKERRQ(ierr); 1876 *flag = PETSC_FALSE; 1877 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 1878 goto theend1; 1879 } 1880 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1881 if (snes->ops->solve != SNESSolveVI_SS) { 1882 ierr = SNESVIComputeInactiveSetFnorm(snes,g,w,gnorm);CHKERRQ(ierr); 1883 } else { 1884 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); 1885 } 1886 if (snes->domainerror) { 1887 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1888 PetscFunctionReturn(0); 1889 } 1890 if (PetscIsInfOrNanReal(*gnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1891 if (vi->lsmonitor) { 1892 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: gnorm after quadratic fit %G\n",*gnorm);CHKERRQ(ierr); 1893 } 1894 if (.5*(*gnorm)*(*gnorm) < .5*fnorm*fnorm + lambda*vi->alpha*initslope) { /* sufficient reduction */ 1895 if (vi->lsmonitor) { 1896 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: Quadratically determined step, lambda=%18.16e\n",lambda);CHKERRQ(ierr); 1897 } 1898 goto theend1; 1899 } 1900 1901 /* Fit points with cubic */ 1902 count = 1; 1903 while (PETSC_TRUE) { 1904 if (lambda <= minlambda) { 1905 if (vi->lsmonitor) { 1906 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: unable to find good step length! After %D tries \n",count);CHKERRQ(ierr); 1907 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, minlambda=%18.16e, lambda=%18.16e, initial slope=%18.16e\n",fnorm,*gnorm,*ynorm,minlambda,lambda,initslope);CHKERRQ(ierr); 1908 } 1909 *flag = PETSC_FALSE; 1910 break; 1911 } 1912 t1 = .5*((*gnorm)*(*gnorm) - fnorm*fnorm) - lambda*initslope; 1913 t2 = .5*(gnormprev*gnormprev - fnorm*fnorm) - lambdaprev*initslope; 1914 a = (t1/(lambda*lambda) - t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev); 1915 b = (-lambdaprev*t1/(lambda*lambda) + lambda*t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev); 1916 d = b*b - 3*a*initslope; 1917 if (d < 0.0) d = 0.0; 1918 if (a == 0.0) { 1919 lambdatemp = -initslope/(2.0*b); 1920 } else { 1921 lambdatemp = (-b + sqrt(d))/(3.0*a); 1922 } 1923 lambdaprev = lambda; 1924 gnormprev = *gnorm; 1925 if (lambdatemp > .5*lambda) lambdatemp = .5*lambda; 1926 if (lambdatemp <= .1*lambda) lambda = .1*lambda; 1927 else lambda = lambdatemp; 1928 ierr = VecWAXPY(w,-lambda,y,x);CHKERRQ(ierr); 1929 ierr = SNESVIProjectOntoBounds(snes,w);CHKERRQ(ierr); 1930 if (snes->nfuncs >= snes->max_funcs) { 1931 ierr = PetscInfo1(snes,"Exceeded maximum function evaluations, while looking for good step length! %D \n",count);CHKERRQ(ierr); 1932 ierr = PetscInfo5(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lambda=%18.16e, initial slope=%18.16e\n",fnorm,*gnorm,*ynorm,lambda,initslope);CHKERRQ(ierr); 1933 *flag = PETSC_FALSE; 1934 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 1935 break; 1936 } 1937 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1938 if (snes->ops->solve != SNESSolveVI_SS) { 1939 ierr = SNESVIComputeInactiveSetFnorm(snes,g,w,gnorm);CHKERRQ(ierr); 1940 } else { 1941 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); 1942 } 1943 if (snes->domainerror) { 1944 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1945 PetscFunctionReturn(0); 1946 } 1947 if (PetscIsInfOrNanReal(*gnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1948 if (.5*(*gnorm)*(*gnorm) < .5*fnorm*fnorm + lambda*vi->alpha*initslope) { /* is reduction enough? */ 1949 if (vi->lsmonitor) { 1950 ierr = PetscPrintf(comm," Line search: Cubically determined step, current gnorm %G lambda=%18.16e\n",*gnorm,lambda);CHKERRQ(ierr); 1951 } 1952 break; 1953 } else { 1954 if (vi->lsmonitor) { 1955 ierr = PetscPrintf(comm," Line search: Cubic step no good, shrinking lambda, current gnorm %G lambda=%18.16e\n",*gnorm,lambda);CHKERRQ(ierr); 1956 } 1957 } 1958 count++; 1959 } 1960 theend1: 1961 /* Optional user-defined check for line search step validity */ 1962 if (vi->postcheckstep && *flag) { 1963 ierr = (*vi->postcheckstep)(snes,x,y,w,vi->postcheck,&changed_y,&changed_w);CHKERRQ(ierr); 1964 if (changed_y) { 1965 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); 1966 ierr = SNESVIProjectOntoBounds(snes,w);CHKERRQ(ierr); 1967 } 1968 if (changed_y || changed_w) { /* recompute the function if the step has changed */ 1969 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1970 if (snes->ops->solve != SNESSolveVI_SS) { 1971 ierr = SNESVIComputeInactiveSetFnorm(snes,g,w,gnorm);CHKERRQ(ierr); 1972 } else { 1973 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); 1974 } 1975 if (snes->domainerror) { 1976 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1977 PetscFunctionReturn(0); 1978 } 1979 if (PetscIsInfOrNanReal(*gnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1980 ierr = VecNormBegin(y,NORM_2,ynorm);CHKERRQ(ierr); 1981 ierr = VecNormEnd(y,NORM_2,ynorm);CHKERRQ(ierr); 1982 } 1983 } 1984 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1985 PetscFunctionReturn(0); 1986 } 1987 1988 #undef __FUNCT__ 1989 #define __FUNCT__ "SNESLineSearchQuadratic_VI" 1990 /* 1991 This routine does a quadratic line search while keeping the iterates within the variable bounds 1992 */ 1993 PetscErrorCode SNESLineSearchQuadratic_VI(SNES snes,void *lsctx,Vec x,Vec f,Vec g,Vec y,Vec w,PetscReal fnorm,PetscReal xnorm,PetscReal *ynorm,PetscReal *gnorm,PetscBool *flag) 1994 { 1995 /* 1996 Note that for line search purposes we work with with the related 1997 minimization problem: 1998 min z(x): R^n -> R, 1999 where z(x) = .5 * fnorm*fnorm,and fnorm = || f ||_2. 2000 */ 2001 PetscReal initslope,minlambda,lambda,lambdatemp,rellength; 2002 #if defined(PETSC_USE_COMPLEX) 2003 PetscScalar cinitslope; 2004 #endif 2005 PetscErrorCode ierr; 2006 PetscInt count; 2007 SNES_VI *vi = (SNES_VI*)snes->data; 2008 PetscBool changed_w = PETSC_FALSE,changed_y = PETSC_FALSE; 2009 2010 PetscFunctionBegin; 2011 ierr = PetscLogEventBegin(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 2012 *flag = PETSC_TRUE; 2013 2014 ierr = VecNorm(y,NORM_2,ynorm);CHKERRQ(ierr); 2015 if (*ynorm == 0.0) { 2016 if (vi->lsmonitor) { 2017 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor,"Line search: Direction and size is 0\n");CHKERRQ(ierr); 2018 } 2019 *gnorm = fnorm; 2020 ierr = VecCopy(x,w);CHKERRQ(ierr); 2021 ierr = VecCopy(f,g);CHKERRQ(ierr); 2022 *flag = PETSC_FALSE; 2023 goto theend2; 2024 } 2025 if (*ynorm > vi->maxstep) { /* Step too big, so scale back */ 2026 ierr = VecScale(y,vi->maxstep/(*ynorm));CHKERRQ(ierr); 2027 *ynorm = vi->maxstep; 2028 } 2029 ierr = VecMaxPointwiseDivide(y,x,&rellength);CHKERRQ(ierr); 2030 minlambda = vi->minlambda/rellength; 2031 ierr = MatMult(snes->jacobian,y,w);CHKERRQ(ierr); 2032 #if defined(PETSC_USE_COMPLEX) 2033 ierr = VecDot(f,w,&cinitslope);CHKERRQ(ierr); 2034 initslope = PetscRealPart(cinitslope); 2035 #else 2036 ierr = VecDot(f,w,&initslope);CHKERRQ(ierr); 2037 #endif 2038 if (initslope > 0.0) initslope = -initslope; 2039 if (initslope == 0.0) initslope = -1.0; 2040 ierr = PetscInfo1(snes,"Initslope %G \n",initslope);CHKERRQ(ierr); 2041 2042 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); 2043 ierr = SNESVIProjectOntoBounds(snes,w);CHKERRQ(ierr); 2044 if (snes->nfuncs >= snes->max_funcs) { 2045 ierr = PetscInfo(snes,"Exceeded maximum function evaluations, while checking full step length!\n");CHKERRQ(ierr); 2046 *flag = PETSC_FALSE; 2047 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 2048 goto theend2; 2049 } 2050 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 2051 if (snes->ops->solve != SNESSolveVI_SS) { 2052 ierr = SNESVIComputeInactiveSetFnorm(snes,g,w,gnorm);CHKERRQ(ierr); 2053 } else { 2054 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); 2055 } 2056 if (snes->domainerror) { 2057 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 2058 PetscFunctionReturn(0); 2059 } 2060 if (PetscIsInfOrNanReal(*gnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 2061 if (.5*(*gnorm)*(*gnorm) <= .5*fnorm*fnorm + vi->alpha*initslope) { /* Sufficient reduction */ 2062 if (vi->lsmonitor) { 2063 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: Using full step: fnorm %G gnorm %G\n",fnorm,*gnorm);CHKERRQ(ierr); 2064 } 2065 goto theend2; 2066 } 2067 2068 /* Fit points with quadratic */ 2069 lambda = 1.0; 2070 count = 1; 2071 while (PETSC_TRUE) { 2072 if (lambda <= minlambda) { /* bad luck; use full step */ 2073 if (vi->lsmonitor) { 2074 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor,"Line search: Unable to find good step length! %D \n",count);CHKERRQ(ierr); 2075 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor,"Line search: fnorm=%G, gnorm=%G, ynorm=%G, lambda=%G, initial slope=%G\n",fnorm,*gnorm,*ynorm,lambda,initslope);CHKERRQ(ierr); 2076 } 2077 ierr = VecCopy(x,w);CHKERRQ(ierr); 2078 *flag = PETSC_FALSE; 2079 break; 2080 } 2081 lambdatemp = -initslope/((*gnorm)*(*gnorm) - fnorm*fnorm - 2.0*initslope); 2082 if (lambdatemp > .5*lambda) lambdatemp = .5*lambda; 2083 if (lambdatemp <= .1*lambda) lambda = .1*lambda; 2084 else lambda = lambdatemp; 2085 2086 ierr = VecWAXPY(w,-lambda,y,x);CHKERRQ(ierr); 2087 ierr = SNESVIProjectOntoBounds(snes,w);CHKERRQ(ierr); 2088 if (snes->nfuncs >= snes->max_funcs) { 2089 ierr = PetscInfo1(snes,"Exceeded maximum function evaluations, while looking for good step length! %D \n",count);CHKERRQ(ierr); 2090 ierr = PetscInfo5(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lambda=%18.16e, initial slope=%18.16e\n",fnorm,*gnorm,*ynorm,lambda,initslope);CHKERRQ(ierr); 2091 *flag = PETSC_FALSE; 2092 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 2093 break; 2094 } 2095 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 2096 if (snes->domainerror) { 2097 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 2098 PetscFunctionReturn(0); 2099 } 2100 if (snes->ops->solve != SNESSolveVI_SS) { 2101 ierr = SNESVIComputeInactiveSetFnorm(snes,g,w,gnorm);CHKERRQ(ierr); 2102 } else { 2103 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); 2104 } 2105 if (PetscIsInfOrNanReal(*gnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 2106 if (.5*(*gnorm)*(*gnorm) < .5*fnorm*fnorm + lambda*vi->alpha*initslope) { /* sufficient reduction */ 2107 if (vi->lsmonitor) { 2108 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line Search: Quadratically determined step, lambda=%G\n",lambda);CHKERRQ(ierr); 2109 } 2110 break; 2111 } 2112 count++; 2113 } 2114 theend2: 2115 /* Optional user-defined check for line search step validity */ 2116 if (vi->postcheckstep) { 2117 ierr = (*vi->postcheckstep)(snes,x,y,w,vi->postcheck,&changed_y,&changed_w);CHKERRQ(ierr); 2118 if (changed_y) { 2119 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); 2120 ierr = SNESVIProjectOntoBounds(snes,w);CHKERRQ(ierr); 2121 } 2122 if (changed_y || changed_w) { /* recompute the function if the step has changed */ 2123 ierr = SNESComputeFunction(snes,w,g); 2124 if (snes->domainerror) { 2125 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 2126 PetscFunctionReturn(0); 2127 } 2128 if (snes->ops->solve != SNESSolveVI_SS) { 2129 ierr = SNESVIComputeInactiveSetFnorm(snes,g,w,gnorm);CHKERRQ(ierr); 2130 } else { 2131 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); 2132 } 2133 2134 ierr = VecNormBegin(y,NORM_2,ynorm);CHKERRQ(ierr); 2135 ierr = VecNormEnd(y,NORM_2,ynorm);CHKERRQ(ierr); 2136 if (PetscIsInfOrNanReal(*gnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 2137 } 2138 } 2139 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 2140 PetscFunctionReturn(0); 2141 } 2142 2143 typedef PetscErrorCode (*FCN2)(SNES,void*,Vec,Vec,Vec,Vec,Vec,PetscReal,PetscReal,PetscReal*,PetscReal*,PetscBool*); /* force argument to next function to not be extern C*/ 2144 /* -------------------------------------------------------------------------- */ 2145 EXTERN_C_BEGIN 2146 #undef __FUNCT__ 2147 #define __FUNCT__ "SNESLineSearchSet_VI" 2148 PetscErrorCode SNESLineSearchSet_VI(SNES snes,FCN2 func,void *lsctx) 2149 { 2150 PetscFunctionBegin; 2151 ((SNES_VI *)(snes->data))->LineSearch = func; 2152 ((SNES_VI *)(snes->data))->lsP = lsctx; 2153 PetscFunctionReturn(0); 2154 } 2155 EXTERN_C_END 2156 2157 /* -------------------------------------------------------------------------- */ 2158 EXTERN_C_BEGIN 2159 #undef __FUNCT__ 2160 #define __FUNCT__ "SNESLineSearchSetMonitor_VI" 2161 PetscErrorCode SNESLineSearchSetMonitor_VI(SNES snes,PetscBool flg) 2162 { 2163 SNES_VI *vi = (SNES_VI*)snes->data; 2164 PetscErrorCode ierr; 2165 2166 PetscFunctionBegin; 2167 if (flg && !vi->lsmonitor) { 2168 ierr = PetscViewerASCIIMonitorCreate(((PetscObject)snes)->comm,"stdout",((PetscObject)snes)->tablevel,&vi->lsmonitor);CHKERRQ(ierr); 2169 } else if (!flg && vi->lsmonitor) { 2170 ierr = PetscViewerASCIIMonitorDestroy(&vi->lsmonitor);CHKERRQ(ierr); 2171 } 2172 PetscFunctionReturn(0); 2173 } 2174 EXTERN_C_END 2175 2176 /* 2177 SNESView_VI - Prints info from the SNESVI data structure. 2178 2179 Input Parameters: 2180 . SNES - the SNES context 2181 . viewer - visualization context 2182 2183 Application Interface Routine: SNESView() 2184 */ 2185 #undef __FUNCT__ 2186 #define __FUNCT__ "SNESView_VI" 2187 static PetscErrorCode SNESView_VI(SNES snes,PetscViewer viewer) 2188 { 2189 SNES_VI *vi = (SNES_VI *)snes->data; 2190 const char *cstr,*tstr; 2191 PetscErrorCode ierr; 2192 PetscBool iascii; 2193 2194 PetscFunctionBegin; 2195 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 2196 if (iascii) { 2197 if (vi->LineSearch == SNESLineSearchNo_VI) cstr = "SNESLineSearchNo"; 2198 else if (vi->LineSearch == SNESLineSearchQuadratic_VI) cstr = "SNESLineSearchQuadratic"; 2199 else if (vi->LineSearch == SNESLineSearchCubic_VI) cstr = "SNESLineSearchCubic"; 2200 else cstr = "unknown"; 2201 if (snes->ops->solve == SNESSolveVI_SS) tstr = "Semismooth"; 2202 else if (snes->ops->solve == SNESSolveVI_RS) tstr = "Reduced Space"; 2203 else if (snes->ops->solve == SNESSolveVI_RSAUG) tstr = "Reduced space with augmented variables"; 2204 else tstr = "unknown"; 2205 ierr = PetscViewerASCIIPrintf(viewer," VI algorithm: %s\n",tstr);CHKERRQ(ierr); 2206 ierr = PetscViewerASCIIPrintf(viewer," line search variant: %s\n",cstr);CHKERRQ(ierr); 2207 ierr = PetscViewerASCIIPrintf(viewer," alpha=%G, maxstep=%G, minlambda=%G\n",vi->alpha,vi->maxstep,vi->minlambda);CHKERRQ(ierr); 2208 } else { 2209 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type %s not supported for SNES EQ VI",((PetscObject)viewer)->type_name); 2210 } 2211 PetscFunctionReturn(0); 2212 } 2213 2214 #undef __FUNCT__ 2215 #define __FUNCT__ "SNESVISetVariableBounds" 2216 /*@ 2217 SNESVISetVariableBounds - Sets the lower and upper bounds for the solution vector. xl <= x <= xu. 2218 2219 Input Parameters: 2220 . snes - the SNES context. 2221 . xl - lower bound. 2222 . xu - upper bound. 2223 2224 Notes: 2225 If this routine is not called then the lower and upper bounds are set to 2226 SNES_VI_INF and SNES_VI_NINF respectively during SNESSetUp(). 2227 2228 @*/ 2229 PetscErrorCode SNESVISetVariableBounds(SNES snes, Vec xl, Vec xu) 2230 { 2231 SNES_VI *vi; 2232 PetscErrorCode ierr; 2233 2234 PetscFunctionBegin; 2235 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2236 PetscValidHeaderSpecific(xl,VEC_CLASSID,2); 2237 PetscValidHeaderSpecific(xu,VEC_CLASSID,3); 2238 if (!snes->vec_func) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call SNESSetFunction() first"); 2239 if (xl->map->N != snes->vec_func->map->N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Incompatible vector lengths lower bound = %D solution vector = %D",xl->map->N,snes->vec_func->map->N); 2240 if (xu->map->N != snes->vec_func->map->N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Incompatible vector lengths: upper bound = %D solution vector = %D",xu->map->N,snes->vec_func->map->N); 2241 ierr = SNESSetType(snes,SNESVI);CHKERRQ(ierr); 2242 vi = (SNES_VI*)snes->data; 2243 vi->usersetxbounds = PETSC_TRUE; 2244 vi->xl = xl; 2245 vi->xu = xu; 2246 PetscFunctionReturn(0); 2247 } 2248 2249 /* -------------------------------------------------------------------------- */ 2250 /* 2251 SNESSetFromOptions_VI - Sets various parameters for the SNESVI method. 2252 2253 Input Parameter: 2254 . snes - the SNES context 2255 2256 Application Interface Routine: SNESSetFromOptions() 2257 */ 2258 #undef __FUNCT__ 2259 #define __FUNCT__ "SNESSetFromOptions_VI" 2260 static PetscErrorCode SNESSetFromOptions_VI(SNES snes) 2261 { 2262 SNES_VI *vi = (SNES_VI *)snes->data; 2263 const char *lses[] = {"basic","basicnonorms","quadratic","cubic"}; 2264 const char *vies[] = {"ss","rs","rsaug"}; 2265 PetscErrorCode ierr; 2266 PetscInt indx; 2267 PetscBool flg,set,flg2; 2268 2269 PetscFunctionBegin; 2270 ierr = PetscOptionsHead("SNES semismooth method options");CHKERRQ(ierr); 2271 ierr = PetscOptionsBool("-snes_vi_monitor","Monitor all non-active variables","None",PETSC_FALSE,&flg,0);CHKERRQ(ierr); 2272 if (flg) { 2273 ierr = SNESMonitorSet(snes,SNESMonitorVI,0,0);CHKERRQ(ierr); 2274 } 2275 ierr = PetscOptionsReal("-snes_ls_alpha","Function norm must decrease by","None",vi->alpha,&vi->alpha,0);CHKERRQ(ierr); 2276 ierr = PetscOptionsReal("-snes_ls_maxstep","Step must be less than","None",vi->maxstep,&vi->maxstep,0);CHKERRQ(ierr); 2277 ierr = PetscOptionsReal("-snes_ls_minlambda","Minimum lambda allowed","None",vi->minlambda,&vi->minlambda,0);CHKERRQ(ierr); 2278 ierr = PetscOptionsReal("-snes_vi_const_tol","constraint tolerance","None",vi->const_tol,&vi->const_tol,0);CHKERRQ(ierr); 2279 ierr = PetscOptionsBool("-snes_ls_monitor","Print progress of line searches","SNESLineSearchSetMonitor",vi->lsmonitor ? PETSC_TRUE : PETSC_FALSE,&flg,&set);CHKERRQ(ierr); 2280 if (set) {ierr = SNESLineSearchSetMonitor(snes,flg);CHKERRQ(ierr);} 2281 ierr = PetscOptionsEList("-snes_vi_type","Semismooth algorithm used","",vies,3,"ss",&indx,&flg2);CHKERRQ(ierr); 2282 if (flg2) { 2283 switch (indx) { 2284 case 0: 2285 snes->ops->solve = SNESSolveVI_SS; 2286 break; 2287 case 1: 2288 snes->ops->solve = SNESSolveVI_RS; 2289 break; 2290 case 2: 2291 snes->ops->solve = SNESSolveVI_RSAUG; 2292 } 2293 } 2294 ierr = PetscOptionsEList("-snes_ls","Line search used","SNESLineSearchSet",lses,4,"basic",&indx,&flg);CHKERRQ(ierr); 2295 if (flg) { 2296 switch (indx) { 2297 case 0: 2298 ierr = SNESLineSearchSet(snes,SNESLineSearchNo_VI,PETSC_NULL);CHKERRQ(ierr); 2299 break; 2300 case 1: 2301 ierr = SNESLineSearchSet(snes,SNESLineSearchNoNorms_VI,PETSC_NULL);CHKERRQ(ierr); 2302 break; 2303 case 2: 2304 ierr = SNESLineSearchSet(snes,SNESLineSearchQuadratic_VI,PETSC_NULL);CHKERRQ(ierr); 2305 break; 2306 case 3: 2307 ierr = SNESLineSearchSet(snes,SNESLineSearchCubic_VI,PETSC_NULL);CHKERRQ(ierr); 2308 break; 2309 } 2310 } 2311 ierr = PetscOptionsTail();CHKERRQ(ierr); 2312 PetscFunctionReturn(0); 2313 } 2314 /* -------------------------------------------------------------------------- */ 2315 /*MC 2316 SNESVI - Various solvers for variational inequalities based on Newton's method 2317 2318 Options Database: 2319 + -snes_vi_type <ss,rs,rsaug> a semi-smooth solver, a reduced space active set method and a reduced space active set method that does not eliminate the active constraints from the Jacobian instead augments the Jacobian with 2320 additional variables that enforce the constraints 2321 - -snes_vi_monitor - prints the number of active constraints at each iteration. 2322 2323 2324 Level: beginner 2325 2326 .seealso: SNESCreate(), SNES, SNESSetType(), SNESTR, SNESLineSearchSet(), 2327 SNESLineSearchSetPostCheck(), SNESLineSearchNo(), SNESLineSearchCubic(), SNESLineSearchQuadratic(), 2328 SNESLineSearchSet(), SNESLineSearchNoNorms(), SNESLineSearchSetPreCheck(), SNESLineSearchSetParams(), SNESLineSearchGetParams() 2329 2330 M*/ 2331 EXTERN_C_BEGIN 2332 #undef __FUNCT__ 2333 #define __FUNCT__ "SNESCreate_VI" 2334 PetscErrorCode SNESCreate_VI(SNES snes) 2335 { 2336 PetscErrorCode ierr; 2337 SNES_VI *vi; 2338 2339 PetscFunctionBegin; 2340 snes->ops->setup = SNESSetUp_VI; 2341 snes->ops->solve = SNESSolveVI_RS; 2342 snes->ops->destroy = SNESDestroy_VI; 2343 snes->ops->setfromoptions = SNESSetFromOptions_VI; 2344 snes->ops->view = SNESView_VI; 2345 snes->ops->reset = 0; /* XXX Implement!!! */ 2346 snes->ops->converged = SNESDefaultConverged_VI; 2347 2348 ierr = PetscNewLog(snes,SNES_VI,&vi);CHKERRQ(ierr); 2349 snes->data = (void*)vi; 2350 vi->alpha = 1.e-4; 2351 vi->maxstep = 1.e8; 2352 vi->minlambda = 1.e-12; 2353 vi->LineSearch = SNESLineSearchNo_VI; 2354 vi->lsP = PETSC_NULL; 2355 vi->postcheckstep = PETSC_NULL; 2356 vi->postcheck = PETSC_NULL; 2357 vi->precheckstep = PETSC_NULL; 2358 vi->precheck = PETSC_NULL; 2359 vi->const_tol = 2.2204460492503131e-16; 2360 vi->checkredundancy = PETSC_NULL; 2361 2362 ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESLineSearchSetMonitor_C","SNESLineSearchSetMonitor_VI",SNESLineSearchSetMonitor_VI);CHKERRQ(ierr); 2363 ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESLineSearchSet_C","SNESLineSearchSet_VI",SNESLineSearchSet_VI);CHKERRQ(ierr); 2364 2365 PetscFunctionReturn(0); 2366 } 2367 EXTERN_C_END 2368