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 /////////////////////////////////////////