C++: Tridiagonal Solver


by Saladsamurai
Tags: solver, tridiagonal
Saladsamurai
Saladsamurai is online now
#1
Jan21-12, 06:18 PM
Saladsamurai's Avatar
P: 3,012
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.

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?
Phys.Org News Partner Science news on Phys.org
Going nuts? Turkey looks to pistachios to heat new eco-city
Space-tested fluid flow concept advances infectious disease diagnoses
SpaceX launches supplies to space station (Update)
AlephZero
AlephZero is offline
#2
Jan21-12, 07:10 PM
Engineering
Sci Advisor
HW Helper
Thanks
P: 6,344
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
Saladsamurai is online now
#3
Jan21-12, 07:57 PM
Saladsamurai's Avatar
P: 3,012
Quote Quote by AlephZero View Post
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
jtbell is offline
#4
Jan21-12, 08:12 PM
Mentor
jtbell's Avatar
P: 11,231

C++: Tridiagonal Solver


Quote Quote by Saladsamurai View Post
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:

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
Hurkyl is offline
#5
Jan21-12, 08:37 PM
Emeritus
Sci Advisor
PF Gold
Hurkyl's Avatar
P: 16,101
Typo fixes:

    double *returnedArray
should be
    double *&returnedArray
(I hope I have the ordering on * and & correct)

and

    returnedArray = new double(1000);
should be
    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
Saladsamurai is online now
#6
Jan21-12, 08:48 PM
Saladsamurai's Avatar
P: 3,012
Quote Quote by jtbell View Post
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:

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 Quote by Hurkyl View Post
Typo fixes:

    double *returnedArray
should be
    double *&returnedArray
(I hope I have the ordering on * and & correct)

and

    returnedArray = new double(1000);
should be
    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
AlephZero is offline
#7
Jan21-12, 09:32 PM
Engineering
Sci Advisor
HW Helper
Thanks
P: 6,344
Quote Quote by jtbell View Post
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
Saladsamurai is online now
#8
Jan23-12, 08:55 PM
Saladsamurai's Avatar
P: 3,012
Quote Quote by AlephZero View Post
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 like Serena is offline
#9
Jan24-12, 12:19 PM
HW Helper
I like Serena's Avatar
P: 6,189
I would do something like:

#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:
  1. No memory allocations that you have to track yourself.
  2. Preallocating the result to avoid performance overhead.
  3. Return the solution as an output parameter to avoid the performance penalty that returning an object brings.
  4. Return a result that indicates whether the algorithm was successful, otherwise you have no good way to know if there is a solution.
  5. Pass input parameters as "const" to protect yourself from mistakes and to indicate that they are input parameters
  6. Pass all parameters by reference (&) to avoid the performance penalty that copying an object brings.
  7. Define your type for a vector in one place for easy maintainability.
AlephZero
AlephZero is offline
#10
Jan24-12, 08:05 PM
Engineering
Sci Advisor
HW Helper
Thanks
P: 6,344
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...


Register to reply

Related Discussions
symmetric tridiagonal matrix Calculus & Beyond Homework 3
Diagonalisation a tridiagonal symnmnetric matrix General Math 0
eigenvalues of a tridiagonal matrix Programming & Computer Science 4
Splines/Tridiagonal Matrices General Math 0
Tridiagonal matrix Calculus & Beyond Homework 1