1 /* Defines the basic SNES object */ 2 #include <private/snesimpl.h> 3 4 static PetscErrorCode BuildNGmresSoln(PetscScalar*,Vec,SNES,PetscInt,PetscInt); 5 6 /* Private structure for the Anderson mixing method aka nonlinear Krylov */ 7 typedef struct { 8 PetscScalar *hh_origin; 9 Vec *v, *w, *q; 10 PetscReal *f2; /* 2-norms of function (residual) at each stage */ 11 PetscInt msize; /* maximum size of space */ 12 PetscInt csize; /* current size of space */ 13 PetscScalar beta; /* relaxation parameter */ 14 PetscScalar *nrs; /* temp that holds the coefficients of the Krylov vectors that form the minimum residual solution */ 15 16 } SNES_NGMRES; 17 18 #define HH(a,b) (ngmres->hh_origin + (a)*(ngmres->msize)+(b)) 19 20 #undef __FUNCT__ 21 #define __FUNCT__ "SNESReset_NGMRES" 22 PetscErrorCode SNESReset_NGMRES(SNES snes) 23 { 24 SNES_NGMRES *ngmres = (SNES_NGMRES*) snes->data; 25 PetscErrorCode ierr; 26 27 PetscFunctionBegin; 28 ierr = VecDestroyVecs(ngmres->msize, &ngmres->v);CHKERRQ(ierr); 29 ierr = VecDestroyVecs(ngmres->msize, &ngmres->w);CHKERRQ(ierr); 30 ierr = VecDestroyVecs(ngmres->msize, &ngmres->q);CHKERRQ(ierr); 31 PetscFunctionReturn(0); 32 } 33 34 #undef __FUNCT__ 35 #define __FUNCT__ "SNESDestroy_NGMRES" 36 PetscErrorCode SNESDestroy_NGMRES(SNES snes) 37 { 38 PetscErrorCode ierr; 39 40 PetscFunctionBegin; 41 ierr = SNESReset_NGMRES(snes);CHKERRQ(ierr); 42 if (snes->work) {ierr = VecDestroyVecs(snes->nwork, &snes->work);CHKERRQ(ierr);} 43 PetscFunctionReturn(0); 44 } 45 46 #undef __FUNCT__ 47 #define __FUNCT__ "SNESSetUp_NGMRES" 48 PetscErrorCode SNESSetUp_NGMRES(SNES snes) 49 { 50 SNES_NGMRES *ngmres = (SNES_NGMRES *) snes->data; 51 PetscInt msize,hh; 52 PetscErrorCode ierr; 53 54 PetscFunctionBegin; 55 #if 0 56 if (snes->pc_side != PC_LEFT) {SETERRQ(((PetscObject) snes)->comm, PETSC_ERR_SUP, "Only left preconditioning allowed for SNESNGMRES");} 57 #endif 58 ngmres->beta = 1.0; 59 msize = ngmres->msize; /* restart size */ 60 hh = msize * msize; 61 ierr = PetscMalloc2(hh,PetscScalar,&ngmres->hh_origin,msize,PetscScalar,&ngmres->nrs);CHKERRQ(ierr); 62 63 ierr = PetscMemzero(ngmres->hh_origin,hh*sizeof(PetscScalar));CHKERRQ(ierr); 64 ierr = PetscMemzero(ngmres->nrs,msize*sizeof(PetscScalar));CHKERRQ(ierr); 65 66 // ierr = PetscLogObjectMemory(ksp,(hh+msize)*sizeof(PetscScalar));CHKERRQ(ierr); 67 68 ierr = VecDuplicateVecs(snes->vec_sol, ngmres->msize, &ngmres->v);CHKERRQ(ierr); 69 ierr = VecDuplicateVecs(snes->vec_sol, ngmres->msize, &ngmres->w);CHKERRQ(ierr); 70 ierr = VecDuplicateVecs(snes->vec_sol, ngmres->msize, &ngmres->q);CHKERRQ(ierr); 71 ierr = SNESDefaultGetWork(snes, 2);CHKERRQ(ierr); 72 PetscFunctionReturn(0); 73 } 74 75 #undef __FUNCT__ 76 #define __FUNCT__ "SNESSetFromOptions_NGMRES" 77 PetscErrorCode SNESSetFromOptions_NGMRES(SNES snes) 78 { 79 SNES_NGMRES *ngmres = (SNES_NGMRES *) snes->data; 80 PetscErrorCode ierr; 81 82 PetscFunctionBegin; 83 ierr = PetscOptionsHead("SNES NGMRES options");CHKERRQ(ierr); 84 ierr = PetscOptionsInt("-snes_ngmres_restart", "Number of directions", "SNES", ngmres->msize, &ngmres->msize, PETSC_NULL);CHKERRQ(ierr); 85 ierr = PetscOptionsTail();CHKERRQ(ierr); 86 PetscFunctionReturn(0); 87 } 88 89 #undef __FUNCT__ 90 #define __FUNCT__ "SNESView_NGMRES" 91 PetscErrorCode SNESView_NGMRES(SNES snes, PetscViewer viewer) 92 { 93 SNES_NGMRES *ngmres = (SNES_NGMRES *) snes->data; 94 PetscBool iascii; 95 PetscErrorCode ierr; 96 97 PetscFunctionBegin; 98 ierr = PetscTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);CHKERRQ(ierr); 99 if (iascii) { 100 ierr = PetscViewerASCIIPrintf(viewer, " Size of space %d\n", ngmres->msize);CHKERRQ(ierr); 101 } 102 PetscFunctionReturn(0); 103 } 104 105 #undef __FUNCT__ 106 #define __FUNCT__ "SNESSolve_NGMRES" 107 PetscErrorCode SNESSolve_NGMRES(SNES snes) 108 { 109 SNES pc; 110 SNES_NGMRES *ngmres = (SNES_NGMRES *) snes->data; 111 Vec X,F, Fold, Xold,temp,*dX = ngmres->v,*dF = ngmres->w; 112 PetscScalar *nrs=ngmres->nrs; 113 PetscReal gnorm; 114 PetscInt i, j, k,ivec, l, flag, it; 115 PetscErrorCode ierr; 116 117 PetscFunctionBegin; 118 snes->reason = SNES_CONVERGED_ITERATING; 119 X = snes->vec_sol; 120 F = snes->vec_func; 121 Fold = snes->work[0]; 122 Xold = snes->work[1]; 123 ierr = VecDuplicate(X,&Xold);CHKERRQ(ierr); 124 ierr = VecDuplicate(X,&temp);CHKERRQ(ierr); 125 126 ierr = SNESGetPC(snes, &pc);CHKERRQ(ierr); 127 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 128 snes->iter = 0; 129 snes->norm = 0.; 130 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 131 ierr = SNESComputeFunction(snes, X, temp);CHKERRQ(ierr); /* r = F(x) */ 132 #if 0 133 ierr = SNESSolve(snes->pc, temp, F);CHKERRQ(ierr); /* p = P(r) */ 134 #else 135 ierr = VecCopy(temp, F);CHKERRQ(ierr); /* p = r */ 136 #endif 137 138 if (snes->domainerror) { 139 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 140 PetscFunctionReturn(0); 141 } 142 ierr = VecNorm(F, NORM_2, &gnorm);CHKERRQ(ierr); /* fnorm = ||r|| */ 143 if (PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FP, "Infinite or not-a-number generated in norm"); 144 145 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 146 snes->norm = gnorm; 147 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 148 SNESLogConvHistory(snes, gnorm, 0); 149 ierr = SNESMonitor(snes, 0, gnorm);CHKERRQ(ierr); 150 151 /* set parameter for default relative tolerance convergence test */ 152 snes->ttol = gnorm*snes->rtol; 153 /* test convergence */ 154 ierr = (*snes->ops->converged)(snes,0,0.0,0.0,gnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 155 if (snes->reason) PetscFunctionReturn(0); 156 157 /* for k=0 */ 158 159 k=0; 160 ierr= VecCopy(X,Xold); CHKERRQ(ierr); /* Xbar_0= X_0 */ 161 ierr= VecCopy(F,Fold); CHKERRQ(ierr); /* Fbar_0 = f_0= B(b-Ax_0) */ 162 ierr= VecWAXPY(X, ngmres->beta, Fold, Xold); CHKERRQ(ierr); /*X_1 = X_bar + beta * Fbar */ 163 164 /* to calculate f_1 */ 165 ierr = SNESComputeFunction(snes, X, temp);CHKERRQ(ierr); /* r = F(x) */ 166 #if 0 167 ierr = SNESSolve(snes->pc, temp, F);CHKERRQ(ierr); /* p = P(r) */ 168 #else 169 ierr = VecCopy(temp, F);CHKERRQ(ierr); /* p = r */ 170 #endif 171 172 /* calculate dX and dF for k=0 */ 173 ierr= VecWAXPY(dX[k],-1.0, Xold, X); CHKERRQ(ierr); /* dX= X_1 - X_0 */ 174 ierr= VecWAXPY(dF[k],-1.0, Fold, F); CHKERRQ(ierr); /* dF= f_1 - f_0 */ 175 176 177 ierr= VecCopy(X,Xold); CHKERRQ(ierr); /* Xbar_0= X_0 */ 178 ierr= VecCopy(F,Fold); CHKERRQ(ierr); /* Fbar_0 = f_0= B(b-Ax_0) */ 179 180 flag=0; 181 182 183 184 for (k=1; k<snes->max_its; k += 1) { /* begin the iteration */ 185 186 printf("iteration %d\n",k); 187 l=ngmres->msize; 188 189 if(k<l) { 190 l=k; 191 ivec=k; 192 } 193 else{ 194 ivec=l-1; 195 } 196 197 198 it=l-1; 199 ierr = BuildNGmresSoln(nrs,Fold,snes,it,flag);CHKERRQ(ierr); 200 201 202 /* to obtain the solution at k+1 step */ 203 ierr= VecCopy(Xold,X); CHKERRQ(ierr); /* X=Xold+Fold-(dX + dF) *nrd */ 204 ierr= VecAXPY(X,1.0,Fold); CHKERRQ(ierr); /* X= Xold+Fold */ 205 for(i=0;i<l;i++){ /*X= Xold+Fold- (dX+dF*beta) *innerb */ 206 ierr= VecAXPY(X,-nrs[i], dX[i]);CHKERRQ(ierr); 207 ierr= VecAXPY(X,-nrs[i]*ngmres->beta, dF[i]);CHKERRQ(ierr); 208 } 209 210 211 212 213 /* to calculate f_k+1 */ 214 ierr = SNESComputeFunction(snes, X, temp);CHKERRQ(ierr); /* r = F(x) */ 215 #if 0 216 ierr = SNESSolve(snes->pc, temp, F);CHKERRQ(ierr); /* p = P(r) */ 217 #else 218 ierr = VecCopy(temp, F);CHKERRQ(ierr); /* p = r */ 219 #endif 220 221 222 223 /* check the convegence */ 224 ierr = VecNorm(F, NORM_2, &gnorm);CHKERRQ(ierr); /* fnorm = ||r|| */ 225 SNESLogConvHistory(snes, gnorm, k); 226 ierr = SNESMonitor(snes, k, gnorm);CHKERRQ(ierr); 227 snes->iter =k; 228 ierr = (*snes->ops->converged)(snes,0,0.0,0.0,gnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 229 230 if (snes->reason) PetscFunctionReturn(0); 231 232 233 234 /* calculate dX and dF for k=0 */ 235 if( k>ivec) {/* we need to replace the old vectors */ 236 flag=1; 237 for(i=0;i<it;i++){ 238 ierr= VecCopy(dX[i+1],dX[i]); CHKERRQ(ierr); /* X=Xold+Fold-(dX + dF) *nrd */ 239 ierr= VecCopy(dF[i+1],dF[i]); CHKERRQ(ierr); /* X=Xold+Fold-(dX + dF) *nrd */ 240 for(j=0;j<l;j++) 241 *HH(j,i)=*HH(j,i+1); 242 } 243 } 244 245 ierr= VecWAXPY(dX[ivec],-1.0, Xold, X); CHKERRQ(ierr); /* dX= X_1 - X_0 */ 246 ierr= VecWAXPY(dF[ivec],-1.0, Fold, F); CHKERRQ(ierr); /* dF= f_1 - f_0 */ 247 ierr= VecCopy(X,Xold); CHKERRQ(ierr); 248 ierr= VecCopy(F,Fold); CHKERRQ(ierr); 249 250 } 251 snes->reason = SNES_DIVERGED_MAX_IT; 252 PetscFunctionReturn(0); 253 } 254 255 /*MC 256 SNESNGMRES - The Nonlinear Generalized Minimum Residual (NGMRES) method of Oosterlee and Washio. 257 258 Level: beginner 259 260 Notes: Supports only left preconditioning 261 262 "Krylov Subspace Acceleration of Nonlinear Multigrid with Application to Recirculating Flows", C. W. Oosterlee and T. Washio, 263 SIAM Journal on Scientific Computing, 21(5), 2000. 264 265 .seealso: SNESCreate(), SNES, SNESSetType(), SNESType (for list of available types) 266 M*/ 267 EXTERN_C_BEGIN 268 #undef __FUNCT__ 269 #define __FUNCT__ "SNESCreate_NGMRES" 270 PetscErrorCode SNESCreate_NGMRES(SNES snes) 271 { 272 SNES_NGMRES *ngmres; 273 PetscErrorCode ierr; 274 275 PetscFunctionBegin; 276 snes->ops->destroy = SNESDestroy_NGMRES; 277 snes->ops->setup = SNESSetUp_NGMRES; 278 snes->ops->setfromoptions = SNESSetFromOptions_NGMRES; 279 snes->ops->view = SNESView_NGMRES; 280 snes->ops->solve = SNESSolve_NGMRES; 281 snes->ops->reset = SNESReset_NGMRES; 282 283 ierr = PetscNewLog(snes, SNES_NGMRES, &ngmres);CHKERRQ(ierr); 284 snes->data = (void*) ngmres; 285 ngmres->msize = 30; 286 ngmres->csize = 0; 287 288 ierr = SNESGetPC(snes, &snes->pc);CHKERRQ(ierr); 289 #if 0 290 if (ksp->pc_side != PC_LEFT) {ierr = PetscInfo(ksp,"WARNING! Setting PC_SIDE for NGMRES to left!\n");CHKERRQ(ierr);} 291 snes->pc_side = PC_LEFT; 292 #endif 293 PetscFunctionReturn(0); 294 } 295 EXTERN_C_END 296 297 298 299 /* 300 BuildNGmresSoln - create the solution of the least square problem to determine the coefficients 301 302 Input parameters: 303 nrs - the solution 304 Fold - the RHS vector 305 flag - =1, we need to replace the old vector by new one, =0, we still add new vector 306 This is an internal routine that knows about the NGMRES internals. 307 */ 308 #undef __FUNCT__ 309 #define __FUNCT__ "BuildNGmresSoln" 310 static PetscErrorCode BuildNGmresSoln(PetscScalar* nrs, Vec Fold, SNES snes,PetscInt it, PetscInt flag) 311 { 312 PetscScalar tt,temps; 313 PetscErrorCode ierr; 314 PetscInt i,ii,j,l; 315 SNES_NGMRES *ngmres = (SNES_NGMRES *)(snes->data); 316 Vec *dF=ngmres->w, *Q=ngmres->q,temp; 317 PetscReal gam,areal; 318 PetscScalar a,b,c,s; 319 320 PetscFunctionBegin; 321 ierr = VecDuplicate(Fold,&temp);CHKERRQ(ierr); 322 l=it+1; 323 324 /* Solve for solution vector that minimizes the residual */ 325 326 if(flag==1) { // we need to replace the old vector and need to modify the QR factors, use Givens rotation 327 for(i=0;i<it;i++){ 328 /* calculate the Givens rotation */ 329 a=*HH(i,i); 330 b=*HH(i+1,i); 331 gam=PetscRealPart(PetscSqrtScalar(PetscConj(a)*a+PetscConj(b)*b)); 332 if (gam!= 0.0) { 333 gam=1.0/gam; 334 } else { 335 snes->reason = SNES_DIVERGED_MAX_IT; 336 ierr = PetscInfo2(snes,"Likely your matrix or preconditioner is singular. HH(it,it) is identically zero; it = %D nrs(it) = %G",it,10); 337 PetscFunctionReturn(0); 338 } 339 c= a*gam; 340 s= b*gam; 341 342 #if defined(PETSC_USE_COMPLEX) 343 /* update the Q factor */ 344 ierr= VecCopy(Q[i],temp); CHKERRQ(ierr); 345 ierr = VecAXPBY(temp,s,PetscConj(c),Q[i+1]);CHKERRQ(ierr); /*temp= c*Q[i]+s*Q[i+1] */ 346 ierr = VecAXPBY(Q[i+1],-s,c,Q[i]);CHKERRQ(ierr); /* Q[i+1]= -s*Q[i] + c*Q[i+1] */ 347 ierr= VecCopy(temp,Q[i]); CHKERRQ(ierr); /* Q[i]= c*Q[i] + s*Q[i+1] */ 348 /* update the R factor */ 349 for(j=0;j<l;j++){ 350 a= *HH(i,j); 351 b=*HH(i+1,j); 352 temps=PetscConj(c)* a+s* b; 353 *HH(i+1,j)=-s*a+c*b; 354 *HH(i,j)=temps; 355 } 356 #else 357 /* update the Q factor */ 358 ierr= VecCopy(Q[i],temp); CHKERRQ(ierr); 359 ierr = VecAXPBY(temp,s,c,Q[i+1]);CHKERRQ(ierr); /*temp= c*Q[i]+s*Q[i+1] */ 360 ierr = VecAXPBY(Q[i+1],-s,c,Q[i]);CHKERRQ(ierr); /* Q[i+1]= -s*Q[i] + c*Q[i+1] */ 361 ierr= VecCopy(temp,Q[i]); CHKERRQ(ierr); /* Q[i]= c*Q[i] + s*Q[i+1] */ 362 /* update the R factor */ 363 for(j=0;j<l;j++){ 364 a= *HH(i,j); 365 b=*HH(i+1,j); 366 temps=c* a+s* b; 367 *HH(i+1,j)=-s*a+c*b; 368 *HH(i,j)=temps; 369 } 370 #endif 371 372 } 373 } 374 375 // add a new vector, use modified Gram-Schmidt 376 ierr= VecCopy(dF[it],temp); CHKERRQ(ierr); 377 for(i=0;i<it;i++){ 378 ierr=VecDot(temp,Q[i],HH(i,it));CHKERRQ(ierr); /* h(i,l-1)= dF[l-1]'*Q[i] */ 379 ierr = VecAXPBY(temp,-*HH(i,it),1.0,Q[i]);CHKERRQ(ierr); /* temp= temp- h(i,l-1)*Q[i] */ 380 } 381 382 ierr=VecCopy(temp,Q[it]);CHKERRQ(ierr); 383 384 ierr=VecNormalize(Q[it],&areal);CHKERRQ(ierr); 385 386 *HH(it,it) = areal; 387 388 389 /* modify the RHS with Q'*Fold*/ 390 391 for(i=0;i<l;i++) 392 ierr=VecDot(Fold,Q[i],ngmres->nrs+i);CHKERRQ(ierr); /* nrs= Fold'*Q[i] */ 393 394 /* start backsubstitution to solve the least square problem */ 395 396 if (*HH(it,it) != 0.0) { 397 398 nrs[it] = nrs[it]/ *HH(it,it); 399 } else { 400 snes->reason = SNES_DIVERGED_MAX_IT; 401 ierr = PetscInfo2(snes,"Likely your matrix or preconditioner is singular. HH(it,it) is identically zero; it = %D nrs(it) = %G",it,10); 402 PetscFunctionReturn(0); 403 } 404 for (ii=1; ii<=it; ii++) { 405 i = it - ii; 406 tt = nrs[i]; 407 for (j=i+1; j<=it; j++) tt = tt - *HH(i,j) * nrs[j]; 408 if (*HH(i,i) == 0.0) { 409 snes->reason = SNES_DIVERGED_MAX_IT; 410 ierr = PetscInfo1(snes,"Likely your matrix or preconditioner is singular. HH(k,k) is identically zero; i = %D",i); 411 PetscFunctionReturn(0); 412 } 413 nrs[i] = tt / *HH(i,i); 414 } 415 416 417 PetscFunctionReturn(0); 418 } 419