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