LAPACK 3.3.1 Linear Algebra PACKage

# zptts2.f

Go to the documentation of this file.
```00001       SUBROUTINE ZPTTS2( IUPLO, N, NRHS, D, E, B, LDB )
00002 *
00003 *  -- LAPACK routine (version 3.3.1) --
00004 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00005 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00006 *  -- April 2011                                                      --
00007 *
00008 *     .. Scalar Arguments ..
00009       INTEGER            IUPLO, LDB, N, NRHS
00010 *     ..
00011 *     .. Array Arguments ..
00012       DOUBLE PRECISION   D( * )
00013       COMPLEX*16         B( LDB, * ), E( * )
00014 *     ..
00015 *
00016 *  Purpose
00017 *  =======
00018 *
00019 *  ZPTTS2 solves a tridiagonal system of the form
00020 *     A * X = B
00021 *  using the factorization A = U**H *D*U or A = L*D*L**H computed by ZPTTRF.
00022 *  D is a diagonal matrix specified in the vector D, U (or L) is a unit
00023 *  bidiagonal matrix whose superdiagonal (subdiagonal) is specified in
00024 *  the vector E, and X and B are N by NRHS matrices.
00025 *
00026 *  Arguments
00027 *  =========
00028 *
00029 *  IUPLO   (input) INTEGER
00030 *          Specifies the form of the factorization and whether the
00031 *          vector E is the superdiagonal of the upper bidiagonal factor
00032 *          U or the subdiagonal of the lower bidiagonal factor L.
00033 *          = 1:  A = U**H *D*U, E is the superdiagonal of U
00034 *          = 0:  A = L*D*L**H, E is the subdiagonal of L
00035 *
00036 *  N       (input) INTEGER
00037 *          The order of the tridiagonal matrix A.  N >= 0.
00038 *
00039 *  NRHS    (input) INTEGER
00040 *          The number of right hand sides, i.e., the number of columns
00041 *          of the matrix B.  NRHS >= 0.
00042 *
00043 *  D       (input) DOUBLE PRECISION array, dimension (N)
00044 *          The n diagonal elements of the diagonal matrix D from the
00045 *          factorization A = U**H *D*U or A = L*D*L**H.
00046 *
00047 *  E       (input) COMPLEX*16 array, dimension (N-1)
00048 *          If IUPLO = 1, the (n-1) superdiagonal elements of the unit
00049 *          bidiagonal factor U from the factorization A = U**H*D*U.
00050 *          If IUPLO = 0, the (n-1) subdiagonal elements of the unit
00051 *          bidiagonal factor L from the factorization A = L*D*L**H.
00052 *
00053 *  B       (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
00054 *          On entry, the right hand side vectors B for the system of
00055 *          linear equations.
00056 *          On exit, the solution vectors, X.
00057 *
00058 *  LDB     (input) INTEGER
00059 *          The leading dimension of the array B.  LDB >= max(1,N).
00060 *
00061 *  =====================================================================
00062 *
00063 *     .. Local Scalars ..
00064       INTEGER            I, J
00065 *     ..
00066 *     .. External Subroutines ..
00067       EXTERNAL           ZDSCAL
00068 *     ..
00069 *     .. Intrinsic Functions ..
00070       INTRINSIC          DCONJG
00071 *     ..
00072 *     .. Executable Statements ..
00073 *
00074 *     Quick return if possible
00075 *
00076       IF( N.LE.1 ) THEN
00077          IF( N.EQ.1 )
00078      \$      CALL ZDSCAL( NRHS, 1.D0 / D( 1 ), B, LDB )
00079          RETURN
00080       END IF
00081 *
00082       IF( IUPLO.EQ.1 ) THEN
00083 *
00084 *        Solve A * X = B using the factorization A = U**H *D*U,
00085 *        overwriting each right hand side vector with its solution.
00086 *
00087          IF( NRHS.LE.2 ) THEN
00088             J = 1
00089    10       CONTINUE
00090 *
00091 *           Solve U**H * x = b.
00092 *
00093             DO 20 I = 2, N
00094                B( I, J ) = B( I, J ) - B( I-1, J )*DCONJG( E( I-1 ) )
00095    20       CONTINUE
00096 *
00097 *           Solve D * U * x = b.
00098 *
00099             DO 30 I = 1, N
00100                B( I, J ) = B( I, J ) / D( I )
00101    30       CONTINUE
00102             DO 40 I = N - 1, 1, -1
00103                B( I, J ) = B( I, J ) - B( I+1, J )*E( I )
00104    40       CONTINUE
00105             IF( J.LT.NRHS ) THEN
00106                J = J + 1
00107                GO TO 10
00108             END IF
00109          ELSE
00110             DO 70 J = 1, NRHS
00111 *
00112 *              Solve U**H * x = b.
00113 *
00114                DO 50 I = 2, N
00115                   B( I, J ) = B( I, J ) - B( I-1, J )*DCONJG( E( I-1 ) )
00116    50          CONTINUE
00117 *
00118 *              Solve D * U * x = b.
00119 *
00120                B( N, J ) = B( N, J ) / D( N )
00121                DO 60 I = N - 1, 1, -1
00122                   B( I, J ) = B( I, J ) / D( I ) - B( I+1, J )*E( I )
00123    60          CONTINUE
00124    70       CONTINUE
00125          END IF
00126       ELSE
00127 *
00128 *        Solve A * X = B using the factorization A = L*D*L**H,
00129 *        overwriting each right hand side vector with its solution.
00130 *
00131          IF( NRHS.LE.2 ) THEN
00132             J = 1
00133    80       CONTINUE
00134 *
00135 *           Solve L * x = b.
00136 *
00137             DO 90 I = 2, N
00138                B( I, J ) = B( I, J ) - B( I-1, J )*E( I-1 )
00139    90       CONTINUE
00140 *
00141 *           Solve D * L**H * x = b.
00142 *
00143             DO 100 I = 1, N
00144                B( I, J ) = B( I, J ) / D( I )
00145   100       CONTINUE
00146             DO 110 I = N - 1, 1, -1
00147                B( I, J ) = B( I, J ) - B( I+1, J )*DCONJG( E( I ) )
00148   110       CONTINUE
00149             IF( J.LT.NRHS ) THEN
00150                J = J + 1
00151                GO TO 80
00152             END IF
00153          ELSE
00154             DO 140 J = 1, NRHS
00155 *
00156 *              Solve L * x = b.
00157 *
00158                DO 120 I = 2, N
00159                   B( I, J ) = B( I, J ) - B( I-1, J )*E( I-1 )
00160   120          CONTINUE
00161 *
00162 *              Solve D * L**H * x = b.
00163 *
00164                B( N, J ) = B( N, J ) / D( N )
00165                DO 130 I = N - 1, 1, -1
00166                   B( I, J ) = B( I, J ) / D( I ) -
00167      \$                        B( I+1, J )*DCONJG( E( I ) )
00168   130          CONTINUE
00169   140       CONTINUE
00170          END IF
00171       END IF
00172 *
00173       RETURN
00174 *
00175 *     End of ZPTTS2
00176 *
00177       END
```