Oh my god... sorry another typo, ill just post the whole code that isnt typo ridden, those 2 loops dont really use the same variable.
Code:
/****************************************************************
* Sequential program to solve Ax=b using Jacobi's method *
****************************************************************
* A - Compact Storage Vector
* C - Column vector for A
* R - Row position vector for A
* x(n) - solution vector *
* b(n) - RHS vector *
****************************************************************/
#include <stdio.h>
#include <math.h>
#include <time.h>
#define n 4000
main()
{
int iter,maxit,rcheck,s;
int i,j;
double xnew[n],xold[n];
double b[n],a[n*n],aii,sum;
int r[n+1],c[(n*n)],l,m,k;
double toler,diff,maxdif;
clock_t t1;
toler = .000001;
iter=0;
maxit=50;
maxdif=1.0;
s=1;
rcheck=1;
/* generate sample data for compact storage vector A and right- hand-side-vector (RHS) b */
for (i = 0;i<=n-1;i++){
rcheck = 1;
for(j=0;j<=n-1;j++){
if (i==j) {
a[s] = 4.0;
c[s] = j;
s=s+1;
if (rcheck==1){
r[i] = s-1;
rcheck = 0;
}
}
if (i==j-1||i==j+1){
a[s] = -1.0;
c[s] = j;
s=s+1;
if (rcheck==1){
r[i] = s-1;
rcheck = 0;
}
}
}
b[i] = 2.0;
}
/*FOR TESTING THE COMPACT STORAGE SCHEME*/
/* for(i=1;i<=s;i++){
printf("\n%4.2f %d %d\n",a[i],c[i],r[i]);
}
*/
b[0] = 3.0;
b[n-1] = 3.0;
/* The N+1 element in R must be set to the size of A+1 to give the sum loop
work to do on the last iteration*/
r[n] = s+1;
/* time section of code to solve Ax = b*/
t1 = clock();
//set initial solution vector to RHS Vector b
for (i = 0;i<=n-1;i++){
xold[i] = b[i];
}
// iterate until |xold - xnew| < toler
// or max number of iterations exceeded
// R(I) represents the first non zero in a row
// so it works to sum 1 row at a time, C(J) is the
// corresponding column that element was in in the matrix
while ((iter<maxit)&&(maxdif>=toler)){
for (i=0;i<=n-1;i++){
sum = b[i];
for (j = r[i];j<=r[i+1]-1;j++){
if (c[j]!=i) sum = sum + (-a[j]*xold[c[j]]);
}
//Finds the A(I,I) element in the Compact Storage Vector
m=r[i+1]-1;
l=r[i];
k=l;
if (c[k]==i) aii=a[k];
while ((c[k]!=i)&&(k<=m)){
k=k+1;
if (c[k]==i) aii=a[k];
}
xnew[i] = sum / aii;
}
// determine maximum vector norm
maxdif = fabs(xold[0]-xnew[0]);
for (i=1;i<=n-1;i++){
diff = fabs(xold[i]-xnew[i]);
if (diff>maxdif) maxdif = diff;
}
// reset recent approximation to the solution vector x to the newly generated x values
for(i=0;i<=n-1;i++){
xold[i]=xnew[i];
}
// keep track of the number of iterations required to reach the given tolerance
iter=iter+1;
}
printf("\nElapsed time = %f",clock()-t1);
printf("\nSystem size = %d x %d",n,n);
printf("\nNumber of iterations = %d",iter);
printf("\nMaximum Error = %f\n",maxdif);
// used for testing solution
// printf("\n");
// for(i=0;i<=n-1;i++){
// printf("%f\n",xold[i]);
// }
}