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 PetscFunctionReturn(0); 141 } 142 143 extern PetscErrorCode DMSetVI(DM,Vec,Vec,Vec,Vec,IS); 144 145 #undef __FUNCT__ 146 #define __FUNCT__ "DMCoarsen_SNESVI" 147 /* 148 DMCoarsen_SNESVI - Computes the regular coarsened DM then computes additional information about its inactive set 149 150 */ 151 PetscErrorCode DMCoarsen_SNESVI(DM dm1,MPI_Comm comm,DM *dm2) 152 { 153 PetscErrorCode ierr; 154 PetscContainer isnes; 155 DM_SNESVI *dmsnesvi1; 156 Vec upper,lower,values,F; 157 IS inactive; 158 VecScatter inject; 159 160 PetscFunctionBegin; 161 ierr = PetscObjectQuery((PetscObject)dm1,"VI",(PetscObject *)&isnes);CHKERRQ(ierr); 162 if (!isnes) SETERRQ(((PetscObject)dm1)->comm,PETSC_ERR_PLIB,"Composed VI data structure is missing"); 163 ierr = PetscContainerGetPointer(isnes,(void**)&dmsnesvi1);CHKERRQ(ierr); 164 165 /* get the original coarsen */ 166 ierr = (*dmsnesvi1->coarsen)(dm1,comm,dm2);CHKERRQ(ierr); 167 168 /* not sure why this extra reference is needed, but without the dm2 disappears too early */ 169 ierr = PetscObjectReference((PetscObject)*dm2);CHKERRQ(ierr); 170 171 /* need to set back global vectors in order to use the original injection */ 172 ierr = DMClearGlobalVectors(dm1);CHKERRQ(ierr); 173 dm1->ops->createglobalvector = dmsnesvi1->createglobalvector; 174 ierr = DMGetInjection(*dm2,dm1,&inject);CHKERRQ(ierr); 175 ierr = DMCreateGlobalVector(*dm2,&upper);CHKERRQ(ierr); 176 ierr = VecScatterBegin(inject,dmsnesvi1->upper,upper,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 177 ierr = VecScatterEnd(inject,dmsnesvi1->upper,upper,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 178 ierr = DMCreateGlobalVector(*dm2,&lower);CHKERRQ(ierr); 179 ierr = VecScatterBegin(inject,dmsnesvi1->lower,lower,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 180 ierr = VecScatterEnd(inject,dmsnesvi1->lower,lower,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 181 ierr = DMCreateGlobalVector(*dm2,&values);CHKERRQ(ierr); 182 ierr = VecScatterBegin(inject,dmsnesvi1->values,values,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 183 ierr = VecScatterEnd(inject,dmsnesvi1->values,values,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 184 ierr = DMCreateGlobalVector(*dm2,&F);CHKERRQ(ierr); 185 ierr = VecScatterBegin(inject,dmsnesvi1->F,F,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 186 ierr = VecScatterEnd(inject,dmsnesvi1->F,F,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 187 ierr = VecScatterDestroy(&inject);CHKERRQ(ierr); 188 ierr = DMClearGlobalVectors(dm1);CHKERRQ(ierr); 189 dm1->ops->createglobalvector = DMCreateGlobalVector_SNESVI; 190 ierr = SNESVIGetInActiveSetIS(upper,lower,values,F,&inactive);CHKERRQ(ierr); 191 ierr = DMSetVI(*dm2,upper,lower,values,F,inactive);CHKERRQ(ierr); 192 ierr = VecDestroy(&upper);CHKERRQ(ierr); 193 ierr = VecDestroy(&lower);CHKERRQ(ierr); 194 ierr = VecDestroy(&values);CHKERRQ(ierr); 195 ierr = VecDestroy(&F);CHKERRQ(ierr); 196 ierr = ISDestroy(&inactive);CHKERRQ(ierr); 197 PetscFunctionReturn(0); 198 } 199 200 #undef __FUNCT__ 201 #define __FUNCT__ "DMDestroy_SNESVI" 202 PetscErrorCode DMDestroy_SNESVI(DM_SNESVI *dmsnesvi) 203 { 204 PetscErrorCode ierr; 205 206 PetscFunctionBegin; 207 /* reset the base methods in the DM object that were changed when the DM_SNESVI was reset */ 208 dmsnesvi->dm->ops->getinterpolation = dmsnesvi->getinterpolation; 209 dmsnesvi->dm->ops->coarsen = dmsnesvi->coarsen; 210 dmsnesvi->dm->ops->createglobalvector = dmsnesvi->createglobalvector; 211 212 ierr = VecDestroy(&dmsnesvi->upper);CHKERRQ(ierr); 213 ierr = VecDestroy(&dmsnesvi->lower);CHKERRQ(ierr); 214 ierr = VecDestroy(&dmsnesvi->values);CHKERRQ(ierr); 215 ierr = VecDestroy(&dmsnesvi->F);CHKERRQ(ierr); 216 ierr = ISDestroy(&dmsnesvi->inactive);CHKERRQ(ierr); 217 ierr = PetscFree(dmsnesvi);CHKERRQ(ierr); 218 PetscFunctionReturn(0); 219 } 220 221 #undef __FUNCT__ 222 #define __FUNCT__ "DMSetVI" 223 /* 224 DMSetVI - Marks a DM as associated with a VI problem. This causes the interpolation/restriction operators to 225 be restricted to only those variables NOT associated with active constraints. 226 227 */ 228 PetscErrorCode DMSetVI(DM dm,Vec upper,Vec lower,Vec values,Vec F,IS inactive) 229 { 230 PetscErrorCode ierr; 231 PetscContainer isnes; 232 DM_SNESVI *dmsnesvi; 233 PetscInt Ni,Nu; 234 235 PetscFunctionBegin; 236 if (!dm) PetscFunctionReturn(0); 237 ierr = ISGetSize(inactive,&Ni);CHKERRQ(ierr); 238 ierr = VecGetSize(upper,&Nu);CHKERRQ(ierr); 239 if (Ni == Nu) 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 } 1186 1187 /* ierr = ISView(IS_inact,0);CHKERRQ(ierr); */ 1188 /* ierr = ISView(IS_act,0);CHKERRQ(ierr);*/ 1189 /* ierr = MatView(snes->jacobian_pre,0); */ 1190 1191 ierr = KSPSetOperators(snes->ksp,jac_inact_inact,prejac_inact_inact,flg);CHKERRQ(ierr); 1192 ierr = KSPSetUp(snes->ksp);CHKERRQ(ierr); 1193 { 1194 PC pc; 1195 PetscBool flg; 1196 ierr = KSPGetPC(snes->ksp,&pc);CHKERRQ(ierr); 1197 ierr = PetscTypeCompare((PetscObject)pc,PCFIELDSPLIT,&flg);CHKERRQ(ierr); 1198 if (flg) { 1199 KSP *subksps; 1200 ierr = PCFieldSplitGetSubKSP(pc,PETSC_NULL,&subksps);CHKERRQ(ierr); 1201 ierr = KSPGetPC(subksps[0],&pc);CHKERRQ(ierr); 1202 ierr = PetscFree(subksps);CHKERRQ(ierr); 1203 ierr = PetscTypeCompare((PetscObject)pc,PCBJACOBI,&flg);CHKERRQ(ierr); 1204 if (flg) { 1205 PetscInt n,N = 101*101,j,cnts[3] = {0,0,0}; 1206 const PetscInt *ii; 1207 1208 ierr = ISGetSize(IS_inact,&n);CHKERRQ(ierr); 1209 ierr = ISGetIndices(IS_inact,&ii);CHKERRQ(ierr); 1210 for (j=0; j<n; j++) { 1211 if (ii[j] < N) cnts[0]++; 1212 else if (ii[j] < 2*N) cnts[1]++; 1213 else if (ii[j] < 3*N) cnts[2]++; 1214 } 1215 ierr = ISRestoreIndices(IS_inact,&ii);CHKERRQ(ierr); 1216 1217 ierr = PCBJacobiSetTotalBlocks(pc,3,cnts);CHKERRQ(ierr); 1218 } 1219 } 1220 } 1221 1222 ierr = SNES_KSPSolve(snes,snes->ksp,F_inact,Y_inact);CHKERRQ(ierr); 1223 ierr = KSPGetConvergedReason(snes->ksp,&kspreason);CHKERRQ(ierr); 1224 if (kspreason < 0) { 1225 if (++snes->numLinearSolveFailures >= snes->maxLinearSolveFailures) { 1226 ierr = PetscInfo2(snes,"iter=%D, number linear solve failures %D greater than current SNES allowed, stopping solve\n",snes->iter,snes->numLinearSolveFailures);CHKERRQ(ierr); 1227 snes->reason = SNES_DIVERGED_LINEAR_SOLVE; 1228 break; 1229 } 1230 } 1231 1232 ierr = VecScatterBegin(scat_act,Y_act,Y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1233 ierr = VecScatterEnd(scat_act,Y_act,Y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1234 ierr = VecScatterBegin(scat_inact,Y_inact,Y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1235 ierr = VecScatterEnd(scat_inact,Y_inact,Y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1236 1237 ierr = VecDestroy(&F_inact);CHKERRQ(ierr); 1238 ierr = VecDestroy(&Y_act);CHKERRQ(ierr); 1239 ierr = VecDestroy(&Y_inact);CHKERRQ(ierr); 1240 ierr = VecScatterDestroy(&scat_act);CHKERRQ(ierr); 1241 ierr = VecScatterDestroy(&scat_inact);CHKERRQ(ierr); 1242 ierr = ISDestroy(&IS_act);CHKERRQ(ierr); 1243 if (!isequal) { 1244 ierr = ISDestroy(&vi->IS_inact_prev);CHKERRQ(ierr); 1245 ierr = ISDuplicate(IS_inact,&vi->IS_inact_prev);CHKERRQ(ierr); 1246 } 1247 ierr = ISDestroy(&IS_inact);CHKERRQ(ierr); 1248 ierr = MatDestroy(&jac_inact_inact);CHKERRQ(ierr); 1249 if (snes->jacobian != snes->jacobian_pre) { 1250 ierr = MatDestroy(&prejac_inact_inact);CHKERRQ(ierr); 1251 } 1252 ierr = KSPGetIterationNumber(snes->ksp,&lits);CHKERRQ(ierr); 1253 snes->linear_its += lits; 1254 ierr = PetscInfo2(snes,"iter=%D, linear solve iterations=%D\n",snes->iter,lits);CHKERRQ(ierr); 1255 /* 1256 if (vi->precheckstep) { 1257 PetscBool changed_y = PETSC_FALSE; 1258 ierr = (*vi->precheckstep)(snes,X,Y,vi->precheck,&changed_y);CHKERRQ(ierr); 1259 } 1260 1261 if (PetscLogPrintInfo){ 1262 ierr = SNESVICheckResidual_Private(snes,snes->jacobian,F,Y,G,W);CHKERRQ(ierr); 1263 } 1264 */ 1265 /* Compute a (scaled) negative update in the line search routine: 1266 Y <- X - lambda*Y 1267 and evaluate G = function(Y) (depends on the line search). 1268 */ 1269 ierr = VecCopy(Y,snes->vec_sol_update);CHKERRQ(ierr); 1270 ynorm = 1; gnorm = fnorm; 1271 ierr = (*vi->LineSearch)(snes,vi->lsP,X,F,G,Y,W,fnorm,xnorm,&ynorm,&gnorm,&lssucceed);CHKERRQ(ierr); 1272 ierr = PetscInfo4(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lssucceed=%d\n",fnorm,gnorm,ynorm,(int)lssucceed);CHKERRQ(ierr); 1273 if (snes->reason == SNES_DIVERGED_FUNCTION_COUNT) break; 1274 if (snes->domainerror) { 1275 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 1276 ierr = DMDestroyVI(snes->dm);CHKERRQ(ierr); 1277 PetscFunctionReturn(0); 1278 } 1279 if (!lssucceed) { 1280 if (++snes->numFailures >= snes->maxFailures) { 1281 PetscBool ismin; 1282 snes->reason = SNES_DIVERGED_LINE_SEARCH; 1283 ierr = SNESVICheckLocalMin_Private(snes,snes->jacobian,G,W,gnorm,&ismin);CHKERRQ(ierr); 1284 if (ismin) snes->reason = SNES_DIVERGED_LOCAL_MIN; 1285 break; 1286 } 1287 } 1288 /* Update function and solution vectors */ 1289 fnorm = gnorm; 1290 ierr = VecCopy(G,F);CHKERRQ(ierr); 1291 ierr = VecCopy(W,X);CHKERRQ(ierr); 1292 /* Monitor convergence */ 1293 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 1294 snes->iter = i+1; 1295 snes->norm = fnorm; 1296 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 1297 SNESLogConvHistory(snes,snes->norm,lits); 1298 ierr = SNESMonitor(snes,snes->iter,snes->norm);CHKERRQ(ierr); 1299 /* Test for convergence, xnorm = || X || */ 1300 if (snes->ops->converged != SNESSkipConverged) { ierr = VecNorm(X,NORM_2,&xnorm);CHKERRQ(ierr); } 1301 ierr = (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 1302 if (snes->reason) break; 1303 } 1304 ierr = DMDestroyVI(snes->dm);CHKERRQ(ierr); 1305 if (i == maxits) { 1306 ierr = PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",maxits);CHKERRQ(ierr); 1307 if(!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT; 1308 } 1309 PetscFunctionReturn(0); 1310 } 1311 1312 #undef __FUNCT__ 1313 #define __FUNCT__ "SNESVISetRedundancyCheck" 1314 PetscErrorCode SNESVISetRedundancyCheck(SNES snes,PetscErrorCode (*func)(SNES,IS,IS*,void*),void *ctx) 1315 { 1316 SNES_VI *vi = (SNES_VI*)snes->data; 1317 1318 PetscFunctionBegin; 1319 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1320 vi->checkredundancy = func; 1321 vi->ctxP = ctx; 1322 PetscFunctionReturn(0); 1323 } 1324 1325 #if defined(PETSC_HAVE_MATLAB_ENGINE) 1326 #include <engine.h> 1327 #include <mex.h> 1328 typedef struct {char *funcname; mxArray *ctx;} SNESMatlabContext; 1329 1330 #undef __FUNCT__ 1331 #define __FUNCT__ "SNESVIRedundancyCheck_Matlab" 1332 PetscErrorCode SNESVIRedundancyCheck_Matlab(SNES snes,IS is_act,IS* is_redact,void* ctx) 1333 { 1334 PetscErrorCode ierr; 1335 SNESMatlabContext *sctx = (SNESMatlabContext*)ctx; 1336 int nlhs = 1, nrhs = 5; 1337 mxArray *plhs[1], *prhs[5]; 1338 long long int l1 = 0, l2 = 0, ls = 0; 1339 PetscInt *indices=PETSC_NULL; 1340 1341 PetscFunctionBegin; 1342 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1343 PetscValidHeaderSpecific(is_act,IS_CLASSID,2); 1344 PetscValidPointer(is_redact,3); 1345 PetscCheckSameComm(snes,1,is_act,2); 1346 1347 /* Create IS for reduced active set of size 0, its size and indices will 1348 bet set by the Matlab function */ 1349 ierr = ISCreateGeneral(((PetscObject)snes)->comm,0,indices,PETSC_OWN_POINTER,is_redact);CHKERRQ(ierr); 1350 /* call Matlab function in ctx */ 1351 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 1352 ierr = PetscMemcpy(&l1,&is_act,sizeof(is_act));CHKERRQ(ierr); 1353 ierr = PetscMemcpy(&l2,is_redact,sizeof(is_act));CHKERRQ(ierr); 1354 prhs[0] = mxCreateDoubleScalar((double)ls); 1355 prhs[1] = mxCreateDoubleScalar((double)l1); 1356 prhs[2] = mxCreateDoubleScalar((double)l2); 1357 prhs[3] = mxCreateString(sctx->funcname); 1358 prhs[4] = sctx->ctx; 1359 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESVIRedundancyCheckInternal");CHKERRQ(ierr); 1360 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 1361 mxDestroyArray(prhs[0]); 1362 mxDestroyArray(prhs[1]); 1363 mxDestroyArray(prhs[2]); 1364 mxDestroyArray(prhs[3]); 1365 mxDestroyArray(plhs[0]); 1366 PetscFunctionReturn(0); 1367 } 1368 1369 #undef __FUNCT__ 1370 #define __FUNCT__ "SNESVISetRedundancyCheckMatlab" 1371 PetscErrorCode SNESVISetRedundancyCheckMatlab(SNES snes,const char* func,mxArray* ctx) 1372 { 1373 PetscErrorCode ierr; 1374 SNESMatlabContext *sctx; 1375 1376 PetscFunctionBegin; 1377 /* currently sctx is memory bleed */ 1378 ierr = PetscMalloc(sizeof(SNESMatlabContext),&sctx);CHKERRQ(ierr); 1379 ierr = PetscStrallocpy(func,&sctx->funcname);CHKERRQ(ierr); 1380 sctx->ctx = mxDuplicateArray(ctx); 1381 ierr = SNESVISetRedundancyCheck(snes,SNESVIRedundancyCheck_Matlab,sctx);CHKERRQ(ierr); 1382 PetscFunctionReturn(0); 1383 } 1384 1385 #endif 1386 1387 1388 /* Variational Inequality solver using augmented space method. It does the opposite of the 1389 reduced space method i.e. it identifies the active set variables and instead of discarding 1390 them it augments the original system by introducing additional equality 1391 constraint equations for active set variables. The user can optionally provide an IS for 1392 a subset of 'redundant' active set variables which will treated as free variables. 1393 Specific implementation for Allen-Cahn problem 1394 */ 1395 #undef __FUNCT__ 1396 #define __FUNCT__ "SNESSolveVI_RSAUG" 1397 PetscErrorCode SNESSolveVI_RSAUG(SNES snes) 1398 { 1399 SNES_VI *vi = (SNES_VI*)snes->data; 1400 PetscErrorCode ierr; 1401 PetscInt maxits,i,lits; 1402 PetscBool lssucceed; 1403 MatStructure flg = DIFFERENT_NONZERO_PATTERN; 1404 PetscReal fnorm,gnorm,xnorm=0,ynorm; 1405 Vec Y,X,F,G,W; 1406 KSPConvergedReason kspreason; 1407 1408 PetscFunctionBegin; 1409 snes->numFailures = 0; 1410 snes->numLinearSolveFailures = 0; 1411 snes->reason = SNES_CONVERGED_ITERATING; 1412 1413 maxits = snes->max_its; /* maximum number of iterations */ 1414 X = snes->vec_sol; /* solution vector */ 1415 F = snes->vec_func; /* residual vector */ 1416 Y = snes->work[0]; /* work vectors */ 1417 G = snes->work[1]; 1418 W = snes->work[2]; 1419 1420 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 1421 snes->iter = 0; 1422 snes->norm = 0.0; 1423 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 1424 1425 ierr = SNESVIProjectOntoBounds(snes,X);CHKERRQ(ierr); 1426 ierr = SNESComputeFunction(snes,X,F);CHKERRQ(ierr); 1427 if (snes->domainerror) { 1428 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 1429 PetscFunctionReturn(0); 1430 } 1431 ierr = SNESVIComputeInactiveSetFnorm(snes,F,X,&fnorm);CHKERRQ(ierr); 1432 ierr = VecNormBegin(X,NORM_2,&xnorm);CHKERRQ(ierr); /* xnorm <- ||x|| */ 1433 ierr = VecNormEnd(X,NORM_2,&xnorm);CHKERRQ(ierr); 1434 if (PetscIsInfOrNanReal(fnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1435 1436 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 1437 snes->norm = fnorm; 1438 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 1439 SNESLogConvHistory(snes,fnorm,0); 1440 ierr = SNESMonitor(snes,0,fnorm);CHKERRQ(ierr); 1441 1442 /* set parameter for default relative tolerance convergence test */ 1443 snes->ttol = fnorm*snes->rtol; 1444 /* test convergence */ 1445 ierr = (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 1446 if (snes->reason) PetscFunctionReturn(0); 1447 1448 for (i=0; i<maxits; i++) { 1449 IS IS_act,IS_inact; /* _act -> active set _inact -> inactive set */ 1450 IS IS_redact; /* redundant active set */ 1451 Mat J_aug,Jpre_aug; 1452 Vec F_aug,Y_aug; 1453 PetscInt nis_redact,nis_act; 1454 const PetscInt *idx_redact,*idx_act; 1455 PetscInt k; 1456 PetscInt *idx_actkept=PETSC_NULL,nkept=0; /* list of kept active set */ 1457 PetscScalar *f,*f2; 1458 PetscBool isequal; 1459 PetscInt i1=0,j1=0; 1460 1461 /* Call general purpose update function */ 1462 if (snes->ops->update) { 1463 ierr = (*snes->ops->update)(snes, snes->iter);CHKERRQ(ierr); 1464 } 1465 ierr = SNESComputeJacobian(snes,X,&snes->jacobian,&snes->jacobian_pre,&flg);CHKERRQ(ierr); 1466 1467 /* Create active and inactive index sets */ 1468 ierr = SNESVICreateIndexSets_RS(snes,X,F,&IS_act,&IS_inact);CHKERRQ(ierr); 1469 1470 /* Get local active set size */ 1471 ierr = ISGetLocalSize(IS_act,&nis_act);CHKERRQ(ierr); 1472 if (nis_act) { 1473 ierr = ISGetIndices(IS_act,&idx_act);CHKERRQ(ierr); 1474 IS_redact = PETSC_NULL; 1475 if(vi->checkredundancy) { 1476 (*vi->checkredundancy)(snes,IS_act,&IS_redact,vi->ctxP); 1477 } 1478 1479 if(!IS_redact) { 1480 /* User called checkredundancy function but didn't create IS_redact because 1481 there were no redundant active set variables */ 1482 /* Copy over all active set indices to the list */ 1483 ierr = PetscMalloc(nis_act*sizeof(PetscInt),&idx_actkept);CHKERRQ(ierr); 1484 for(k=0;k < nis_act;k++) idx_actkept[k] = idx_act[k]; 1485 nkept = nis_act; 1486 } else { 1487 ierr = ISGetLocalSize(IS_redact,&nis_redact);CHKERRQ(ierr); 1488 ierr = PetscMalloc((nis_act-nis_redact)*sizeof(PetscInt),&idx_actkept);CHKERRQ(ierr); 1489 1490 /* Create reduced active set list */ 1491 ierr = ISGetIndices(IS_act,&idx_act);CHKERRQ(ierr); 1492 ierr = ISGetIndices(IS_redact,&idx_redact);CHKERRQ(ierr); 1493 j1 = 0;nkept = 0; 1494 for(k=0;k<nis_act;k++) { 1495 if(j1 < nis_redact && idx_act[k] == idx_redact[j1]) j1++; 1496 else idx_actkept[nkept++] = idx_act[k]; 1497 } 1498 ierr = ISRestoreIndices(IS_act,&idx_act);CHKERRQ(ierr); 1499 ierr = ISRestoreIndices(IS_redact,&idx_redact);CHKERRQ(ierr); 1500 1501 ierr = ISDestroy(&IS_redact);CHKERRQ(ierr); 1502 } 1503 1504 /* Create augmented F and Y */ 1505 ierr = VecCreate(((PetscObject)snes)->comm,&F_aug);CHKERRQ(ierr); 1506 ierr = VecSetSizes(F_aug,F->map->n+nkept,PETSC_DECIDE);CHKERRQ(ierr); 1507 ierr = VecSetFromOptions(F_aug);CHKERRQ(ierr); 1508 ierr = VecDuplicate(F_aug,&Y_aug);CHKERRQ(ierr); 1509 1510 /* Copy over F to F_aug in the first n locations */ 1511 ierr = VecGetArray(F,&f);CHKERRQ(ierr); 1512 ierr = VecGetArray(F_aug,&f2);CHKERRQ(ierr); 1513 ierr = PetscMemcpy(f2,f,F->map->n*sizeof(PetscScalar));CHKERRQ(ierr); 1514 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr); 1515 ierr = VecRestoreArray(F_aug,&f2);CHKERRQ(ierr); 1516 1517 /* Create the augmented jacobian matrix */ 1518 ierr = MatCreate(((PetscObject)snes)->comm,&J_aug);CHKERRQ(ierr); 1519 ierr = MatSetSizes(J_aug,X->map->n+nkept,X->map->n+nkept,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr); 1520 ierr = MatSetFromOptions(J_aug);CHKERRQ(ierr); 1521 1522 1523 { /* local vars */ 1524 /* Preallocate augmented matrix and set values in it...Doing only sequential case first*/ 1525 PetscInt ncols; 1526 const PetscInt *cols; 1527 const PetscScalar *vals; 1528 PetscScalar value[2]; 1529 PetscInt row,col[2]; 1530 PetscInt *d_nnz; 1531 value[0] = 1.0; value[1] = 0.0; 1532 ierr = PetscMalloc((X->map->n+nkept)*sizeof(PetscInt),&d_nnz);CHKERRQ(ierr); 1533 ierr = PetscMemzero(d_nnz,(X->map->n+nkept)*sizeof(PetscInt));CHKERRQ(ierr); 1534 for(row=0;row<snes->jacobian->rmap->n;row++) { 1535 ierr = MatGetRow(snes->jacobian,row,&ncols,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 1536 d_nnz[row] += ncols; 1537 ierr = MatRestoreRow(snes->jacobian,row,&ncols,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 1538 } 1539 1540 for(k=0;k<nkept;k++) { 1541 d_nnz[idx_actkept[k]] += 1; 1542 d_nnz[snes->jacobian->rmap->n+k] += 2; 1543 } 1544 ierr = MatSeqAIJSetPreallocation(J_aug,PETSC_NULL,d_nnz);CHKERRQ(ierr); 1545 1546 ierr = PetscFree(d_nnz);CHKERRQ(ierr); 1547 1548 /* Set the original jacobian matrix in J_aug */ 1549 for(row=0;row<snes->jacobian->rmap->n;row++) { 1550 ierr = MatGetRow(snes->jacobian,row,&ncols,&cols,&vals);CHKERRQ(ierr); 1551 ierr = MatSetValues(J_aug,1,&row,ncols,cols,vals,INSERT_VALUES);CHKERRQ(ierr); 1552 ierr = MatRestoreRow(snes->jacobian,row,&ncols,&cols,&vals);CHKERRQ(ierr); 1553 } 1554 /* Add the augmented part */ 1555 for(k=0;k<nkept;k++) { 1556 row = snes->jacobian->rmap->n+k; 1557 col[0] = idx_actkept[k]; col[1] = snes->jacobian->rmap->n+k; 1558 ierr = MatSetValues(J_aug,1,&row,1,col,value,INSERT_VALUES);CHKERRQ(ierr); 1559 ierr = MatSetValues(J_aug,1,&col[0],1,&row,&value[0],INSERT_VALUES);CHKERRQ(ierr); 1560 } 1561 ierr = MatAssemblyBegin(J_aug,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1562 ierr = MatAssemblyEnd(J_aug,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1563 /* Only considering prejac = jac for now */ 1564 Jpre_aug = J_aug; 1565 } /* local vars*/ 1566 ierr = ISRestoreIndices(IS_act,&idx_act);CHKERRQ(ierr); 1567 ierr = ISDestroy(&IS_act);CHKERRQ(ierr); 1568 } else { 1569 F_aug = F; J_aug = snes->jacobian; Y_aug = Y; Jpre_aug = snes->jacobian_pre; 1570 } 1571 1572 ierr = ISEqual(vi->IS_inact_prev,IS_inact,&isequal);CHKERRQ(ierr); 1573 if (!isequal) { 1574 ierr = SNESVIResetPCandKSP(snes,J_aug,Jpre_aug);CHKERRQ(ierr); 1575 } 1576 ierr = KSPSetOperators(snes->ksp,J_aug,Jpre_aug,flg);CHKERRQ(ierr); 1577 ierr = KSPSetUp(snes->ksp);CHKERRQ(ierr); 1578 /* { 1579 PC pc; 1580 PetscBool flg; 1581 ierr = KSPGetPC(snes->ksp,&pc);CHKERRQ(ierr); 1582 ierr = PetscTypeCompare((PetscObject)pc,PCFIELDSPLIT,&flg);CHKERRQ(ierr); 1583 if (flg) { 1584 KSP *subksps; 1585 ierr = PCFieldSplitGetSubKSP(pc,PETSC_NULL,&subksps);CHKERRQ(ierr); 1586 ierr = KSPGetPC(subksps[0],&pc);CHKERRQ(ierr); 1587 ierr = PetscFree(subksps);CHKERRQ(ierr); 1588 ierr = PetscTypeCompare((PetscObject)pc,PCBJACOBI,&flg);CHKERRQ(ierr); 1589 if (flg) { 1590 PetscInt n,N = 101*101,j,cnts[3] = {0,0,0}; 1591 const PetscInt *ii; 1592 1593 ierr = ISGetSize(IS_inact,&n);CHKERRQ(ierr); 1594 ierr = ISGetIndices(IS_inact,&ii);CHKERRQ(ierr); 1595 for (j=0; j<n; j++) { 1596 if (ii[j] < N) cnts[0]++; 1597 else if (ii[j] < 2*N) cnts[1]++; 1598 else if (ii[j] < 3*N) cnts[2]++; 1599 } 1600 ierr = ISRestoreIndices(IS_inact,&ii);CHKERRQ(ierr); 1601 1602 ierr = PCBJacobiSetTotalBlocks(pc,3,cnts);CHKERRQ(ierr); 1603 } 1604 } 1605 } 1606 */ 1607 ierr = SNES_KSPSolve(snes,snes->ksp,F_aug,Y_aug);CHKERRQ(ierr); 1608 ierr = KSPGetConvergedReason(snes->ksp,&kspreason);CHKERRQ(ierr); 1609 if (kspreason < 0) { 1610 if (++snes->numLinearSolveFailures >= snes->maxLinearSolveFailures) { 1611 ierr = PetscInfo2(snes,"iter=%D, number linear solve failures %D greater than current SNES allowed, stopping solve\n",snes->iter,snes->numLinearSolveFailures);CHKERRQ(ierr); 1612 snes->reason = SNES_DIVERGED_LINEAR_SOLVE; 1613 break; 1614 } 1615 } 1616 1617 if(nis_act) { 1618 PetscScalar *y1,*y2; 1619 ierr = VecGetArray(Y,&y1);CHKERRQ(ierr); 1620 ierr = VecGetArray(Y_aug,&y2);CHKERRQ(ierr); 1621 /* Copy over inactive Y_aug to Y */ 1622 j1 = 0; 1623 for(i1=Y->map->rstart;i1 < Y->map->rend;i1++) { 1624 if(j1 < nkept && idx_actkept[j1] == i1) j1++; 1625 else y1[i1-Y->map->rstart] = y2[i1-Y->map->rstart]; 1626 } 1627 ierr = VecRestoreArray(Y,&y1);CHKERRQ(ierr); 1628 ierr = VecRestoreArray(Y_aug,&y2);CHKERRQ(ierr); 1629 1630 ierr = VecDestroy(&F_aug);CHKERRQ(ierr); 1631 ierr = VecDestroy(&Y_aug);CHKERRQ(ierr); 1632 ierr = MatDestroy(&J_aug);CHKERRQ(ierr); 1633 ierr = PetscFree(idx_actkept);CHKERRQ(ierr); 1634 } 1635 1636 if (!isequal) { 1637 ierr = ISDestroy(&vi->IS_inact_prev);CHKERRQ(ierr); 1638 ierr = ISDuplicate(IS_inact,&vi->IS_inact_prev);CHKERRQ(ierr); 1639 } 1640 ierr = ISDestroy(&IS_inact);CHKERRQ(ierr); 1641 1642 ierr = KSPGetIterationNumber(snes->ksp,&lits);CHKERRQ(ierr); 1643 snes->linear_its += lits; 1644 ierr = PetscInfo2(snes,"iter=%D, linear solve iterations=%D\n",snes->iter,lits);CHKERRQ(ierr); 1645 /* 1646 if (vi->precheckstep) { 1647 PetscBool changed_y = PETSC_FALSE; 1648 ierr = (*vi->precheckstep)(snes,X,Y,vi->precheck,&changed_y);CHKERRQ(ierr); 1649 } 1650 1651 if (PetscLogPrintInfo){ 1652 ierr = SNESVICheckResidual_Private(snes,snes->jacobian,F,Y,G,W);CHKERRQ(ierr); 1653 } 1654 */ 1655 /* Compute a (scaled) negative update in the line search routine: 1656 Y <- X - lambda*Y 1657 and evaluate G = function(Y) (depends on the line search). 1658 */ 1659 ierr = VecCopy(Y,snes->vec_sol_update);CHKERRQ(ierr); 1660 ynorm = 1; gnorm = fnorm; 1661 ierr = (*vi->LineSearch)(snes,vi->lsP,X,F,G,Y,W,fnorm,xnorm,&ynorm,&gnorm,&lssucceed);CHKERRQ(ierr); 1662 ierr = PetscInfo4(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lssucceed=%d\n",fnorm,gnorm,ynorm,(int)lssucceed);CHKERRQ(ierr); 1663 if (snes->reason == SNES_DIVERGED_FUNCTION_COUNT) break; 1664 if (snes->domainerror) { 1665 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 1666 PetscFunctionReturn(0); 1667 } 1668 if (!lssucceed) { 1669 if (++snes->numFailures >= snes->maxFailures) { 1670 PetscBool ismin; 1671 snes->reason = SNES_DIVERGED_LINE_SEARCH; 1672 ierr = SNESVICheckLocalMin_Private(snes,snes->jacobian,G,W,gnorm,&ismin);CHKERRQ(ierr); 1673 if (ismin) snes->reason = SNES_DIVERGED_LOCAL_MIN; 1674 break; 1675 } 1676 } 1677 /* Update function and solution vectors */ 1678 fnorm = gnorm; 1679 ierr = VecCopy(G,F);CHKERRQ(ierr); 1680 ierr = VecCopy(W,X);CHKERRQ(ierr); 1681 /* Monitor convergence */ 1682 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 1683 snes->iter = i+1; 1684 snes->norm = fnorm; 1685 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 1686 SNESLogConvHistory(snes,snes->norm,lits); 1687 ierr = SNESMonitor(snes,snes->iter,snes->norm);CHKERRQ(ierr); 1688 /* Test for convergence, xnorm = || X || */ 1689 if (snes->ops->converged != SNESSkipConverged) { ierr = VecNorm(X,NORM_2,&xnorm);CHKERRQ(ierr); } 1690 ierr = (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 1691 if (snes->reason) break; 1692 } 1693 if (i == maxits) { 1694 ierr = PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",maxits);CHKERRQ(ierr); 1695 if(!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT; 1696 } 1697 PetscFunctionReturn(0); 1698 } 1699 1700 /* -------------------------------------------------------------------------- */ 1701 /* 1702 SNESSetUp_VI - Sets up the internal data structures for the later use 1703 of the SNESVI nonlinear solver. 1704 1705 Input Parameter: 1706 . snes - the SNES context 1707 . x - the solution vector 1708 1709 Application Interface Routine: SNESSetUp() 1710 1711 Notes: 1712 For basic use of the SNES solvers, the user need not explicitly call 1713 SNESSetUp(), since these actions will automatically occur during 1714 the call to SNESSolve(). 1715 */ 1716 #undef __FUNCT__ 1717 #define __FUNCT__ "SNESSetUp_VI" 1718 PetscErrorCode SNESSetUp_VI(SNES snes) 1719 { 1720 PetscErrorCode ierr; 1721 SNES_VI *vi = (SNES_VI*) snes->data; 1722 PetscInt i_start[3],i_end[3]; 1723 1724 PetscFunctionBegin; 1725 1726 ierr = SNESDefaultGetWork(snes,3);CHKERRQ(ierr); 1727 1728 if (vi->computevariablebounds) { 1729 ierr = VecDestroy(&vi->xl);CHKERRQ(ierr); 1730 ierr = VecDestroy(&vi->xu);CHKERRQ(ierr); 1731 ierr = (*vi->computevariablebounds)(snes,&vi->xl,&vi->xu);CHKERRQ(ierr); 1732 } else if (!vi->xl && !vi->xu) { 1733 /* If the lower and upper bound on variables are not set, set it to -Inf and Inf */ 1734 ierr = VecDuplicate(snes->vec_sol, &vi->xl);CHKERRQ(ierr); 1735 ierr = VecSet(vi->xl,SNES_VI_NINF);CHKERRQ(ierr); 1736 ierr = VecDuplicate(snes->vec_sol, &vi->xu);CHKERRQ(ierr); 1737 ierr = VecSet(vi->xu,SNES_VI_INF);CHKERRQ(ierr); 1738 } else { 1739 /* Check if lower bound, upper bound and solution vector distribution across the processors is identical */ 1740 ierr = VecGetOwnershipRange(snes->vec_sol,i_start,i_end);CHKERRQ(ierr); 1741 ierr = VecGetOwnershipRange(vi->xl,i_start+1,i_end+1);CHKERRQ(ierr); 1742 ierr = VecGetOwnershipRange(vi->xu,i_start+2,i_end+2);CHKERRQ(ierr); 1743 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])) 1744 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."); 1745 } 1746 if (snes->ops->solve != SNESSolveVI_SS) { 1747 /* Set up previous active index set for the first snes solve 1748 vi->IS_inact_prev = 0,1,2,....N */ 1749 PetscInt *indices; 1750 PetscInt i,n,rstart,rend; 1751 1752 ierr = VecGetOwnershipRange(snes->vec_sol,&rstart,&rend);CHKERRQ(ierr); 1753 ierr = VecGetLocalSize(snes->vec_sol,&n);CHKERRQ(ierr); 1754 ierr = PetscMalloc(n*sizeof(PetscInt),&indices);CHKERRQ(ierr); 1755 for(i=0;i < n; i++) indices[i] = rstart + i; 1756 ierr = ISCreateGeneral(((PetscObject)snes)->comm,n,indices,PETSC_OWN_POINTER,&vi->IS_inact_prev); 1757 } 1758 1759 if (snes->ops->solve == SNESSolveVI_SS) { 1760 ierr = VecDuplicate(snes->vec_sol, &vi->dpsi);CHKERRQ(ierr); 1761 ierr = VecDuplicate(snes->vec_sol, &vi->phi);CHKERRQ(ierr); 1762 ierr = VecDuplicate(snes->vec_sol, &vi->Da);CHKERRQ(ierr); 1763 ierr = VecDuplicate(snes->vec_sol, &vi->Db);CHKERRQ(ierr); 1764 ierr = VecDuplicate(snes->vec_sol, &vi->z);CHKERRQ(ierr); 1765 ierr = VecDuplicate(snes->vec_sol, &vi->t);CHKERRQ(ierr); 1766 } 1767 PetscFunctionReturn(0); 1768 } 1769 /* -------------------------------------------------------------------------- */ 1770 #undef __FUNCT__ 1771 #define __FUNCT__ "SNESReset_VI" 1772 PetscErrorCode SNESReset_VI(SNES snes) 1773 { 1774 SNES_VI *vi = (SNES_VI*) snes->data; 1775 PetscErrorCode ierr; 1776 1777 PetscFunctionBegin; 1778 ierr = VecDestroy(&vi->xl);CHKERRQ(ierr); 1779 ierr = VecDestroy(&vi->xu);CHKERRQ(ierr); 1780 if (snes->ops->solve != SNESSolveVI_SS) { 1781 ierr = ISDestroy(&vi->IS_inact_prev);CHKERRQ(ierr); 1782 } 1783 PetscFunctionReturn(0); 1784 } 1785 1786 /* 1787 SNESDestroy_VI - Destroys the private SNES_VI context that was created 1788 with SNESCreate_VI(). 1789 1790 Input Parameter: 1791 . snes - the SNES context 1792 1793 Application Interface Routine: SNESDestroy() 1794 */ 1795 #undef __FUNCT__ 1796 #define __FUNCT__ "SNESDestroy_VI" 1797 PetscErrorCode SNESDestroy_VI(SNES snes) 1798 { 1799 SNES_VI *vi = (SNES_VI*) snes->data; 1800 PetscErrorCode ierr; 1801 1802 PetscFunctionBegin; 1803 1804 if (snes->ops->solve == SNESSolveVI_SS) { 1805 /* clear vectors */ 1806 ierr = VecDestroy(&vi->dpsi);CHKERRQ(ierr); 1807 ierr = VecDestroy(&vi->phi);CHKERRQ(ierr); 1808 ierr = VecDestroy(&vi->Da);CHKERRQ(ierr); 1809 ierr = VecDestroy(&vi->Db);CHKERRQ(ierr); 1810 ierr = VecDestroy(&vi->z);CHKERRQ(ierr); 1811 ierr = VecDestroy(&vi->t);CHKERRQ(ierr); 1812 } 1813 1814 ierr = PetscViewerASCIIMonitorDestroy(&vi->lsmonitor);CHKERRQ(ierr); 1815 ierr = PetscFree(snes->data);CHKERRQ(ierr); 1816 1817 /* clear composed functions */ 1818 ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESLineSearchSet_C","",PETSC_NULL);CHKERRQ(ierr); 1819 ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESLineSearchSetMonitor_C","",PETSC_NULL);CHKERRQ(ierr); 1820 PetscFunctionReturn(0); 1821 } 1822 1823 /* -------------------------------------------------------------------------- */ 1824 #undef __FUNCT__ 1825 #define __FUNCT__ "SNESLineSearchNo_VI" 1826 1827 /* 1828 This routine does not actually do a line search but it takes a full newton 1829 step while ensuring that the new iterates remain within the constraints. 1830 1831 */ 1832 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) 1833 { 1834 PetscErrorCode ierr; 1835 SNES_VI *vi = (SNES_VI*)snes->data; 1836 PetscBool changed_w = PETSC_FALSE,changed_y = PETSC_FALSE; 1837 1838 PetscFunctionBegin; 1839 *flag = PETSC_TRUE; 1840 ierr = PetscLogEventBegin(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1841 ierr = VecNorm(y,NORM_2,ynorm);CHKERRQ(ierr); /* ynorm = || y || */ 1842 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); /* w <- x - y */ 1843 ierr = SNESVIProjectOntoBounds(snes,w);CHKERRQ(ierr); 1844 if (vi->postcheckstep) { 1845 ierr = (*vi->postcheckstep)(snes,x,y,w,vi->postcheck,&changed_y,&changed_w);CHKERRQ(ierr); 1846 } 1847 if (changed_y) { 1848 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); /* w <- x - y */ 1849 ierr = SNESVIProjectOntoBounds(snes,w);CHKERRQ(ierr); 1850 } 1851 ierr = SNESVIProjectOntoBounds(snes,w);CHKERRQ(ierr); 1852 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1853 if (!snes->domainerror) { 1854 if (snes->ops->solve != SNESSolveVI_SS) { 1855 ierr = SNESVIComputeInactiveSetFnorm(snes,g,w,gnorm);CHKERRQ(ierr); 1856 } else { 1857 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); /* gnorm = || g || */ 1858 } 1859 if (PetscIsInfOrNanReal(*gnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1860 } 1861 if (vi->lsmonitor) { 1862 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: Using full step: fnorm %G gnorm %G\n",fnorm,*gnorm);CHKERRQ(ierr); 1863 } 1864 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1865 PetscFunctionReturn(0); 1866 } 1867 1868 /* -------------------------------------------------------------------------- */ 1869 #undef __FUNCT__ 1870 #define __FUNCT__ "SNESLineSearchNoNorms_VI" 1871 1872 /* 1873 This routine is a copy of SNESLineSearchNoNorms in snes/impls/ls/ls.c 1874 */ 1875 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) 1876 { 1877 PetscErrorCode ierr; 1878 SNES_VI *vi = (SNES_VI*)snes->data; 1879 PetscBool changed_w = PETSC_FALSE,changed_y = PETSC_FALSE; 1880 1881 PetscFunctionBegin; 1882 *flag = PETSC_TRUE; 1883 ierr = PetscLogEventBegin(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1884 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); /* w <- x - y */ 1885 ierr = SNESVIProjectOntoBounds(snes,w);CHKERRQ(ierr); 1886 if (vi->postcheckstep) { 1887 ierr = (*vi->postcheckstep)(snes,x,y,w,vi->postcheck,&changed_y,&changed_w);CHKERRQ(ierr); 1888 } 1889 if (changed_y) { 1890 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); /* w <- x - y */ 1891 ierr = SNESVIProjectOntoBounds(snes,w);CHKERRQ(ierr); 1892 } 1893 1894 /* don't evaluate function the last time through */ 1895 if (snes->iter < snes->max_its-1) { 1896 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1897 } 1898 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1899 PetscFunctionReturn(0); 1900 } 1901 1902 /* -------------------------------------------------------------------------- */ 1903 #undef __FUNCT__ 1904 #define __FUNCT__ "SNESLineSearchCubic_VI" 1905 /* 1906 This routine implements a cubic line search while doing a projection on the variable bounds 1907 */ 1908 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) 1909 { 1910 PetscReal initslope,lambdaprev,gnormprev,a,b,d,t1,t2,rellength; 1911 PetscReal minlambda,lambda,lambdatemp; 1912 #if defined(PETSC_USE_COMPLEX) 1913 PetscScalar cinitslope; 1914 #endif 1915 PetscErrorCode ierr; 1916 PetscInt count; 1917 SNES_VI *vi = (SNES_VI*)snes->data; 1918 PetscBool changed_w = PETSC_FALSE,changed_y = PETSC_FALSE; 1919 MPI_Comm comm; 1920 1921 PetscFunctionBegin; 1922 ierr = PetscObjectGetComm((PetscObject)snes,&comm);CHKERRQ(ierr); 1923 ierr = PetscLogEventBegin(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1924 *flag = PETSC_TRUE; 1925 1926 ierr = VecNorm(y,NORM_2,ynorm);CHKERRQ(ierr); 1927 if (*ynorm == 0.0) { 1928 if (vi->lsmonitor) { 1929 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: Initial direction and size is 0\n");CHKERRQ(ierr); 1930 } 1931 *gnorm = fnorm; 1932 ierr = VecCopy(x,w);CHKERRQ(ierr); 1933 ierr = VecCopy(f,g);CHKERRQ(ierr); 1934 *flag = PETSC_FALSE; 1935 goto theend1; 1936 } 1937 if (*ynorm > vi->maxstep) { /* Step too big, so scale back */ 1938 if (vi->lsmonitor) { 1939 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: Scaling step by %G old ynorm %G\n",vi->maxstep/(*ynorm),*ynorm);CHKERRQ(ierr); 1940 } 1941 ierr = VecScale(y,vi->maxstep/(*ynorm));CHKERRQ(ierr); 1942 *ynorm = vi->maxstep; 1943 } 1944 ierr = VecMaxPointwiseDivide(y,x,&rellength);CHKERRQ(ierr); 1945 minlambda = vi->minlambda/rellength; 1946 ierr = MatMult(snes->jacobian,y,w);CHKERRQ(ierr); 1947 #if defined(PETSC_USE_COMPLEX) 1948 ierr = VecDot(f,w,&cinitslope);CHKERRQ(ierr); 1949 initslope = PetscRealPart(cinitslope); 1950 #else 1951 ierr = VecDot(f,w,&initslope);CHKERRQ(ierr); 1952 #endif 1953 if (initslope > 0.0) initslope = -initslope; 1954 if (initslope == 0.0) initslope = -1.0; 1955 1956 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); 1957 ierr = SNESVIProjectOntoBounds(snes,w);CHKERRQ(ierr); 1958 if (snes->nfuncs >= snes->max_funcs) { 1959 ierr = PetscInfo(snes,"Exceeded maximum function evaluations, while checking full step length!\n");CHKERRQ(ierr); 1960 *flag = PETSC_FALSE; 1961 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 1962 goto theend1; 1963 } 1964 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1965 if (snes->ops->solve != SNESSolveVI_SS) { 1966 ierr = SNESVIComputeInactiveSetFnorm(snes,g,w,gnorm);CHKERRQ(ierr); 1967 } else { 1968 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); 1969 } 1970 if (snes->domainerror) { 1971 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1972 PetscFunctionReturn(0); 1973 } 1974 if (PetscIsInfOrNanReal(*gnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1975 ierr = PetscInfo2(snes,"Initial fnorm %G gnorm %G\n",fnorm,*gnorm);CHKERRQ(ierr); 1976 if (.5*(*gnorm)*(*gnorm) <= .5*fnorm*fnorm + vi->alpha*initslope) { /* Sufficient reduction */ 1977 if (vi->lsmonitor) { 1978 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: Using full step: fnorm %G gnorm %G\n",fnorm,*gnorm);CHKERRQ(ierr); 1979 } 1980 goto theend1; 1981 } 1982 1983 /* Fit points with quadratic */ 1984 lambda = 1.0; 1985 lambdatemp = -initslope/((*gnorm)*(*gnorm) - fnorm*fnorm - 2.0*initslope); 1986 lambdaprev = lambda; 1987 gnormprev = *gnorm; 1988 if (lambdatemp > .5*lambda) lambdatemp = .5*lambda; 1989 if (lambdatemp <= .1*lambda) lambda = .1*lambda; 1990 else lambda = lambdatemp; 1991 1992 ierr = VecWAXPY(w,-lambda,y,x);CHKERRQ(ierr); 1993 ierr = SNESVIProjectOntoBounds(snes,w);CHKERRQ(ierr); 1994 if (snes->nfuncs >= snes->max_funcs) { 1995 ierr = PetscInfo1(snes,"Exceeded maximum function evaluations, while attempting quadratic backtracking! %D \n",snes->nfuncs);CHKERRQ(ierr); 1996 *flag = PETSC_FALSE; 1997 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 1998 goto theend1; 1999 } 2000 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 2001 if (snes->ops->solve != SNESSolveVI_SS) { 2002 ierr = SNESVIComputeInactiveSetFnorm(snes,g,w,gnorm);CHKERRQ(ierr); 2003 } else { 2004 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); 2005 } 2006 if (snes->domainerror) { 2007 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 2008 PetscFunctionReturn(0); 2009 } 2010 if (PetscIsInfOrNanReal(*gnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 2011 if (vi->lsmonitor) { 2012 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: gnorm after quadratic fit %G\n",*gnorm);CHKERRQ(ierr); 2013 } 2014 if (.5*(*gnorm)*(*gnorm) < .5*fnorm*fnorm + lambda*vi->alpha*initslope) { /* sufficient reduction */ 2015 if (vi->lsmonitor) { 2016 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: Quadratically determined step, lambda=%18.16e\n",lambda);CHKERRQ(ierr); 2017 } 2018 goto theend1; 2019 } 2020 2021 /* Fit points with cubic */ 2022 count = 1; 2023 while (PETSC_TRUE) { 2024 if (lambda <= minlambda) { 2025 if (vi->lsmonitor) { 2026 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: unable to find good step length! After %D tries \n",count);CHKERRQ(ierr); 2027 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); 2028 } 2029 *flag = PETSC_FALSE; 2030 break; 2031 } 2032 t1 = .5*((*gnorm)*(*gnorm) - fnorm*fnorm) - lambda*initslope; 2033 t2 = .5*(gnormprev*gnormprev - fnorm*fnorm) - lambdaprev*initslope; 2034 a = (t1/(lambda*lambda) - t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev); 2035 b = (-lambdaprev*t1/(lambda*lambda) + lambda*t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev); 2036 d = b*b - 3*a*initslope; 2037 if (d < 0.0) d = 0.0; 2038 if (a == 0.0) { 2039 lambdatemp = -initslope/(2.0*b); 2040 } else { 2041 lambdatemp = (-b + sqrt(d))/(3.0*a); 2042 } 2043 lambdaprev = lambda; 2044 gnormprev = *gnorm; 2045 if (lambdatemp > .5*lambda) lambdatemp = .5*lambda; 2046 if (lambdatemp <= .1*lambda) lambda = .1*lambda; 2047 else lambda = lambdatemp; 2048 ierr = VecWAXPY(w,-lambda,y,x);CHKERRQ(ierr); 2049 ierr = SNESVIProjectOntoBounds(snes,w);CHKERRQ(ierr); 2050 if (snes->nfuncs >= snes->max_funcs) { 2051 ierr = PetscInfo1(snes,"Exceeded maximum function evaluations, while looking for good step length! %D \n",count);CHKERRQ(ierr); 2052 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); 2053 *flag = PETSC_FALSE; 2054 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 2055 break; 2056 } 2057 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 2058 if (snes->ops->solve != SNESSolveVI_SS) { 2059 ierr = SNESVIComputeInactiveSetFnorm(snes,g,w,gnorm);CHKERRQ(ierr); 2060 } else { 2061 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); 2062 } 2063 if (snes->domainerror) { 2064 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 2065 PetscFunctionReturn(0); 2066 } 2067 if (PetscIsInfOrNanReal(*gnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 2068 if (.5*(*gnorm)*(*gnorm) < .5*fnorm*fnorm + lambda*vi->alpha*initslope) { /* is reduction enough? */ 2069 if (vi->lsmonitor) { 2070 ierr = PetscPrintf(comm," Line search: Cubically determined step, current gnorm %G lambda=%18.16e\n",*gnorm,lambda);CHKERRQ(ierr); 2071 } 2072 break; 2073 } else { 2074 if (vi->lsmonitor) { 2075 ierr = PetscPrintf(comm," Line search: Cubic step no good, shrinking lambda, current gnorm %G lambda=%18.16e\n",*gnorm,lambda);CHKERRQ(ierr); 2076 } 2077 } 2078 count++; 2079 } 2080 theend1: 2081 /* Optional user-defined check for line search step validity */ 2082 if (vi->postcheckstep && *flag) { 2083 ierr = (*vi->postcheckstep)(snes,x,y,w,vi->postcheck,&changed_y,&changed_w);CHKERRQ(ierr); 2084 if (changed_y) { 2085 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); 2086 ierr = SNESVIProjectOntoBounds(snes,w);CHKERRQ(ierr); 2087 } 2088 if (changed_y || changed_w) { /* recompute the function if the step has changed */ 2089 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 2090 if (snes->ops->solve != SNESSolveVI_SS) { 2091 ierr = SNESVIComputeInactiveSetFnorm(snes,g,w,gnorm);CHKERRQ(ierr); 2092 } else { 2093 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); 2094 } 2095 if (snes->domainerror) { 2096 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 2097 PetscFunctionReturn(0); 2098 } 2099 if (PetscIsInfOrNanReal(*gnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 2100 ierr = VecNormBegin(y,NORM_2,ynorm);CHKERRQ(ierr); 2101 ierr = VecNormEnd(y,NORM_2,ynorm);CHKERRQ(ierr); 2102 } 2103 } 2104 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 2105 PetscFunctionReturn(0); 2106 } 2107 2108 #undef __FUNCT__ 2109 #define __FUNCT__ "SNESLineSearchQuadratic_VI" 2110 /* 2111 This routine does a quadratic line search while keeping the iterates within the variable bounds 2112 */ 2113 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) 2114 { 2115 /* 2116 Note that for line search purposes we work with with the related 2117 minimization problem: 2118 min z(x): R^n -> R, 2119 where z(x) = .5 * fnorm*fnorm,and fnorm = || f ||_2. 2120 */ 2121 PetscReal initslope,minlambda,lambda,lambdatemp,rellength; 2122 #if defined(PETSC_USE_COMPLEX) 2123 PetscScalar cinitslope; 2124 #endif 2125 PetscErrorCode ierr; 2126 PetscInt count; 2127 SNES_VI *vi = (SNES_VI*)snes->data; 2128 PetscBool changed_w = PETSC_FALSE,changed_y = PETSC_FALSE; 2129 2130 PetscFunctionBegin; 2131 ierr = PetscLogEventBegin(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 2132 *flag = PETSC_TRUE; 2133 2134 ierr = VecNorm(y,NORM_2,ynorm);CHKERRQ(ierr); 2135 if (*ynorm == 0.0) { 2136 if (vi->lsmonitor) { 2137 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor,"Line search: Direction and size is 0\n");CHKERRQ(ierr); 2138 } 2139 *gnorm = fnorm; 2140 ierr = VecCopy(x,w);CHKERRQ(ierr); 2141 ierr = VecCopy(f,g);CHKERRQ(ierr); 2142 *flag = PETSC_FALSE; 2143 goto theend2; 2144 } 2145 if (*ynorm > vi->maxstep) { /* Step too big, so scale back */ 2146 ierr = VecScale(y,vi->maxstep/(*ynorm));CHKERRQ(ierr); 2147 *ynorm = vi->maxstep; 2148 } 2149 ierr = VecMaxPointwiseDivide(y,x,&rellength);CHKERRQ(ierr); 2150 minlambda = vi->minlambda/rellength; 2151 ierr = MatMult(snes->jacobian,y,w);CHKERRQ(ierr); 2152 #if defined(PETSC_USE_COMPLEX) 2153 ierr = VecDot(f,w,&cinitslope);CHKERRQ(ierr); 2154 initslope = PetscRealPart(cinitslope); 2155 #else 2156 ierr = VecDot(f,w,&initslope);CHKERRQ(ierr); 2157 #endif 2158 if (initslope > 0.0) initslope = -initslope; 2159 if (initslope == 0.0) initslope = -1.0; 2160 ierr = PetscInfo1(snes,"Initslope %G \n",initslope);CHKERRQ(ierr); 2161 2162 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); 2163 ierr = SNESVIProjectOntoBounds(snes,w);CHKERRQ(ierr); 2164 if (snes->nfuncs >= snes->max_funcs) { 2165 ierr = PetscInfo(snes,"Exceeded maximum function evaluations, while checking full step length!\n");CHKERRQ(ierr); 2166 *flag = PETSC_FALSE; 2167 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 2168 goto theend2; 2169 } 2170 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 2171 if (snes->ops->solve != SNESSolveVI_SS) { 2172 ierr = SNESVIComputeInactiveSetFnorm(snes,g,w,gnorm);CHKERRQ(ierr); 2173 } else { 2174 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); 2175 } 2176 if (snes->domainerror) { 2177 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 2178 PetscFunctionReturn(0); 2179 } 2180 if (PetscIsInfOrNanReal(*gnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 2181 if (.5*(*gnorm)*(*gnorm) <= .5*fnorm*fnorm + vi->alpha*initslope) { /* Sufficient reduction */ 2182 if (vi->lsmonitor) { 2183 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: Using full step: fnorm %G gnorm %G\n",fnorm,*gnorm);CHKERRQ(ierr); 2184 } 2185 goto theend2; 2186 } 2187 2188 /* Fit points with quadratic */ 2189 lambda = 1.0; 2190 count = 1; 2191 while (PETSC_TRUE) { 2192 if (lambda <= minlambda) { /* bad luck; use full step */ 2193 if (vi->lsmonitor) { 2194 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor,"Line search: Unable to find good step length! %D \n",count);CHKERRQ(ierr); 2195 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); 2196 } 2197 ierr = VecCopy(x,w);CHKERRQ(ierr); 2198 *flag = PETSC_FALSE; 2199 break; 2200 } 2201 lambdatemp = -initslope/((*gnorm)*(*gnorm) - fnorm*fnorm - 2.0*initslope); 2202 if (lambdatemp > .5*lambda) lambdatemp = .5*lambda; 2203 if (lambdatemp <= .1*lambda) lambda = .1*lambda; 2204 else lambda = lambdatemp; 2205 2206 ierr = VecWAXPY(w,-lambda,y,x);CHKERRQ(ierr); 2207 ierr = SNESVIProjectOntoBounds(snes,w);CHKERRQ(ierr); 2208 if (snes->nfuncs >= snes->max_funcs) { 2209 ierr = PetscInfo1(snes,"Exceeded maximum function evaluations, while looking for good step length! %D \n",count);CHKERRQ(ierr); 2210 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); 2211 *flag = PETSC_FALSE; 2212 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 2213 break; 2214 } 2215 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 2216 if (snes->domainerror) { 2217 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 2218 PetscFunctionReturn(0); 2219 } 2220 if (snes->ops->solve != SNESSolveVI_SS) { 2221 ierr = SNESVIComputeInactiveSetFnorm(snes,g,w,gnorm);CHKERRQ(ierr); 2222 } else { 2223 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); 2224 } 2225 if (PetscIsInfOrNanReal(*gnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 2226 if (.5*(*gnorm)*(*gnorm) < .5*fnorm*fnorm + lambda*vi->alpha*initslope) { /* sufficient reduction */ 2227 if (vi->lsmonitor) { 2228 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line Search: Quadratically determined step, lambda=%G\n",lambda);CHKERRQ(ierr); 2229 } 2230 break; 2231 } 2232 count++; 2233 } 2234 theend2: 2235 /* Optional user-defined check for line search step validity */ 2236 if (vi->postcheckstep) { 2237 ierr = (*vi->postcheckstep)(snes,x,y,w,vi->postcheck,&changed_y,&changed_w);CHKERRQ(ierr); 2238 if (changed_y) { 2239 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); 2240 ierr = SNESVIProjectOntoBounds(snes,w);CHKERRQ(ierr); 2241 } 2242 if (changed_y || changed_w) { /* recompute the function if the step has changed */ 2243 ierr = SNESComputeFunction(snes,w,g); 2244 if (snes->domainerror) { 2245 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 2246 PetscFunctionReturn(0); 2247 } 2248 if (snes->ops->solve != SNESSolveVI_SS) { 2249 ierr = SNESVIComputeInactiveSetFnorm(snes,g,w,gnorm);CHKERRQ(ierr); 2250 } else { 2251 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); 2252 } 2253 2254 ierr = VecNormBegin(y,NORM_2,ynorm);CHKERRQ(ierr); 2255 ierr = VecNormEnd(y,NORM_2,ynorm);CHKERRQ(ierr); 2256 if (PetscIsInfOrNanReal(*gnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 2257 } 2258 } 2259 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 2260 PetscFunctionReturn(0); 2261 } 2262 2263 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*/ 2264 /* -------------------------------------------------------------------------- */ 2265 EXTERN_C_BEGIN 2266 #undef __FUNCT__ 2267 #define __FUNCT__ "SNESLineSearchSet_VI" 2268 PetscErrorCode SNESLineSearchSet_VI(SNES snes,FCN2 func,void *lsctx) 2269 { 2270 PetscFunctionBegin; 2271 ((SNES_VI *)(snes->data))->LineSearch = func; 2272 ((SNES_VI *)(snes->data))->lsP = lsctx; 2273 PetscFunctionReturn(0); 2274 } 2275 EXTERN_C_END 2276 2277 /* -------------------------------------------------------------------------- */ 2278 EXTERN_C_BEGIN 2279 #undef __FUNCT__ 2280 #define __FUNCT__ "SNESLineSearchSetMonitor_VI" 2281 PetscErrorCode SNESLineSearchSetMonitor_VI(SNES snes,PetscBool flg) 2282 { 2283 SNES_VI *vi = (SNES_VI*)snes->data; 2284 PetscErrorCode ierr; 2285 2286 PetscFunctionBegin; 2287 if (flg && !vi->lsmonitor) { 2288 ierr = PetscViewerASCIIMonitorCreate(((PetscObject)snes)->comm,"stdout",((PetscObject)snes)->tablevel,&vi->lsmonitor);CHKERRQ(ierr); 2289 } else if (!flg && vi->lsmonitor) { 2290 ierr = PetscViewerASCIIMonitorDestroy(&vi->lsmonitor);CHKERRQ(ierr); 2291 } 2292 PetscFunctionReturn(0); 2293 } 2294 EXTERN_C_END 2295 2296 /* 2297 SNESView_VI - Prints info from the SNESVI data structure. 2298 2299 Input Parameters: 2300 . SNES - the SNES context 2301 . viewer - visualization context 2302 2303 Application Interface Routine: SNESView() 2304 */ 2305 #undef __FUNCT__ 2306 #define __FUNCT__ "SNESView_VI" 2307 static PetscErrorCode SNESView_VI(SNES snes,PetscViewer viewer) 2308 { 2309 SNES_VI *vi = (SNES_VI *)snes->data; 2310 const char *cstr,*tstr; 2311 PetscErrorCode ierr; 2312 PetscBool iascii; 2313 2314 PetscFunctionBegin; 2315 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 2316 if (iascii) { 2317 if (vi->LineSearch == SNESLineSearchNo_VI) cstr = "SNESLineSearchNo"; 2318 else if (vi->LineSearch == SNESLineSearchQuadratic_VI) cstr = "SNESLineSearchQuadratic"; 2319 else if (vi->LineSearch == SNESLineSearchCubic_VI) cstr = "SNESLineSearchCubic"; 2320 else cstr = "unknown"; 2321 if (snes->ops->solve == SNESSolveVI_SS) tstr = "Semismooth"; 2322 else if (snes->ops->solve == SNESSolveVI_RS) tstr = "Reduced Space"; 2323 else if (snes->ops->solve == SNESSolveVI_RSAUG) tstr = "Reduced space with augmented variables"; 2324 else tstr = "unknown"; 2325 ierr = PetscViewerASCIIPrintf(viewer," VI algorithm: %s\n",tstr);CHKERRQ(ierr); 2326 ierr = PetscViewerASCIIPrintf(viewer," line search variant: %s\n",cstr);CHKERRQ(ierr); 2327 ierr = PetscViewerASCIIPrintf(viewer," alpha=%G, maxstep=%G, minlambda=%G\n",vi->alpha,vi->maxstep,vi->minlambda);CHKERRQ(ierr); 2328 } else { 2329 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type %s not supported for SNES EQ VI",((PetscObject)viewer)->type_name); 2330 } 2331 PetscFunctionReturn(0); 2332 } 2333 2334 #undef __FUNCT__ 2335 #define __FUNCT__ "SNESVISetVariableBounds" 2336 /*@ 2337 SNESVISetVariableBounds - Sets the lower and upper bounds for the solution vector. xl <= x <= xu. 2338 2339 Input Parameters: 2340 . snes - the SNES context. 2341 . xl - lower bound. 2342 . xu - upper bound. 2343 2344 Notes: 2345 If this routine is not called then the lower and upper bounds are set to 2346 SNES_VI_INF and SNES_VI_NINF respectively during SNESSetUp(). 2347 2348 @*/ 2349 PetscErrorCode SNESVISetVariableBounds(SNES snes, Vec xl, Vec xu) 2350 { 2351 SNES_VI *vi; 2352 PetscErrorCode ierr; 2353 2354 PetscFunctionBegin; 2355 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2356 PetscValidHeaderSpecific(xl,VEC_CLASSID,2); 2357 PetscValidHeaderSpecific(xu,VEC_CLASSID,3); 2358 if (!snes->vec_func) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call SNESSetFunction() first"); 2359 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); 2360 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); 2361 ierr = SNESSetType(snes,SNESVI);CHKERRQ(ierr); 2362 vi = (SNES_VI*)snes->data; 2363 ierr = PetscObjectReference((PetscObject)xl);CHKERRQ(ierr); 2364 ierr = PetscObjectReference((PetscObject)xu);CHKERRQ(ierr); 2365 ierr = VecDestroy(&vi->xl);CHKERRQ(ierr); 2366 ierr = VecDestroy(&vi->xu);CHKERRQ(ierr); 2367 vi->xl = xl; 2368 vi->xu = xu; 2369 PetscFunctionReturn(0); 2370 } 2371 2372 /* -------------------------------------------------------------------------- */ 2373 /* 2374 SNESSetFromOptions_VI - Sets various parameters for the SNESVI method. 2375 2376 Input Parameter: 2377 . snes - the SNES context 2378 2379 Application Interface Routine: SNESSetFromOptions() 2380 */ 2381 #undef __FUNCT__ 2382 #define __FUNCT__ "SNESSetFromOptions_VI" 2383 static PetscErrorCode SNESSetFromOptions_VI(SNES snes) 2384 { 2385 SNES_VI *vi = (SNES_VI *)snes->data; 2386 const char *lses[] = {"basic","basicnonorms","quadratic","cubic"}; 2387 const char *vies[] = {"ss","rs","rsaug"}; 2388 PetscErrorCode ierr; 2389 PetscInt indx; 2390 PetscBool flg,set,flg2; 2391 2392 PetscFunctionBegin; 2393 ierr = PetscOptionsHead("SNES semismooth method options");CHKERRQ(ierr); 2394 ierr = PetscOptionsBool("-snes_vi_monitor","Monitor all non-active variables","None",PETSC_FALSE,&flg,0);CHKERRQ(ierr); 2395 if (flg) { 2396 ierr = SNESMonitorSet(snes,SNESMonitorVI,0,0);CHKERRQ(ierr); 2397 } 2398 ierr = PetscOptionsReal("-snes_ls_alpha","Function norm must decrease by","None",vi->alpha,&vi->alpha,0);CHKERRQ(ierr); 2399 ierr = PetscOptionsReal("-snes_ls_maxstep","Step must be less than","None",vi->maxstep,&vi->maxstep,0);CHKERRQ(ierr); 2400 ierr = PetscOptionsReal("-snes_ls_minlambda","Minimum lambda allowed","None",vi->minlambda,&vi->minlambda,0);CHKERRQ(ierr); 2401 ierr = PetscOptionsReal("-snes_vi_const_tol","constraint tolerance","None",vi->const_tol,&vi->const_tol,0);CHKERRQ(ierr); 2402 ierr = PetscOptionsBool("-snes_ls_monitor","Print progress of line searches","SNESLineSearchSetMonitor",vi->lsmonitor ? PETSC_TRUE : PETSC_FALSE,&flg,&set);CHKERRQ(ierr); 2403 if (set) {ierr = SNESLineSearchSetMonitor(snes,flg);CHKERRQ(ierr);} 2404 ierr = PetscOptionsEList("-snes_vi_type","Semismooth algorithm used","",vies,3,"ss",&indx,&flg2);CHKERRQ(ierr); 2405 if (flg2) { 2406 switch (indx) { 2407 case 0: 2408 snes->ops->solve = SNESSolveVI_SS; 2409 break; 2410 case 1: 2411 snes->ops->solve = SNESSolveVI_RS; 2412 break; 2413 case 2: 2414 snes->ops->solve = SNESSolveVI_RSAUG; 2415 } 2416 } 2417 ierr = PetscOptionsEList("-snes_ls","Line search used","SNESLineSearchSet",lses,4,"basic",&indx,&flg);CHKERRQ(ierr); 2418 if (flg) { 2419 switch (indx) { 2420 case 0: 2421 ierr = SNESLineSearchSet(snes,SNESLineSearchNo_VI,PETSC_NULL);CHKERRQ(ierr); 2422 break; 2423 case 1: 2424 ierr = SNESLineSearchSet(snes,SNESLineSearchNoNorms_VI,PETSC_NULL);CHKERRQ(ierr); 2425 break; 2426 case 2: 2427 ierr = SNESLineSearchSet(snes,SNESLineSearchQuadratic_VI,PETSC_NULL);CHKERRQ(ierr); 2428 break; 2429 case 3: 2430 ierr = SNESLineSearchSet(snes,SNESLineSearchCubic_VI,PETSC_NULL);CHKERRQ(ierr); 2431 break; 2432 } 2433 } 2434 ierr = PetscOptionsTail();CHKERRQ(ierr); 2435 PetscFunctionReturn(0); 2436 } 2437 /* -------------------------------------------------------------------------- */ 2438 /*MC 2439 SNESVI - Various solvers for variational inequalities based on Newton's method 2440 2441 Options Database: 2442 + -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 2443 additional variables that enforce the constraints 2444 - -snes_vi_monitor - prints the number of active constraints at each iteration. 2445 2446 2447 Level: beginner 2448 2449 .seealso: SNESCreate(), SNES, SNESSetType(), SNESTR, SNESLineSearchSet(), 2450 SNESLineSearchSetPostCheck(), SNESLineSearchNo(), SNESLineSearchCubic(), SNESLineSearchQuadratic(), 2451 SNESLineSearchSet(), SNESLineSearchNoNorms(), SNESLineSearchSetPreCheck(), SNESLineSearchSetParams(), SNESLineSearchGetParams() 2452 2453 M*/ 2454 EXTERN_C_BEGIN 2455 #undef __FUNCT__ 2456 #define __FUNCT__ "SNESCreate_VI" 2457 PetscErrorCode SNESCreate_VI(SNES snes) 2458 { 2459 PetscErrorCode ierr; 2460 SNES_VI *vi; 2461 2462 PetscFunctionBegin; 2463 snes->ops->reset = SNESReset_VI; 2464 snes->ops->setup = SNESSetUp_VI; 2465 snes->ops->solve = SNESSolveVI_RS; 2466 snes->ops->destroy = SNESDestroy_VI; 2467 snes->ops->setfromoptions = SNESSetFromOptions_VI; 2468 snes->ops->view = SNESView_VI; 2469 snes->ops->converged = SNESDefaultConverged_VI; 2470 2471 ierr = PetscNewLog(snes,SNES_VI,&vi);CHKERRQ(ierr); 2472 snes->data = (void*)vi; 2473 vi->alpha = 1.e-4; 2474 vi->maxstep = 1.e8; 2475 vi->minlambda = 1.e-12; 2476 vi->LineSearch = SNESLineSearchNo_VI; 2477 vi->lsP = PETSC_NULL; 2478 vi->postcheckstep = PETSC_NULL; 2479 vi->postcheck = PETSC_NULL; 2480 vi->precheckstep = PETSC_NULL; 2481 vi->precheck = PETSC_NULL; 2482 vi->const_tol = 2.2204460492503131e-16; 2483 vi->checkredundancy = PETSC_NULL; 2484 2485 ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESLineSearchSetMonitor_C","SNESLineSearchSetMonitor_VI",SNESLineSearchSetMonitor_VI);CHKERRQ(ierr); 2486 ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESLineSearchSet_C","SNESLineSearchSet_VI",SNESLineSearchSet_VI);CHKERRQ(ierr); 2487 2488 PetscFunctionReturn(0); 2489 } 2490 EXTERN_C_END 2491