Like Tree1Likes
  • 1 Post By tabstop

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. #1
    Registered User
    Join Date
    Jul 2011
    Location
    Ireland
    Posts
    5

    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. #2
    Banned
    Join Date
    Aug 2010
    Location
    Ontario Canada
    Posts
    9,547
    So...Was there a question in there someplace?

  3. #3
    and the Hat of Guessing tabstop's Avatar
    Join Date
    Nov 2007
    Posts
    14,185
    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.
    AndrewHunter likes this.

  4. #4
    Registered User
    Join Date
    Jul 2011
    Location
    Ireland
    Posts
    5
    Quote Originally Posted by CommonTater View Post
    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. #5
    Registered User
    Join Date
    Jul 2011
    Location
    Ireland
    Posts
    5
    Quote Originally Posted by tabstop View Post
    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. #6
    Banned
    Join Date
    Aug 2010
    Location
    Ontario Canada
    Posts
    9,547
    Quote Originally Posted by roubalard View Post
    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. #7
    and the Hat of Guessing tabstop's Avatar
    Join Date
    Nov 2007
    Posts
    14,185
    Quote Originally Posted by roubalard View Post
    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. #8
    Registered User
    Join Date
    Sep 2006
    Posts
    8,868
    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. #9
    Registered User
    Join Date
    Jul 2011
    Location
    Ireland
    Posts
    5
    Quote Originally Posted by tabstop View Post
    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. #10
    Registered User
    Join Date
    Jul 2011
    Location
    Ireland
    Posts
    5
    Quote Originally Posted by Adak View Post
    ** You DID review what tabstop said about missing the final assignment, right? **
    Thanks. I'm just a bit confused. Can you expand please.

  11. #11
    C++ Witch laserlight's Avatar
    Join Date
    Oct 2003
    Location
    Singapore
    Posts
    21,660
    Quote 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.
    C + C++ Compiler: MinGW port of GCC
    Version Control System: Bazaar

    Look up a C++ Reference and learn How To Ask Questions The Smart Way

  12. #12
    Algorithm Dissector iMalc's Avatar
    Join Date
    Dec 2005
    Location
    New Zealand
    Posts
    6,302
    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.
    My homepage
    Advice: Take only as directed - If symptoms persist, please see your debugger

    Linus Torvalds: "But it clearly is the only right way. The fact that everybody else does it some other way only means that they are wrong"

  13. #13
    Registered User
    Join Date
    Sep 2006
    Posts
    8,868
    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 subscribe to a feed

Similar Threads

  1. Implementing Library Algorithms
    By jewelz in forum C++ Programming
    Replies: 9
    Last Post: 04-06-2009, 12:44 AM
  2. Iterative functions
    By russel1013 in forum C Programming
    Replies: 7
    Last Post: 07-21-2008, 03:14 AM
  3. iterative combinations
    By antreas_z in forum C++ Programming
    Replies: 3
    Last Post: 06-02-2003, 06:39 AM
  4. iterative, recursion, always possible?
    By Vber in forum C Programming
    Replies: 7
    Last Post: 02-28-2003, 12:59 AM

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21