#include <math.h> #define MAXITS 1000 #define EPS 1.0e-5 void sor(double **a, double **b, double **c, double **d, double **e, double **f, double **u, int jmax, double rjac) { void nrerror(char error_text[]); int ipass,j,jsw,l,lsw,n; double anorm,anormf=0.0,omega=1.0,resid; for (j=2;j<jmax;j++) for (l=2;l<jmax;l++) anormf += fabs(f[j][l]); for (n=1;n<=MAXITS;n++) { anorm=0.0; jsw=1; for (ipass=1;ipass<=2;ipass++) { lsw=jsw; for (j=2;j<jmax;j++) { for (l=lsw+1;l<jmax;l+=2) { resid=a[j][l]*u[j+1][l] +b[j][l]*u[j-1][l] +c[j][l]*u[j][l+1] +d[j][l]*u[j][l-1] +e[j][l]*u[j][l] -f[j][l]; anorm += fabs(resid); u[j][l] -= omega*resid/e[j][l]; } lsw=3-lsw; } jsw=3-jsw; omega=(n == 1 && ipass == 1 ? 1.0/(1.0-0.5*rjac*rjac) : 1.0/(1.0-0.25*rjac*rjac*omega)); } if (anorm < EPS*anormf) return; } nrerror("MAXITS exceeded"); } #undef MAXITS #undef EPS