C++: Tridiagonal Solver

Saladsamurai

I am writing a Tridiagonal solver using c++. I know that they exist already, but I want to go through the motions for the sake of learning. I am a novice C++ user, however, even before I start writing the program, I get the feeling that the manner in which I am going to write it is not the best way. This is how I was about to start my function, but I want some input on alternative ways to do this.

Code:

void triSolve(double oUp[],double oLow[], double oMain[], 
                                       double oRHS[], int n, double soln[])
{
    // Thomas' Algorithm starts here

    return;
}

where:
oUp[] = upper diagonal of tridiagonal system
oLow[] = lower diagonal of tridiagonal system
oMain[] = Main diagonal of tridiagonal system
n = size of system
oRHS[] = right hand side of tridiagonal system
soln[] = solution to the problem

As you can see, this function returns nothing; it modifies the value of the solution vector in the main() program by replacing its contents with the contents of the formal parameter soln[]. I feel like this isn't the best way. Should I instead try to return the solution vector to the main program?

What are your thoughts?

AlephZero

Just my 2 cents worth. Others may disagree!

1. For efficiency, the less dynamic allocation and freeing of memory the better.
2. For programming sanity, it's rarely a good idea for a function to allocate a significant amount of memory, and then rely on the caller of the function to free it when it's no longer required. The consequences of getting the memory management logic wrong are either memory leaks or crashes.
3. Many library equation solvers return the solution in the same array that you supplied the right hand side. That way, if the caller wants to remember what the right hand side was, it's his/her responsibility - and mostly, you don't want to remember it anyway.

Saladsamurai

Quote by AlephZero

Just my 2 cents worth. Others may disagree!

1. For efficiency, the less dynamic allocation and freeing of memory the better.
2. For programming sanity, it's rarely a good idea for a function to allocate a significant amount of memory, and then rely on the caller of the function to free it when it's no longer required. The consequences of getting the memory management logic wrong are either memory leaks or crashes.
3. Many library equation solvers return the solution in the same array that you supplied the right hand side. That way, if the caller wants to remember what the right hand side was, it's his/her responsibility - and mostly, you don't want to remember it anyway.

Hi AlphaZero

1. For efficiency, the less dynamic allocation and freeing of memory the better.

Thanks for your input.While I understand what you are saying at face value, I am not sure what impact is has on my code as written. How do I avoid 'dynamically allocating memory'? Do you mean, that the less amount of formal parameters that I declare, the better? Or is there more to it than that?

2. For programming sanity, it's rarely a good idea for a function to allocate a significant amount of memory, and then rely on the caller of the function to free it when it's no longer required. The consequences of getting the memory management logic wrong are either memory leaks or crashes.

What does it mean "to rely on the caller to free it"? (I am new at all this.) So if within a function I declare a local object double someArray[1000] and then I call the function from within main(), it allocates enough memory for the array. When it returns from the function, does it not automatically deallocate the memory?

3. Many library equation solvers return the solution in the same array that you supplied the right hand side. That way, if the caller wants to remember what the right hand side was, it's his/her responsibility - and mostly, you don't want to remember it anyway.

I like it. You are right that the RHS is no longer needed within the function and so it could be overwritten.

jtbell

Quote by Saladsamurai

What does it mean "to rely on the caller to free it"? (I am new at all this.) So if within a function I declare a local object double someArray[1000] and then I call the function from within main(), it allocates enough memory for the array. When it returns from the function, does it not automatically deallocate the memory?

It does, and you don't have to worry about it. AlephZero is referring to this sort of thing, where the function allocates the memory for the array dynamically and then returns a pointer to it:

Code:

void Whatever (double *returnedArray /* and other parameters */)
{
    returnedArray = new double(1000);
    /* fill the array */
}

In this case it's up to the calling program to deallocate the memory using free(), in order to avoid a memory leak.

Hurkyl

Typo fixes:

Code:

    double *returnedArray

should be

Code:

    double *&returnedArray

(I hope I have the ordering on * and & correct)

and

Code:

    returnedArray = new double(1000);

should be

Code:

    returnedArray = new double[1000];

and

free()

should be

delete[]

Of course, since nobody's said so yet, I feel like I should point out that you generally shouldn't be using this paradigm in C++; usually one should be using standard containers (or custom containers/views as appropriate) instead of dynamically allocating arrays by hand.

Saladsamurai

Quote by jtbell

It does, and you don't have to worry about it. AlephZero is referring to this sort of thing, where the function allocates the memory for the array dynamically and then returns a pointer to it:

