44 lines
948 B
C
44 lines
948 B
C
|
|
#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
|