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

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

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.

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.

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++'.

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.)

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.

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? **

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.

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.

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".