serial version of Jacobi's method for solving * the linear system Ax = b. *

/* serial_jacobi.c -- serial version of Jacobi's method for solving
* the linear system Ax = b.
*
* Input:
* n: order of system
* tol: convergence tolerance
* max_iter: maximum number of iterations
* A: coefficient matrix
* b: right-hand side of system
*
* Output:
* x: the solution if the method converges
* max_iter: if the method fails to converge
*
* Notes:
* 1. A should be strongly diagonally dominant in
* order to insure convergence.
* 2. A, x, and b are statically allocated.
*
* See Chap 10, pp. 218 & ff in PPMPI.
*/
#include
#include

#define MAX_DIM 12

typedef float MATRIX_T[MAX_DIM][MAX_DIM];

int Jacobi(MATRIX_T A, float x[], float b[], int n,
float tol, int max_iter);
void Read_matrix(char* prompt, MATRIX_T A, int n);
void Read_vector(char* prompt, float x[], int n);
void Print_vector(char* title, float x[], int n);
void Print_matrix(char* title, MATRIX_T A, int n);

main(int argc, char* argv) {
MATRIX_T A;
float x[MAX_DIM];
float b[MAX_DIM];
int n;
float tol;
int max_iter;
int converged;

printf("Enter n, tolerance, and max number of iterations\n");
scanf("%d %f %d", &n, &tol, &max_iter);
Read_matrix("Enter the matrix", A, n);
Read_vector("Enter the right-hand side", b, n);

converged = Jacobi(A, x, b, n, tol, max_iter);

if (converged)
Print_vector("The solution is", x, n);
else
printf("Failed to converge in %d iterations\n", max_iter);
} /* main */


/*********************************************************************/
/* Return 1 if iteration converged, 0 otherwise */
/* MATRIX_T is just a 2-dimensional array */
int Jacobi(
MATRIX_T A /* in */,
float x[] /* out */,
float b[] /* in */,
int n /* in */,
float tol /* in */,
int max_iter /* in */) {
int i, j;
int iter_num;
float x_old[MAX_DIM];

float Distance(float x[], float y[], int n);

/* Initialize x */
for (i = 0; i < n; i++)
x[i] = b[i];

iter_num = 0;
do {
iter_num++;

for (i = 0; i < n; i++)
x_old[i] = x[i];

for (i = 0; i < n; i++){
x[i] = b[i];
for (j = 0; j < i; j++)
x[i] = x[i] - A[i][j]*x_old[j];
for (j = i+1; j < n; j++)
x[i] = x[i] - A[i][j]*x_old[j];
x[i] = x[i]/A[i][i];
}
} while ((iter_num < max_iter) &&
(Distance(x,x_old,n) >= tol));

if (Distance(x,x_old,n) < tol)
return 1;
else
return 0;
} /* Jacobi */


/*********************************************************************/
float Distance(float x[], float y[], int n) {
int i;
float sum = 0.0;

for (i = 0; i < n; i++) {
sum = sum + (x[i] - y[i])*(x[i] - y[i]);
}
return sqrt(sum);
} /* Distance */


/*********************************************************************/
void Read_matrix(
char* prompt /* in */,
MATRIX_T A /* out */,
int n /* in */) {
int i, j;

printf("%s\n", prompt);
for (i = 0; i < n; i++)
for (j = 0; j < n; j++)
scanf("%f", &A[i][j]);
} /* Read_matrix */


/*********************************************************************/
void Read_vector(
char* prompt /* in */,
float x[] /* out */,
int n /* in */) {
int i;

printf("%s\n", prompt);
for (i = 0; i < n; i++)
scanf("%f", &x[i]);
} /* Read_vector */


/*********************************************************************/
void Print_vector(
char* title /* in */,
float x[] /* in */,
int n /* in */) {
int i;

printf("%s\n", title);
for (i = 0; i < n; i++)
printf("%4.1f ", x[i]);
printf("\n");
} /* Print_vector */

/*********************************************************************/
void Print_matrix(
char* title /* in */,
MATRIX_T A /* in */,
int n /* in */) {
int i,j;

printf("%s\n", title);
for (i = 0; i < n; i++) {
for (j = 0; j < n; j++)
printf("%f", A[i][j]);
printf("\n");
}
} /* Print_matrix */