Internal: How to Use LA_TRISOL to Solve Multiple Systems of Linear Equations

The documentation for LA_TRISOL states that "The LA_TRISOL function may also be used to solve for multiple systems of linear equations, with each column of B representing a different set of equations", but it is not clear how this should be used in practice. Where do the k left-hand-sides of the equations go to, i.e. the k A matrices in AX=B? There is an example code in the docs of LA_TRIMPROVE, but this is only for one set of equations (k=1).

LA_TRISOL will solve multiple systems of linear equations involving the same tridiagonal array A. For instance, suppose I have the following two systems of linear equations:

-4t + u = 6
2t - 4u + v = -8
2u - 4v + w = -5
2v -4w = 8

-4t + u = 1
2t - 4u + v = -2
2u - 4v + w = -3
2v -4w = -8

These can be represented as Ax=B1 and Ax=B2 where
         
   |-4 1 0 0|
A = | 2 -4 1 0|
   | 0 2 -4 1|
   | 0 0 2 -4|

B1 =    | 6|
         |-8|
         |-5|
         | 8|

B2 =    | 1|
         |-2|
         |-3|
         |-8|

Both systems can be solved simultaneously using the equation Ax=B where

         | 6 1|
B =       |-8 -2|
         |-5 -3|
         | 8 -8|


Here is an example that uses LA_TRISOL to solve both of these systems at once:

pro test

; Solve the two following tridiagonal system of equations:
; -4t + u = 6
; 2t - 4u + v = -8
; 2u - 4v + w = -5
; 2v -4w = 8
;
; -4t + u = 1
; 2t - 4u + v = -2
; 2u - 4v + w = -3
; 2v -4w = -8


; Define array A:
aUpper = [1, 1, 1]
adiag = [-4, -4, -4, -4]
aLower = [2, 2, 2]

; Define right-hand side vector for first system of equations
b1 = transpose([6, -8, -5, 8])

; Define right-hand side vector for second system of equations
b2 = transpose([1, -2, -3, -8])

; Decompose A:
dLower = aLower
dArray = adiag
dUpper = aUpper
la_tridc, dLower, dArray, dUpper, u2, index

; Compute the solution for the first system of equations
x = la_trisol(dLower, dArray, dUpper, u2, index, b1)
print, 'Solution for first system:'
print, x

; Compute the solution for the second system of equations
x = la_trisol(dLower, dArray, dUpper, u2, index, b2)
print, 'Solution for first system:'
print, x

; Solve both systems at once...

; Define B
b = [b1,b2]
x = la_trisol(dLower, dArray, dUpper, u2, index, b)
print, 'Solution for both:'
print, x

End