# Implementing iterative algorithms in C

This is a discussion on Implementing iterative algorithms in C within the C Programming forums, part of the General Programming Boards category; Hi All, Please excuse any newbee errors. I recently purchased 'Parallel Iterative Algorithms' (From Sequential to Grid Computing) by Bahi ...

1. ## Implementing iterative algorithms in C

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

2. So...Was there a question in there someplace?

3. It's customary to write for loops like
Code:
`for (i = 0; i < SIZE; i++)`
although what you have is a pretty direct translation of your algorithm, which is also good. I don't know what you have against [bracket notation], but you should go to a doctor for it.

You never put your answer into your x array at the end of Jacobi.

4. Originally Posted by CommonTater
So...Was there a question in there someplace?
Thanks for the reply. One thing I was getting confused with were pointers. So some possible questions might be; Do mine make sense? Am I using them correctly? I don't want to get into a habit of using pointers without fully understanding when and where they should be used. Eventually I want to be able to write clean efficient code.

5. Originally Posted by tabstop
It's customary to write for loops like
Code:
`for (i = 0; i < SIZE; i++)`
although what you have is a pretty direct translation of your algorithm, which is also good. I don't know what you have against [bracket notation], but you should go to a doctor for it.

You never put your answer into your x array at the end of Jacobi.
Thanks for reply. I'm just getting to grips with c. The code itself is a work in progress. Is there any significant difference between using '++i' and 'i++'.

6. Originally Posted by roubalard
Thanks for the reply. One thing I was getting confused with were pointers. So some possible questions might be; Do mine make sense? Am I using them correctly? I don't want to get into a habit of using pointers without fully understanding when and where they should be used. Eventually I want to be able to write clean efficient code.
I have a simple rule about pointers: Don't use them unless you have to.

I use the array[i] notation almost exclusively for arrays. I use the struct.var notation in structs.

It's that old saw... "Just because you CAN do something doesn't mean you SHOULD"

Pointers themselves are pretty simple... they're like post office boxes... the pointer is the box number and you need to dereference (unlock) it to get at the contents... int *p = pointer.... &addr = address.... *p = contents. (For some patently silly reason the C gods decided to make the creation and the use syntactically idential... causing lord only knows how much confusion over the years.)

7. Originally Posted by roubalard
Thanks for reply. I'm just getting to grips with c. The code itself is a work in progress. Is there any significant difference between using '++i' and 'i++'.
No, and I apologize for switching it because that's not the difference you are supposed to focus on. Look in the middle.

8. Yes, there is an efficiency difference in most C compilers, between using pointers, and using the array[index] notation. Of course, some cpu's are so optimized to work with integers that the index notation is actually faster.

For the rest, you'll need to pull one hair out of your head. Split it lengthwise, four to 10 times. Now hold it up to the light, and you can clearly see the difference pointer notation makes, for yourself, in efficiency.

** You DID review what tabstop said about missing the final assignment, right? **

9. Originally Posted by tabstop
No, and I apologize for switching it because that's not the difference you are supposed to focus on. Look in the middle.
Thanks for pointing that out. It's clicked now what you getting at. That's one of the reasons why I decided to post - there's nothing like a fresh pair of eyes.

10. Originally Posted by Adak
** You DID review what tabstop said about missing the final assignment, right? **
Thanks. I'm just a bit confused. Can you expand please.

11. Originally Posted by Adak
Yes, there is an efficiency difference in most C compilers, between using pointers, and using the array[index] notation. Of course, some cpu's are so optimized to work with integers that the index notation is actually faster.
The small efficiency difference only comes into play if you are incrementing the pointer rather than computing with an offset on each iteration (i.e., instead of both incrementing the index and computing the resulting pointer, you just increment the pointer). Besides, optimising compilers can often recognise this and produce code as if you wrote the loop in a more efficient manner by incrementing the pointer.

12. Use the arrays directly in this instance. Copying their start address into pointers and then using he pointers instead means that the compiler will not be able to work out the ailasing rules for that piece of code and will have to assume that a write through one pointer could affect a read through another. As a result, using pointers cannot result in as efficient code.

Really for this code, there's no reason that a pointer should appear anywhere at all. All the things you are doing with pointers can be done with the arrays themselves instead.

13. I'm just echoing tabstop's observation, here.

You never put your answer into your x array at the end of Jacobi
With rare exceptions, I think programmers who prefer to use pointers to work with arrays, have lost their guide on the trail, and are now swerving back to where the snow is 3 feet deep.

To learn about arrays, and how pointers can work with them; for simple walking through an array, say checking for palindromes or reversing strings -- great. For general work with arrays, *absolutely* not.

You've gained NOTHING for efficiency, 99.9% of the time, (may have lost some if you aren't careful), and you've given up so much in readibility and clarity.

Why?

I hope it's not just so you can look "clever".

Popular pages Recent additions