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