```01:       SUBROUTINE DPTTS2( N, NRHS, D, E, B, LDB )
02: *
03: *  -- LAPACK routine (version 3.2) --
04: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
05: *     November 2006
06: *
07: *     .. Scalar Arguments ..
08:       INTEGER            LDB, N, NRHS
09: *     ..
10: *     .. Array Arguments ..
11:       DOUBLE PRECISION   B( LDB, * ), D( * ), E( * )
12: *     ..
13: *
14: *  Purpose
15: *  =======
16: *
17: *  DPTTS2 solves a tridiagonal system of the form
18: *     A * X = B
19: *  using the L*D*L' factorization of A computed by DPTTRF.  D is a
20: *  diagonal matrix specified in the vector D, L is a unit bidiagonal
21: *  matrix whose subdiagonal is specified in the vector E, and X and B
22: *  are N by NRHS matrices.
23: *
24: *  Arguments
25: *  =========
26: *
27: *  N       (input) INTEGER
28: *          The order of the tridiagonal matrix A.  N >= 0.
29: *
30: *  NRHS    (input) INTEGER
31: *          The number of right hand sides, i.e., the number of columns
32: *          of the matrix B.  NRHS >= 0.
33: *
34: *  D       (input) DOUBLE PRECISION array, dimension (N)
35: *          The n diagonal elements of the diagonal matrix D from the
36: *          L*D*L' factorization of A.
37: *
38: *  E       (input) DOUBLE PRECISION array, dimension (N-1)
39: *          The (n-1) subdiagonal elements of the unit bidiagonal factor
40: *          L from the L*D*L' factorization of A.  E can also be regarded
41: *          as the superdiagonal of the unit bidiagonal factor U from the
42: *          factorization A = U'*D*U.
43: *
44: *  B       (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
45: *          On entry, the right hand side vectors B for the system of
46: *          linear equations.
47: *          On exit, the solution vectors, X.
48: *
49: *  LDB     (input) INTEGER
50: *          The leading dimension of the array B.  LDB >= max(1,N).
51: *
52: *  =====================================================================
53: *
54: *     .. Local Scalars ..
55:       INTEGER            I, J
56: *     ..
57: *     .. External Subroutines ..
58:       EXTERNAL           DSCAL
59: *     ..
60: *     .. Executable Statements ..
61: *
62: *     Quick return if possible
63: *
64:       IF( N.LE.1 ) THEN
65:          IF( N.EQ.1 )
66:      \$      CALL DSCAL( NRHS, 1.D0 / D( 1 ), B, LDB )
67:          RETURN
68:       END IF
69: *
70: *     Solve A * X = B using the factorization A = L*D*L',
71: *     overwriting each right hand side vector with its solution.
72: *
73:       DO 30 J = 1, NRHS
74: *
75: *           Solve L * x = b.
76: *
77:          DO 10 I = 2, N
78:             B( I, J ) = B( I, J ) - B( I-1, J )*E( I-1 )
79:    10    CONTINUE
80: *
81: *           Solve D * L' * x = b.
82: *
83:          B( N, J ) = B( N, J ) / D( N )
84:          DO 20 I = N - 1, 1, -1
85:             B( I, J ) = B( I, J ) / D( I ) - B( I+1, J )*E( I )
86:    20    CONTINUE
87:    30 CONTINUE
88: *
89:       RETURN
90: *
91: *     End of DPTTS2
92: *
93:       END
94: ```