Code:

void Whatever (double *returnedArray /* and other parameters */)
{
    returnedArray = new double(1000);
    /* fill the array */
}

In this case it's up to the calling program to deallocate the memory using free(), in order to avoid a memory leak.

Quote by Hurkyl

Typo fixes:

Code:

    double *returnedArray

should be

Code:

    double *&returnedArray

(I hope I have the ordering on * and & correct)

and

Code:

    returnedArray = new double(1000);

should be

Code:

    returnedArray = new double[1000];

and

free()

should be

delete[]

Of course, since nobody's said so yet, I feel like I should point out that you generally shouldn't be using this paradigm in C++; usually one should be using standard containers (or custom containers/views as appropriate) instead of dynamically allocating arrays by hand.

Thanks for the responses folks

Like I said, I am pretty new to this. So I still have to learn about things like pointers, containers, structs, objects ...

Anyone know of any good free references where I can learn about c++ programming for engineers & scientists? I have been looking for something that assumes little prior knowledge. It can be fast paced, but I just don't want something that assumes I am an expert programmer already (like this NIST Course which is too advanced for me).

AlephZero

Quote by jtbell

It does, and you don't have to worry about it. AlephZero is referring to this sort of thing, where the function allocates the memory for the array dynamically and then returns a pointer to it:

Sure, creating local objects on the stack isn't an issue (it's an essential part of implementing the language, and the programmer doesn't have any control over it, so it's GOT to be done efficiently).

But explicitly or implicitly creating objects on the heap can get very expensive if you overdose on it.

I agree about using standard library container objects (though they don't add much functionality here) - but the same memory access issues apply either way, even if the containers make them less obviously visible.

"Premature optimisation is the root of all evil" is also a good motto - but it's a fair bet that something like an equation solving routine will be worth optimising. Letting a container class expand itself incrementally to the right size, when you already know what the right size is, just being lazy IMO.

Saladsamurai

Quote by AlephZero

Sure, creating local objects on the stack isn't an issue (it's an essential part of implementing the language, and the programmer doesn't have any control over it, so it's GOT to be done efficiently).

But explicitly or implicitly creating objects on the heap can get very expensive if you overdose on it.

I agree about using standard library container objects (though they don't add much functionality here) - but the same memory access issues apply either way, even if the containers make them less obviously visible.

"Premature optimisation is the root of all evil" is also a good motto - but it's a fair bet that something like an equation solving routine will be worth optimising. Letting a container class expand itself incrementally to the right size, when you already know what the right size is, just being lazy IMO.

Hi AlphaZero

Sorry, but is there any way for you to dumb this down a bit? In particular, is there something about the way I am writing this function, that you would do differently in addition to simply returning the solution vector to oRHS? I am wide open to ideas as I am here to learn. Thanks you.

I like Serena

I would do something like:

Code:

#include <vector>

typedef std::vector<double> MyVector;

bool triSolve(const MyVector& oUp, const MyVector& oLow, const MyVector& oMain, 
              const MyVector& oRHS, MyVector& soln)
{
    soln.size(oUp.size());
    bool success = false;

    // Thomas' Algorithm starts here

    return success;
}

Highlights:

No memory allocations that you have to track yourself.
Preallocating the result to avoid performance overhead.
Return the solution as an output parameter to avoid the performance penalty that returning an object brings.
Return a result that indicates whether the algorithm was successful, otherwise you have no good way to know if there is a solution.
Pass input parameters as "const" to protect yourself from mistakes and to indicate that they are input parameters
Pass all parameters by reference (&) to avoid the performance penalty that copying an object brings.
Define your type for a vector in one place for easy maintainability.

AlephZero

What ILS said. You might also think about making a class to represent a "tridiagonal matrix" object. Passing 3 separate parameters that have to be correctly related to each other is longwinded to write, and error prone, compared with passing just one. For example you might get the order of the 3 vectors mixed up, or accidentally create them with the wrong lengths, or whatever. And in future you might want to work with both symmetric and non-symmetric tridiagonal matrices, and/or a matrix of complex numbers instead of reals...

Similar Threads for: C++: Tridiagonal Solver
Thread	Forum	Replies
symmetric tridiagonal matrix	Calculus & Beyond Homework	3
Diagonalisation a tridiagonal symnmnetric matrix	General Math	0
eigenvalues of a tridiagonal matrix	Programming & Comp Sci	4
Splines/Tridiagonal Matrices	General Math	0
Tridiagonal matrix	Calculus & Beyond Homework	1