Implementing Jacobi and Gauss Seidel methods (format / structure)

[cswiercz]

From a student's private repo:

So I'm having some real trouble understanding where to start in implementing the jacobi and gauss_seidel methods. I have tried to make new methods within the solvers.c file, which include the decompose and jacobi_step method and then call these methods within the jacobi method. However, I'm not so sure as to how we are supposed to decompose our matrices. Are we using the python functions we created in homework1 or are we writing new ones from scratch. I understand that we can use the methods we have written in linalg.c to matrix multiply and matrix vector multiply and so on... but I'm confused on where to start within implementing these algorithms.

I had originally thought about returning by reference the D, L and U matrices if I made a method called decompose, and work from there, but I'm not sure if I need to dynamically allocate memory and also I'm having some trouble forming the D, L and U matrices. Is there a nice way to do this, like in python.

One last thing, I'm still extremely confused on pointers and their size, especially when using them within methods and passing them in as parameters. Would you mind putting up any helpful links.resources on this?

[cswiercz]

So I'm having some real trouble understanding where to start in implementing the jacobi and gauss_seidel methods. I have tried to make new methods within the solvers.c file, which include the decompose and jacobi_step method and then call these methods within the jacobi method.

First of all, let me say that there are many, many ways to implement Jacobi and Gauss-Seidel in C. Because you have low-level access to the data there are tricks you can use to minimize the number of additional memory allocations you have to make for storing temporary data.

Several students have taken their Python code from Homework 1 and rewrote it in C, which is a good start. So suppose you wanted to write a decompose function similar to Homework 1. My first question to you is,

what is the function prototype of decompose?

That is, what are the arguments? Should it return things by value or by reference? Where should you allocate memory for the matrices, D, L, and U? Do you need to create space for all of them?

However, I'm not so sure as to how we are supposed to decompose our matrices. Are we using the python functions we created in homework1 or are we writing new ones from scratch.

Refer to Homework 1 for the definitions of Jacobi and Gauss-Seidel methods. Although it is possible(?) to call Python functions from C it is...difficult. So you will definitely need to write a completely C version of the linear system solvers.

I understand that we can use the methods we have written in linalg.c to matrix multiply and matrix vector multiply and so on... but I'm confused on where to start within implementing these algorithms.

The definitions for these functions should be supplied in linalg.c. A good starting point is to read the (already provided to you) code written in vec_add. The function prototypes for these are given in linalg.h which you can import into solvers.c. See the lecture on multi-file C programs for a refresher on header files if this part is confusing. Once you've imported the functionality of linalg.h into solvers.c you can use all of that functionality in your Jacobi and Gauss-Seidel code, if you feel like it's necessary. (Hint: most likely.)

Any additional "helper" methods you need to write you can write directly into solvers.c. For example, if you want to use decompose function to "extract" D, L, and U you can write the definition of this function before jacobi(). Remember that this definition should occur before where it's used or else you'll get an error, possibly a warning depending on which compiler you're using.

Finally, take a look at the Makefile to see how the linalg.c and solvers.c code are compiled. See the lecture on multi-file C programs for a refresher on shared-object libraries.

[cswiercz]

I had originally thought about returning by reference the D, L and U matrices if I made a method called decompose, and work from there, but I'm not sure if I need to dynamically allocate memory

I would recommend returning these by reference. Where to allocate these arrays is a very good question. Do you allocate them on the stack? If so, the stack location of which function? jacobi? jacobi_step? decompose? Or, do you allocate them on the heap? If so, within which function? jacobi? jacobi_step? decompose? Some choices are better than others.

If you heap-allocate your temporary storage then remember to free at the end?

also I'm having some trouble forming the D, L and U matrices. Is there a nice way to do this, like in python.

What do you mean by "nice"? In Python, Numpy provides many useful helper functions. We don't have that in C so we would have to write them ourselves. But once we write them we can use them all over the place!

[cswiercz]

One last thing, I'm still extremely confused on pointers and their size, especially when using them within methods and passing them in as parameters. Would you mind putting up any helpful links.resources on this?

I have been writing many useful comments in the lecture code examples. Start there.