1       SUBROUTINE DPPTRS( UPLO, N, NRHS, AP, B, LDB, INFO )
  2 *
  3 *  -- LAPACK routine (version 3.3.1) --
  4 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  5 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  6 *  -- April 2011                                                      --
  7 *
  8 *     .. Scalar Arguments ..
  9       CHARACTER          UPLO
 10       INTEGER            INFO, LDB, N, NRHS
 11 *     ..
 12 *     .. Array Arguments ..
 13       DOUBLE PRECISION   AP( * ), B( LDB, * )
 14 *     ..
 15 *
 16 *  Purpose
 17 *  =======
 18 *
 19 *  DPPTRS solves a system of linear equations A*X = B with a symmetric
 20 *  positive definite matrix A in packed storage using the Cholesky
 21 *  factorization A = U**T*U or A = L*L**T computed by DPPTRF.
 22 *
 23 *  Arguments
 24 *  =========
 25 *
 26 *  UPLO    (input) CHARACTER*1
 27 *          = 'U':  Upper triangle of A is stored;
 28 *          = 'L':  Lower triangle of A is stored.
 29 *
 30 *  N       (input) INTEGER
 31 *          The order of the matrix A.  N >= 0.
 32 *
 33 *  NRHS    (input) INTEGER
 34 *          The number of right hand sides, i.e., the number of columns
 35 *          of the matrix B.  NRHS >= 0.
 36 *
 37 *  AP      (input) DOUBLE PRECISION array, dimension (N*(N+1)/2)
 38 *          The triangular factor U or L from the Cholesky factorization
 39 *          A = U**T*U or A = L*L**T, packed columnwise in a linear
 40 *          array.  The j-th column of U or L is stored in the array AP
 41 *          as follows:
 42 *          if UPLO = 'U', AP(i + (j-1)*j/2) = U(i,j) for 1<=i<=j;
 43 *          if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = L(i,j) for j<=i<=n.
 44 *
 45 *  B       (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
 46 *          On entry, the right hand side matrix B.
 47 *          On exit, the solution matrix X.
 48 *
 49 *  LDB     (input) INTEGER
 50 *          The leading dimension of the array B.  LDB >= max(1,N).
 51 *
 52 *  INFO    (output) INTEGER
 53 *          = 0:  successful exit
 54 *          < 0:  if INFO = -i, the i-th argument had an illegal value
 55 *
 56 *  =====================================================================
 57 *
 58 *     .. Local Scalars ..
 59       LOGICAL            UPPER
 60       INTEGER            I
 61 *     ..
 62 *     .. External Functions ..
 63       LOGICAL            LSAME
 64       EXTERNAL           LSAME
 65 *     ..
 66 *     .. External Subroutines ..
 67       EXTERNAL           DTPSV, XERBLA
 68 *     ..
 69 *     .. Intrinsic Functions ..
 70       INTRINSIC          MAX
 71 *     ..
 72 *     .. Executable Statements ..
 73 *
 74 *     Test the input parameters.
 75 *
 76       INFO = 0
 77       UPPER = LSAME( UPLO, 'U' )
 78       IF.NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
 79          INFO = -1
 80       ELSE IF( N.LT.0 ) THEN
 81          INFO = -2
 82       ELSE IF( NRHS.LT.0 ) THEN
 83          INFO = -3
 84       ELSE IF( LDB.LT.MAX1, N ) ) THEN
 85          INFO = -6
 86       END IF
 87       IF( INFO.NE.0 ) THEN
 88          CALL XERBLA( 'DPPTRS'-INFO )
 89          RETURN
 90       END IF
 91 *
 92 *     Quick return if possible
 93 *
 94       IF( N.EQ.0 .OR. NRHS.EQ.0 )
 95      $   RETURN
 96 *
 97       IF( UPPER ) THEN
 98 *
 99 *        Solve A*X = B where A = U**T * U.
100 *
101          DO 10 I = 1, NRHS
102 *
103 *           Solve U**T *X = B, overwriting B with X.
104 *
105             CALL DTPSV( 'Upper''Transpose''Non-unit', N, AP,
106      $                  B( 1, I ), 1 )
107 *
108 *           Solve U*X = B, overwriting B with X.
109 *
110             CALL DTPSV( 'Upper''No transpose''Non-unit', N, AP,
111      $                  B( 1, I ), 1 )
112    10    CONTINUE
113       ELSE
114 *
115 *        Solve A*X = B where A = L * L**T.
116 *
117          DO 20 I = 1, NRHS
118 *
119 *           Solve L*Y = B, overwriting B with X.
120 *
121             CALL DTPSV( 'Lower''No transpose''Non-unit', N, AP,
122      $                  B( 1, I ), 1 )
123 *
124 *           Solve L**T *X = Y, overwriting B with X.
125 *
126             CALL DTPSV( 'Lower''Transpose''Non-unit', N, AP,
127      $                  B( 1, I ), 1 )
128    20    CONTINUE
129       END IF
130 *
131       RETURN
132 *
133 *     End of DPPTRS
134 *
135       END