I implemented an $LU$ factorization algorithm on c++ for a $n$ by $n$ matrix $A$, and I wanted to now write a routine that solves the systems $Lv = f$ and $Uv = f$. How would you recommend writing such a routine. I am aware of using a linear algebra package on c++ in which I could take the inverse of $L$ and $U$ and solve for $v$ to make the expressions equal to. Although, I was wondering if there is a simpler way of solving a system like this in c++. Any suggestions is greatly appreciated. I made a post to this here

  • $\begingroup$ Check out en.wikipedia.org/wiki/LU_decomposition and en.wikipedia.org/wiki/… for the standard way to solve linear systems given a LU factorization. $\endgroup$ Feb 22 '16 at 22:40
  • $\begingroup$ Also, you should be solving $L y = f$, and then solving $U x = y$ in order to solve the system $(LU)x = f$. $\endgroup$ Feb 22 '16 at 22:45
  • $\begingroup$ yea I did that, but I fail to see that by solving $Ly = f$ and $U x = y$ then I can solve the system $Lv = f$ and $Uv = f$ $\endgroup$
    – Wolfy
    Feb 22 '16 at 23:31
  • $\begingroup$ The system that you present doesn't make a lot of sense in the context of an LU factorization if you're trying to find a $v$ that satisfies both of those simultaneously (assuming that $f$ is given by the problem and $L$, $U$ obtained by factorization). The systems $Lv=f$ and $Uv=f$ are independent of one another and can be solved by either forward-substitution or backward-substitution. $\endgroup$ Feb 22 '16 at 23:57

The basic strategy here is to take advantage of the known structure of the two matrices (lower and upper triangular respectively). This means that the inversion can be found efficiently through forward/backward substitution. For example, the operation of the first row of $L$ on $v$ gives $L_{00}x_0=f_0$, which (for the cases you've been told to consider in which pivoting is unnecessary) can be solved through simple scalar division. Given your knowledge of $x_0$, you can then generate $\tilde{f}_i=f_i-L_{i0}x_0$ for $i\in [1,N-1]$. At the end of this process you have N-1 unknowns, an N-1 dimensional right hand side and an $N-1\times N-1$ lower triangular matrix to work with. This recursion forms the outer loop of the solve. The key is to store $L$ and $U$ such that the substitution process accesses contiguous static memory.

If you haven't already, I strongly recommend reading the relevant subsection of Numerical Recipes by Press et al. The latest edition is even in C++. I would also suggest adding the homework tag, as you have in your previous question.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.