Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

C++: Tridiagonal Solver

  1. Jan 21, 2012 #1
    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 (Text):

    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?
     
  2. jcsd
  3. Jan 21, 2012 #2

    AlephZero

    User Avatar
    Science Advisor
    Homework Helper

    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.
     
  4. Jan 21, 2012 #3
    Hi AlphaZero :smile:

    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?

    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?

    I like it. You are right that the RHS is no longer needed within the function and so it could be overwritten.
     
  5. Jan 21, 2012 #4

    jtbell

    User Avatar

    Staff: Mentor

    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 (Text):

    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.
     
  6. Jan 21, 2012 #5

    Hurkyl

    User Avatar
    Staff Emeritus
    Science Advisor
    Gold Member

    Typo fixes:

    Code (Text):

        double *returnedArray
     
    should be
    Code (Text):

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

    and

    Code (Text):

        returnedArray = new double(1000);
     
    should be
    Code (Text):

        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.
     
  7. Jan 21, 2012 #6
    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).
     
  8. Jan 21, 2012 #7

    AlephZero

    User Avatar
    Science Advisor
    Homework Helper

    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.
     
  9. Jan 23, 2012 #8
    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.
     
  10. Jan 24, 2012 #9

    I like Serena

    User Avatar
    Homework Helper

    I would do something like:

    Code (Text):
    #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.
     
  11. Jan 24, 2012 #10

    AlephZero

    User Avatar
    Science Advisor
    Homework Helper

    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...
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook




Similar Discussions: C++: Tridiagonal Solver
Loading...