LAPACK 3.3.1
Linear Algebra PACKage

cpftrs.f

Go to the documentation of this file.
00001       SUBROUTINE CPFTRS( TRANSR, UPLO, N, NRHS, A, B, LDB, INFO )
00002 *
00003 *  -- LAPACK routine (version 3.3.1)                                    --
00004 *
00005 *  -- Contributed by Fred Gustavson of the IBM Watson Research Center --
00006 *  -- April 2011                                                      --
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       COMPLEX            A( 0: * ), B( LDB, * )
00017 *     ..
00018 *
00019 *  Purpose
00020 *  =======
00021 *
00022 *  CPFTRS solves a system of linear equations A*X = B with a Hermitian
00023 *  positive definite matrix A using the Cholesky factorization
00024 *  A = U**H*U or A = L*L**H computed by CPFTRF.
00025 *
00026 *  Arguments
00027 *  =========
00028 *
00029 *  TRANSR    (input) CHARACTER*1
00030 *          = 'N':  The Normal TRANSR of RFP A is stored;
00031 *          = 'C':  The Conjugate-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) COMPLEX array, dimension ( N*(N+1)/2 );
00045 *          The triangular factor U or L from the Cholesky factorization
00046 *          of RFP A = U**H*U or RFP A = L*L**H, as computed by CPFTRF.
00047 *          See note below for more details about RFP A.
00048 *
00049 *  B       (input/output) COMPLEX 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 Standard Packed Format when N is even.
00064 *  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 *  conjugate-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 *  conjugate-transpose of the last three columns of AP lower.
00083 *  To denote conjugate we place -- above the element. This covers the
00084 *  case N even and TRANSR = 'N'.
00085 *
00086 *         RFP A                   RFP A
00087 *
00088 *                                -- -- --
00089 *        03 04 05                33 43 53
00090 *                                   -- --
00091 *        13 14 15                00 44 54
00092 *                                      --
00093 *        23 24 25                10 11 55
00094 *
00095 *        33 34 35                20 21 22
00096 *        --
00097 *        00 44 45                30 31 32
00098 *        -- --
00099 *        01 11 55                40 41 42
00100 *        -- -- --
00101 *        02 12 22                50 51 52
00102 *
00103 *  Now let TRANSR = 'C'. RFP A in both UPLO cases is just the conjugate-
00104 *  transpose of RFP A above. One therefore gets:
00105 *
00106 *
00107 *           RFP A                   RFP A
00108 *
00109 *     -- -- -- --                -- -- -- -- -- --
00110 *     03 13 23 33 00 01 02    33 00 10 20 30 40 50
00111 *     -- -- -- -- --                -- -- -- -- --
00112 *     04 14 24 34 44 11 12    43 44 11 21 31 41 51
00113 *     -- -- -- -- -- --                -- -- -- --
00114 *     05 15 25 35 45 55 22    53 54 55 22 32 42 52
00115 *
00116 *
00117 *  We next  consider Standard Packed Format when N is odd.
00118 *  We give an example where N = 5.
00119 *
00120 *     AP is Upper                 AP is Lower
00121 *
00122 *   00 01 02 03 04              00
00123 *      11 12 13 14              10 11
00124 *         22 23 24              20 21 22
00125 *            33 34              30 31 32 33
00126 *               44              40 41 42 43 44
00127 *
00128 *
00129 *  Let TRANSR = 'N'. RFP holds AP as follows:
00130 *  For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last
00131 *  three columns of AP upper. The lower triangle A(3:4,0:1) consists of
00132 *  conjugate-transpose of the first two   columns of AP upper.
00133 *  For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first
00134 *  three columns of AP lower. The upper triangle A(0:1,1:2) consists of
00135 *  conjugate-transpose of the last two   columns of AP lower.
00136 *  To denote conjugate we place -- above the element. This covers the
00137 *  case N odd  and TRANSR = 'N'.
00138 *
00139 *         RFP A                   RFP A
00140 *
00141 *                                   -- --
00142 *        02 03 04                00 33 43
00143 *                                      --
00144 *        12 13 14                10 11 44
00145 *
00146 *        22 23 24                20 21 22
00147 *        --
00148 *        00 33 34                30 31 32
00149 *        -- --
00150 *        01 11 44                40 41 42
00151 *
00152 *  Now let TRANSR = 'C'. RFP A in both UPLO cases is just the conjugate-
00153 *  transpose of RFP A above. One therefore gets:
00154 *
00155 *
00156 *           RFP A                   RFP A
00157 *
00158 *     -- -- --                   -- -- -- -- -- --
00159 *     02 12 22 00 01             00 10 20 30 40 50
00160 *     -- -- -- --                   -- -- -- -- --
00161 *     03 13 23 33 11             33 11 21 31 41 51
00162 *     -- -- -- -- --                   -- -- -- --
00163 *     04 14 24 34 44             43 44 22 32 42 52
00164 *
00165 *  =====================================================================
00166 *
00167 *     .. Parameters ..
00168       COMPLEX            CONE
00169       PARAMETER          ( CONE = ( 1.0E+0, 0.0E+0 ) )
00170 *     ..
00171 *     .. Local Scalars ..
00172       LOGICAL            LOWER, NORMALTRANSR
00173 *     ..
00174 *     .. External Functions ..
00175       LOGICAL            LSAME
00176       EXTERNAL           LSAME
00177 *     ..
00178 *     .. External Subroutines ..
00179       EXTERNAL           XERBLA, CTFSM
00180 *     ..
00181 *     .. Intrinsic Functions ..
00182       INTRINSIC          MAX
00183 *     ..
00184 *     .. Executable Statements ..
00185 *
00186 *     Test the input parameters.
00187 *
00188       INFO = 0
00189       NORMALTRANSR = LSAME( TRANSR, 'N' )
00190       LOWER = LSAME( UPLO, 'L' )
00191       IF( .NOT.NORMALTRANSR .AND. .NOT.LSAME( TRANSR, 'C' ) ) THEN
00192          INFO = -1
00193       ELSE IF( .NOT.LOWER .AND. .NOT.LSAME( UPLO, 'U' ) ) THEN
00194          INFO = -2
00195       ELSE IF( N.LT.0 ) THEN
00196          INFO = -3
00197       ELSE IF( NRHS.LT.0 ) THEN
00198          INFO = -4
00199       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
00200          INFO = -7
00201       END IF
00202       IF( INFO.NE.0 ) THEN
00203          CALL XERBLA( 'CPFTRS', -INFO )
00204          RETURN
00205       END IF
00206 *
00207 *     Quick return if possible
00208 *
00209       IF( N.EQ.0 .OR. NRHS.EQ.0 )
00210      $   RETURN
00211 *
00212 *     start execution: there are two triangular solves
00213 *
00214       IF( LOWER ) THEN
00215          CALL CTFSM( TRANSR, 'L', UPLO, 'N', 'N', N, NRHS, CONE, A, B,
00216      $               LDB )
00217          CALL CTFSM( TRANSR, 'L', UPLO, 'C', 'N', N, NRHS, CONE, A, B,
00218      $               LDB )
00219       ELSE
00220          CALL CTFSM( TRANSR, 'L', UPLO, 'C', 'N', N, NRHS, CONE, A, B,
00221      $               LDB )
00222          CALL CTFSM( TRANSR, 'L', UPLO, 'N', 'N', N, NRHS, CONE, A, B,
00223      $               LDB )
00224       END IF
00225 *
00226       RETURN
00227 *
00228 *     End of CPFTRS
00229 *
00230       END
 All Files Functions