19.6 Multigrid Methods for Boundary Value Problems 881 u[2][2]=-h*h*rhs[2][2]/4.0; void relax(double *u,double **rhs,int n) Red-black Gauss-Seidel relaxation for model problem.Updates the current value of the solution u[1..n][1..n],using the right-hand side function rhs [1..n][1..n]. int i,ipass,isw,j,jsw=1; double h,h2; h=1.0/(n-1); h2=h*h; for (ipass=1;ipass<=2;ipass++,jsw=3-jsw){ Red and black sweeps. isw=isw; 83 for(j2;j<知；j+,isw=3-1sw) granted for 19881992 for(1=isw+1;1<n;1+=2) Gauss-Seidel formula. u[i][j]=0.25*(u[1+1][j]+u[i-1][j]+u[i][+1] (including this one) +u[i][j-1]-h2*rhs[i][j]); 11-800 872 to any by Cambridge NUMERICAL RECIPES IN void resid(double **res,double *u,double **rhs,int n) Returns mins the residual for the model problem.Input quantities are u[1..n][1..n]and (Nort rhs[1..n][1..n],while res[1..n][1..n]is returned. f America server computer, users to make one paper e University Press. THE int i,j; double h,h2i; 是 ART h=1.0/(n-1); Programs h2i=1.0/(h*h); for (j=2;j<n;j++) Interior points for(1=2;1<n;1++) st st res[1][j]=-h2i*(u[i+1][]+u[i-1][j]+u[1][j+1]+u[1][j-1]- 4.0*[i][j])+rhs[i][j]; for (i=1;i<=n;i++) Boundary points. res[i][1]=res[i][n]=res[1][i]=res [n][i]=0.0; void copy(double **aout,double *ain,int n) Copies ain [1..n][1..n]to aout [1..n][1..n]. int i,j; for(1=1;1<=n;1+) email to directcustserv@cambridge.org(outside North America). 1988-1992 by Numerical Recipes OF SCIENTIFIC COMPUTING(ISBN 0-521-43108-5) for (j=1;j<=n;++) aout [j][i]=ain[j][i]; Software. void fill0(double *u,int n) Fills u[1..n][1..n]with zeros. int i,j; ving of machine for (j=1;j<=n;++) for(1=1;i<=n;1+) u[i][j]=0.0: The routine mglin is written for clarity,not maximum efficiency,so that it is easy to modify.Several simple changes will speed up the execution time:19.6 Multigrid Methods for Boundary Value Problems 881 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). u[2][2] = -h*h*rhs[2][2]/4.0; } void relax(double **u, double **rhs, int n) Red-black Gauss-Seidel relaxation for model problem. Updates the current value of the solution u[1..n][1..n], using the right-hand side function rhs[1..n][1..n]. { int i,ipass,isw,j,jsw=1; double h,h2; h=1.0/(n-1); h2=h*h; for (ipass=1;ipass<=2;ipass++,jsw=3-jsw) { Red and black sweeps. isw=jsw; for (j=2;j<n;j++,isw=3-isw) for (i=isw+1;i<n;i+=2) Gauss-Seidel formula. u[i][j]=0.25*(u[i+1][j]+u[i-1][j]+u[i][j+1] +u[i][j-1]-h2*rhs[i][j]); } } void resid(double **res, double **u, double **rhs, int n) Returns minus the residual for the model problem. Input quantities are u[1..n][1..n] and rhs[1..n][1..n], while res[1..n][1..n] is returned. { int i,j; double h,h2i; h=1.0/(n-1); h2i=1.0/(h*h); for (j=2;j<n;j++) Interior points. for (i=2;i<n;i++) res[i][j] = -h2i*(u[i+1][j]+u[i-1][j]+u[i][j+1]+u[i][j-1]- 4.0*u[i][j])+rhs[i][j]; for (i=1;i<=n;i++) Boundary points. res[i][1]=res[i][n]=res[1][i]=res[n][i]=0.0; } void copy(double **aout, double **ain, int n) Copies ain[1..n][1..n] to aout[1..n][1..n]. { int i,j; for (i=1;i<=n;i++) for (j=1;j<=n;j++) aout[j][i]=ain[j][i]; } void fill0(double **u, int n) Fills u[1..n][1..n] with zeros. { int i,j; for (j=1;j<=n;j++) for (i=1;i<=n;i++) u[i][j]=0.0; } The routine mglin is written for clarity, not maximum efficiency, so that it is easy to modify. Several simple changes will speed up the execution time: