There are probably a hundred things old fashioned or inefficient about this C code. But mainly, I just want to know why I am getting a seg fault in the routine "jacobi_solve". The driver code is here:
#include <stdio.h>
#include <stdlib.h>
int main() {
int m=20, n=20;
int i,j;
double A[m][n],Anew[m][n];
for (j=0;j<n;j++){
for (i=0;i<m;i++){
A[j][i] = 0.; Anew[j][i] = 0.;
}
}
A[10][10] = 1.;
printf("%x\n",A);
jacobi_solver(A,Anew,m,n);
}
and jacobi_solver is here:
#include<stdio.h>
#include<math.h>
void jacobi_solver(double **A, double **Anew, int m, int n) {
double err,tol=1.e-6;
int iter=0,iter_max=100;
int i,j;
err = tol*10.;
printf("hello\n");
printf("%f %f %d %d\n",err,tol,iter,iter_max);
printf("%x\n",A);
printf("%f\n",A[10][10]); /* <--- seg fault! */
printf("solving ...\n");
while ( err > tol && iter < iter_max ) {
for( j = 1; j < n-1; j++) {
for( i = 1; i < m-1; i++) {
Anew[j][i] = 0.25 * (A[j][i+1] + A[j][i-1] +
A[j-1][i] + A[j+1][i]);
err = fmax(err, abs(Anew[j][i] - A[j][i]));
}
}
printf("%f\n", err);
for( j = 1; j < n-1; j++) {
for( i = 1; i < m-1; i++ ) {
A[j][i] = Anew[j][i];
}
}
iter++;
}
}
Note that if I take away the print statement, I get a seg fault later, in the main solving loop.