19.6 Multigrid Methods for Boundary Value Problems 879 unsigned int j,jcycle,jj,jpost,jpre,nf,ng=0,ngrid,nn; double **ires [NGMAX+1],**irho [NGMAX+1],**irhs [NGMAX+1],**iu[NGMAX+1]; nn=n; while (nn >>1)ng++; if (n !1+(1L <ng))nrerror("n-1 must be a power of 2 in mglin."); if (ng NGMAX)nrerror("increase NGMAX in mglin.") nn=n/2+1; ngrid=ng-1; irho[ngrid]=dmatrix(1,nn,1,nn); Allocate storage for r.h.s.on grid ng-1, rstrct(irho[ngrid],u,nn); and fill it by restricting from the fine grid. while (nn 3){ Similarly allocate storage and fill r.h.s.on all nn=nn/2+1; coarse grids. irho[--ngrid]=dmatrix(1,nn,1,nn); rstrct(irho[ngrid],irho[ngrid+1],nn); (including this one) granted for 18881992 nn=3: iu[1]=dmatrix(1,nn,1,nn); 1600 irhs[1]=dmatrix(1,nn,1,nn) slvsml(iu[1],irho[1]); Initial solution on coarsest grid. 872 free_dmatrix(irho[1],1,nn,1,nn); ngrid-ng; 7423 for (j=2;j<=ngrid;j++){ Nested iteration loop. nn=2*nn-1; iu[j]=dmatrix(1,nn,1,nn); irhs[i]=dmatrix(1,nn,1,nn); (North America to any server computer, from NUMERICAL RECIPES IN ires[j]=dmatrix(1,nn,1,nn); to make one paper Cambridge University Press. THE interp(iu[j],iu[i-1],nn); 是 ART Interpolate from coarse grid to next finer grid. copy(irhs[j],(j !ngrid irho[j]u),nn); Set up r.h.s. for (jcycle=1;jcycle<=ncycle;jcycle++){ V-cycle loop Programs nf=nn; for (jj=j;jj>=2;jj--){ Downward stoke of the V. for (ipre=1;ipre<=NPRE;ipre++) Pre-smoothing. strictly prohibited relax(iu[jj],irhs[ji],nf); resid(ires[jj],iu[jj],irhs[jj],nf); to dir nf=nf/2+1: rstrct(irhs[jj-1],ires[jj],nf); ectc Restriction of the residual is the next r.h.s 18881920 OF SCIENTIFIC COMPUTING(ISBN f1110(iu[jj-1],nf); Zero for initial guess in next relaxation. slvsml(iu[1],irhs[1]); Bottom of V:solve on coars- nf=3; est grid. 10-621 for (jj=2;jj<=j;jj++) Upward stroke of V. nf=2*nf-1; addint(iu[jj],iu[jj-1],ires[jj],nf); .Further reproduction, Numerical Recipes -43108 Use res for temporary storage inside addint. for (jpost=1;jpost<=NPOST;ipost++) Post-smoothing relax(iu[jj],irhs[jj],nf); (outside Software. copy(u,iu[ngrid],n); Return solution in u. for (nn=n,j=ng;j>=2;j--,nn=nn/2+1){ free_dmatrix(ires[i],1,nn,1,nn); free_dmatrix(irhs[j],1,nn,1,nn); free_dmatrix(iu[j],1,nn,1,nn); if (j !ng)free_dmatrix(irho[j],1,nn,1,nn); 2 free_dmatrix(irhs[1],1,3,1,3); free_dmatrix(iu[1],1,3,1,3);19.6 Multigrid Methods for Boundary Value Problems 879 Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machine￾readable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). unsigned int j,jcycle,jj,jpost,jpre,nf,ng=0,ngrid,nn; double **ires[NGMAX+1],**irho[NGMAX+1],**irhs[NGMAX+1],**iu[NGMAX+1]; nn=n; while (nn >>= 1) ng++; if (n != 1+(1L << ng)) nrerror("n-1 must be a power of 2 in mglin."); if (ng > NGMAX) nrerror("increase NGMAX in mglin."); nn=n/2+1; ngrid=ng-1; irho[ngrid]=dmatrix(1,nn,1,nn); Allocate storage for r.h.s. on grid ng − 1, rstrct(irho[ngrid],u,nn); and fill it by restricting from the fine grid. while (nn > 3) { Similarly allocate storage and fill r.h.s. on all nn=nn/2+1; coarse grids. irho[--ngrid]=dmatrix(1,nn,1,nn); rstrct(irho[ngrid],irho[ngrid+1],nn); } nn=3; iu[1]=dmatrix(1,nn,1,nn); irhs[1]=dmatrix(1,nn,1,nn); slvsml(iu[1],irho[1]); Initial solution on coarsest grid. free_dmatrix(irho[1],1,nn,1,nn); ngrid=ng; for (j=2;j<=ngrid;j++) { Nested iteration loop. nn=2*nn-1; iu[j]=dmatrix(1,nn,1,nn); irhs[j]=dmatrix(1,nn,1,nn); ires[j]=dmatrix(1,nn,1,nn); interp(iu[j],iu[j-1],nn); Interpolate from coarse grid to next finer grid. copy(irhs[j],(j != ngrid ? irho[j] : u),nn); Set up r.h.s. for (jcycle=1;jcycle<=ncycle;jcycle++) { V-cycle loop. nf=nn; for (jj=j;jj>=2;jj--) { Downward stoke of the V. for (jpre=1;jpre<=NPRE;jpre++) Pre-smoothing. relax(iu[jj],irhs[jj],nf); resid(ires[jj],iu[jj],irhs[jj],nf); nf=nf/2+1; rstrct(irhs[jj-1],ires[jj],nf); Restriction of the residual is the next r.h.s. fill0(iu[jj-1],nf); Zero for initial guess in next } relaxation. slvsml(iu[1],irhs[1]); Bottom of V: solve on coars￾nf=3; est grid. for (jj=2;jj<=j;jj++) { Upward stroke of V. nf=2*nf-1; addint(iu[jj],iu[jj-1],ires[jj],nf); Use res for temporary storage inside addint. for (jpost=1;jpost<=NPOST;jpost++) Post-smoothing. relax(iu[jj],irhs[jj],nf); } } } copy(u,iu[ngrid],n); Return solution in u. for (nn=n,j=ng;j>=2;j--,nn=nn/2+1) { free_dmatrix(ires[j],1,nn,1,nn); free_dmatrix(irhs[j],1,nn,1,nn); free_dmatrix(iu[j],1,nn,1,nn); if (j != ng) free_dmatrix(irho[j],1,nn,1,nn); } free_dmatrix(irhs[1],1,3,1,3); free_dmatrix(iu[1],1,3,1,3); }