LAPACK 3.3.0

dpftrs.f

Go to the documentation of this file.
00001       SUBROUTINE DPFTRS( TRANSR, UPLO, N, NRHS, A, B, LDB, INFO )
00002 *
00003 *  -- LAPACK routine (version 3.3.0)                                    --
00004 *
00005 *  -- Contributed by Fred Gustavson of the IBM Watson Research Center --
00006 *     November 2010
00007 *
00008 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00009 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00010 *
00011 *     .. Scalar Arguments ..
00012       CHARACTER          TRANSR, UPLO
00013       INTEGER            INFO, LDB, N, NRHS
00014 *     ..
00015 *     .. Array Arguments ..
00016       DOUBLE PRECISION   A( 0: * ), B( LDB, * )
00017 *     ..
00018 *
00019 *  Purpose
00020 *  =======
00021 *
00022 *  DPFTRS solves a system of linear equations A*X = B with a symmetric
00023 *  positive definite matrix A using the Cholesky factorization
00024 *  A = U**T*U or A = L*L**T computed by DPFTRF.
00025 *
00026 *  Arguments
00027 *  =========
00028 *
00029 *  TRANSR  (input) CHARACTER*1
00030 *          = 'N':  The Normal TRANSR of RFP A is stored;
00031 *          = 'T':  The Transpose TRANSR of RFP A is stored.
00032 *
00033 *  UPLO    (input) CHARACTER*1
00034 *          = 'U':  Upper triangle of RFP A is stored;
00035 *          = 'L':  Lower triangle of RFP A is stored.
00036 *
00037 *  N       (input) INTEGER
00038 *          The order of the matrix A.  N >= 0.
00039 *
00040 *  NRHS    (input) INTEGER
00041 *          The number of right hand sides, i.e., the number of columns
00042 *          of the matrix B.  NRHS >= 0.
00043 *
00044 *  A       (input) DOUBLE PRECISION array, dimension ( N*(N+1)/2 ).
00045 *          The triangular factor U or L from the Cholesky factorization
00046 *          of RFP A = U**T*U or RFP A = L*L**T, as computed by DPFTRF.
00047 *          See note below for more details about RFP A.
00048 *
00049 *  B       (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
00050 *          On entry, the right hand side matrix B.
00051 *          On exit, the solution matrix X.
00052 *
00053 *  LDB     (input) INTEGER
00054 *          The leading dimension of the array B.  LDB >= max(1,N).
00055 *
00056 *  INFO    (output) INTEGER
00057 *          = 0:  successful exit
00058 *          < 0:  if INFO = -i, the i-th argument had an illegal value
00059 *
00060 *  Further Details
00061 *  ===============
00062 *
00063 *  We first consider Rectangular Full Packed (RFP) Format when N is
00064 *  even. We give an example where N = 6.
00065 *
00066 *      AP is Upper             AP is Lower
00067 *
00068 *   00 01 02 03 04 05       00
00069 *      11 12 13 14 15       10 11
00070 *         22 23 24 25       20 21 22
00071 *            33 34 35       30 31 32 33
00072 *               44 45       40 41 42 43 44
00073 *                  55       50 51 52 53 54 55
00074 *
00075 *
00076 *  Let TRANSR = 'N'. RFP holds AP as follows:
00077 *  For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last
00078 *  three columns of AP upper. The lower triangle A(4:6,0:2) consists of
00079 *  the transpose of the first three columns of AP upper.
00080 *  For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first
00081 *  three columns of AP lower. The upper triangle A(0:2,0:2) consists of
00082 *  the transpose of the last three columns of AP lower.
00083 *  This covers the case N even and TRANSR = 'N'.
00084 *
00085 *         RFP A                   RFP A
00086 *
00087 *        03 04 05                33 43 53
00088 *        13 14 15                00 44 54
00089 *        23 24 25                10 11 55
00090 *        33 34 35                20 21 22
00091 *        00 44 45                30 31 32
00092 *        01 11 55                40 41 42
00093 *        02 12 22                50 51 52
00094 *
00095 *  Now let TRANSR = 'T'. RFP A in both UPLO cases is just the
00096 *  transpose of RFP A above. One therefore gets:
00097 *
00098 *
00099 *           RFP A                   RFP A
00100 *
00101 *     03 13 23 33 00 01 02    33 00 10 20 30 40 50
00102 *     04 14 24 34 44 11 12    43 44 11 21 31 41 51
00103 *     05 15 25 35 45 55 22    53 54 55 22 32 42 52
00104 *
00105 *
00106 *  We then consider Rectangular Full Packed (RFP) Format when N is
00107 *  odd. We give an example where N = 5.
00108 *
00109 *     AP is Upper                 AP is Lower
00110 *
00111 *   00 01 02 03 04              00
00112 *      11 12 13 14              10 11
00113 *         22 23 24              20 21 22
00114 *            33 34              30 31 32 33
00115 *               44              40 41 42 43 44
00116 *
00117 *
00118 *  Let TRANSR = 'N'. RFP holds AP as follows:
00119 *  For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last
00120 *  three columns of AP upper. The lower triangle A(3:4,0:1) consists of
00121 *  the transpose of the first two columns of AP upper.
00122 *  For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first
00123 *  three columns of AP lower. The upper triangle A(0:1,1:2) consists of
00124 *  the transpose of the last two columns of AP lower.
00125 *  This covers the case N odd and TRANSR = 'N'.
00126 *
00127 *         RFP A                   RFP A
00128 *
00129 *        02 03 04                00 33 43
00130 *        12 13 14                10 11 44
00131 *        22 23 24                20 21 22
00132 *        00 33 34                30 31 32
00133 *        01 11 44                40 41 42
00134 *
00135 *  Now let TRANSR = 'T'. RFP A in both UPLO cases is just the
00136 *  transpose of RFP A above. One therefore gets:
00137 *
00138 *           RFP A                   RFP A
00139 *
00140 *     02 12 22 00 01             00 10 20 30 40 50
00141 *     03 13 23 33 11             33 11 21 31 41 51
00142 *     04 14 24 34 44             43 44 22 32 42 52
00143 *
00144 *  =====================================================================
00145 *
00146 *     .. Parameters ..
00147       DOUBLE PRECISION   ONE
00148       PARAMETER          ( ONE = 1.0D+0 )
00149 *     ..
00150 *     .. Local Scalars ..
00151       LOGICAL            LOWER, NORMALTRANSR
00152 *     ..
00153 *     .. External Functions ..
00154       LOGICAL            LSAME
00155       EXTERNAL           LSAME
00156 *     ..
00157 *     .. External Subroutines ..
00158       EXTERNAL           XERBLA, DTFSM
00159 *     ..
00160 *     .. Intrinsic Functions ..
00161       INTRINSIC          MAX
00162 *     ..
00163 *     .. Executable Statements ..
00164 *
00165 *     Test the input parameters.
00166 *
00167       INFO = 0
00168       NORMALTRANSR = LSAME( TRANSR, 'N' )
00169       LOWER = LSAME( UPLO, 'L' )
00170       IF( .NOT.NORMALTRANSR .AND. .NOT.LSAME( TRANSR, 'T' ) ) THEN
00171          INFO = -1
00172       ELSE IF( .NOT.LOWER .AND. .NOT.LSAME( UPLO, 'U' ) ) THEN
00173          INFO = -2
00174       ELSE IF( N.LT.0 ) THEN
00175          INFO = -3
00176       ELSE IF( NRHS.LT.0 ) THEN
00177          INFO = -4
00178       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
00179          INFO = -7
00180       END IF
00181       IF( INFO.NE.0 ) THEN
00182          CALL XERBLA( 'DPFTRS', -INFO )
00183          RETURN
00184       END IF
00185 *
00186 *     Quick return if possible
00187 *
00188       IF( N.EQ.0 .OR. NRHS.EQ.0 )
00189      +   RETURN
00190 *
00191 *     start execution: there are two triangular solves
00192 *
00193       IF( LOWER ) THEN
00194          CALL DTFSM( TRANSR, 'L', UPLO, 'N', 'N', N, NRHS, ONE, A, B,
00195      +               LDB )
00196          CALL DTFSM( TRANSR, 'L', UPLO, 'T', 'N', N, NRHS, ONE, A, B,
00197      +               LDB )
00198       ELSE
00199          CALL DTFSM( TRANSR, 'L', UPLO, 'T', 'N', N, NRHS, ONE, A, B,
00200      +               LDB )
00201          CALL DTFSM( TRANSR, 'L', UPLO, 'N', 'N', N, NRHS, ONE, A, B,
00202      +               LDB )
00203       END IF
00204 *
00205       RETURN
00206 *
00207 *     End of DPFTRS
00208 *
00209       END
 All Files Functions