Hi All,
Please excuse any newbee errors. I recently purchased 'Parallel Iterative Algorithms' (From Sequential to Grid Computing) by Bahi et al. I am trying to implement Algorithm 2.1 on page 16 (Jacobi Method) in C.
I have a basic working knowledge of C but am struggling a little. The code below is an initial stab at it. Its not working fully yet but what I would really like is some comment and help if possible to get it working.
My primary motive first and foremost is to learn C (for some future gpu related work) and to be able to implement code from algorithms rather than just working with examples in books.
The algorithm is given below and my code attempt to date follows. Many thanks in advance for your help.
Algorithm 2.1 Jacobi algorithm
Size: size of the matrix
X[Size]: solution vector
XOld[Size]: solution vector at the previous iteration
A[Size][Size]: matrix
B[Size]: right-hand side vector
repeat
for i=0 to Size−1 do
X[i] ← 0
for j=0 to i−1 do
X[i] ← X[i]+A[i][j]×XOld[j]
end for
for j=i+1 to Size−1 do
X[i] ← X[i]+A[i][j]×XOld[j]
end for
end for
for i=0 to Size−1 do
XOld[i] ← (B[i]−X[i])/A[i][i]
end for
until stopping criteria is reached
Code:
#include <stdio.h>
#include <math.h>
#include <iostream>
#define SIZE 3
/* define function prototypes */
void jacobi(float *A_ptr, float *b_ptr, float *x_ptr);
void display_A(float *A_ptr);
void diagonal_A(float *A_ptr);
void swap_x(float *xOld, float *xNew);
/* function main begins program execution */
int main(int argc, char *argv[])
{
/* array definitions */
float A[][3] = {
{6, -2, 1},
{-2,7,2},
{1,2,-5}
};
float b[] = {11,2,-1};
float x[] = {0,0,0};
/* pointer definitions */
float *ptr_A = (float *)&A;
//ptr_A = (float *)&A;
printf("The base address of of Matrix A i.e. A[0]: %p\n\n", ptr_A); /* check */
printf("The base address of ptr_A: %p\n\n", &ptr_A);
float *ptr_b;
ptr_b = (float *)&b;
float *ptr_x;
ptr_x = (float *)&x;
printf("\n"); // start a new line of output
printf("Base address of A before function call: %p \n", ptr_A); // base address of A before calling function
printf("\nJacobi Iteration - Start:\n");
jacobi(ptr_A, ptr_b, ptr_x);
printf("\nJacobi Iteration - Complete:\n");
system("PAUSE");
return EXIT_SUCCESS; // indicate successful termination
}
void jacobi(float *A_ptr, float *b_ptr, float *x_ptr)
{
const int max_iter = 10; // max number of iterations
int i, j, iter;
float xOld[SIZE];
float xNew[SIZE];
float *xOld_ptr;
xOld_ptr = (float *)&xOld;
float *xNew_ptr;
xNew_ptr = (float *)&xNew;
// iteration counter
for(iter = 1; iter <= max_iter; ++iter)
{
swap_x(xOld_ptr, xNew_ptr);
// loop through rows
for (i = 0; i <= SIZE-1; ++i)
{
//xNew_ptr[i] = 0; // method 1: initialising the solution vector
*(xNew_ptr + i) = 0; // method 2: initialising the solution vector
// loop through columns
printf("\nHello 1");
for(j = 0; j <= i-1; ++j) // Lower diagonal??
{
printf("Hello");
xNew_ptr[i] += (*(A_ptr + i * SIZE + j) * xOld_ptr[j]);
}
printf("\nHello 2");
for(j = i+1; j = SIZE-1; ++j) // Upper diagonal??
{
xNew_ptr[i] += (*(A_ptr + i * SIZE + j) * xOld_ptr[j]);
}
}
printf("\nHello 3");
for(i = 0; i = SIZE-1; ++i)
{
for(j = 0; j <= SIZE-1; ++j)
{
xOld_ptr[i] = (b_ptr[i] - xNew_ptr[i]) / *(A_ptr + i * SIZE + i);
}
}
}
printf("Maximun number of iterations reached\n");
//return ;
}
//////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////// Update xOld & xNew ///////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////////////
void swap_x(float *xOld_ptr, float *xNew_ptr)
{
int i;
for(i = 0; i <= SIZE-1; ++i)
{
xOld_ptr[i] = xNew_ptr[i];
}
}
///////////////////////////////////////// End of Function /////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////// Display the element values of matrix A //////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////////////
void display_A(float *A_ptr)
{
int i; /* row counter */
int j; /* column counter */
printf("\nFunction Name: display_A\n");
printf("\nBase address of A at function call: %p \n", A_ptr); // base address of A at function call
// check to see pointer to matrix A is correct
printf("\nThe coefficient Matrix A is:\n");
for(i = 0; i <= SIZE-1; ++i)
{
for (j = 0; j <= SIZE-1; ++j)
{
printf("%2.0f ", *(A_ptr + i * SIZE + j));
/* *(base address + row no. * no. of columns + column no. */
}
printf("\n"); // start a new line of output
}
printf("\nEnd of Function: display_A\n");
}
///////////////////////////////////////// End of Function /////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////// Display the diagonal elements of matrix A ////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////////////////
void diagonal_A(float *A_ptr)
{
int i; /* row counter */
int j; /* column counter */
float diag[SIZE];
float *diag_ptr = (float *)&diag;
printf("\nFunction Name: diagonal_A\n");
printf("\nBase address of A at function call: %p \n", A_ptr); // base address of A at function call
// print the diagonal elements of the matrix
printf("\nThe diagonal elements of matrix A are:\n");
for(i = 0; i <= SIZE-1; ++i)
{
for(j = 0;j <= SIZE-1; ++j)
{
if(i == j)
{
diag[i] = *(A_ptr + i * SIZE + j);
printf("%2.0f \n", *(A_ptr + i * SIZE + j));
}
}
}
printf("\nEnd of Function: diagonal_A\n");
}
///////////////////////////////////////// End of Function /////////////////////////////////////////