Is there a better way to write a Tridiagonal solver in C++?

  • C/C++
  • Thread starter Saladsamurai
  • Start date
In summary: in which case the caller of the function is responsible for freeing the memory when they're done with it.
  • #1
Saladsamurai
3,020
7
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?
 
Technology news on Phys.org
  • #2
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.
 
  • #3
AlephZero said:
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 :smile:

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.
 
  • #4
Saladsamurai said:
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.
 
  • #5
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.
 
  • #6
jtbell said:
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 said:
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 :smile: 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).
 
  • #7
jtbell said:
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.
 
  • #8
AlephZero said:
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 :smile: 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.
 
  • #9
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:
  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.
 
  • #10
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...
 

1. What is a tridiagonal matrix in C++?

A tridiagonal matrix is a square matrix where all the elements outside of the main diagonal and the two adjacent diagonals are zero. It is commonly used in numerical methods and linear algebra problems.

2. How does a tridiagonal solver work in C++?

A tridiagonal solver is a specialized algorithm used to efficiently solve tridiagonal systems of linear equations. It uses a combination of back-substitution and forward elimination to find the solution.

3. What are the advantages of using a tridiagonal solver in C++?

Some advantages of using a tridiagonal solver in C++ include its efficiency, accuracy, and stability. It also requires less memory and computational resources compared to other methods for solving linear equations.

4. Can a tridiagonal solver be used for non-tridiagonal matrices in C++?

No, a tridiagonal solver is specifically designed to solve tridiagonal systems of linear equations. It may not produce accurate results for non-tridiagonal matrices, and other methods should be used instead.

5. How can I implement a tridiagonal solver in my C++ code?

There are various ways to implement a tridiagonal solver in C++, including using built-in libraries or writing your own algorithm. It is recommended to do some research and understand the steps involved in solving a tridiagonal system before attempting to implement it in your code.

Similar threads

  • Programming and Computer Science
Replies
5
Views
2K
  • Programming and Computer Science
Replies
23
Views
2K
  • Programming and Computer Science
Replies
6
Views
8K
  • Programming and Computer Science
Replies
2
Views
1K
  • Programming and Computer Science
2
Replies
35
Views
2K
  • Programming and Computer Science
Replies
23
Views
1K
  • Programming and Computer Science
Replies
2
Views
4K
  • Programming and Computer Science
Replies
17
Views
1K
  • Programming and Computer Science
Replies
1
Views
2K
  • Programming and Computer Science
Replies
1
Views
898
Back
Top