LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
zsytrs.f
Go to the documentation of this file.
00001 *> \brief \b ZSYTRS
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download ZSYTRS + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zsytrs.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zsytrs.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zsytrs.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE ZSYTRS( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, INFO )
00022 * 
00023 *       .. Scalar Arguments ..
00024 *       CHARACTER          UPLO
00025 *       INTEGER            INFO, LDA, LDB, N, NRHS
00026 *       ..
00027 *       .. Array Arguments ..
00028 *       INTEGER            IPIV( * )
00029 *       COMPLEX*16         A( LDA, * ), B( LDB, * )
00030 *       ..
00031 *  
00032 *
00033 *> \par Purpose:
00034 *  =============
00035 *>
00036 *> \verbatim
00037 *>
00038 *> ZSYTRS solves a system of linear equations A*X = B with a complex
00039 *> symmetric matrix A using the factorization A = U*D*U**T or
00040 *> A = L*D*L**T computed by ZSYTRF.
00041 *> \endverbatim
00042 *
00043 *  Arguments:
00044 *  ==========
00045 *
00046 *> \param[in] UPLO
00047 *> \verbatim
00048 *>          UPLO is CHARACTER*1
00049 *>          Specifies whether the details of the factorization are stored
00050 *>          as an upper or lower triangular matrix.
00051 *>          = 'U':  Upper triangular, form is A = U*D*U**T;
00052 *>          = 'L':  Lower triangular, form is A = L*D*L**T.
00053 *> \endverbatim
00054 *>
00055 *> \param[in] N
00056 *> \verbatim
00057 *>          N is INTEGER
00058 *>          The order of the matrix A.  N >= 0.
00059 *> \endverbatim
00060 *>
00061 *> \param[in] NRHS
00062 *> \verbatim
00063 *>          NRHS is INTEGER
00064 *>          The number of right hand sides, i.e., the number of columns
00065 *>          of the matrix B.  NRHS >= 0.
00066 *> \endverbatim
00067 *>
00068 *> \param[in] A
00069 *> \verbatim
00070 *>          A is COMPLEX*16 array, dimension (LDA,N)
00071 *>          The block diagonal matrix D and the multipliers used to
00072 *>          obtain the factor U or L as computed by ZSYTRF.
00073 *> \endverbatim
00074 *>
00075 *> \param[in] LDA
00076 *> \verbatim
00077 *>          LDA is INTEGER
00078 *>          The leading dimension of the array A.  LDA >= max(1,N).
00079 *> \endverbatim
00080 *>
00081 *> \param[in] IPIV
00082 *> \verbatim
00083 *>          IPIV is INTEGER array, dimension (N)
00084 *>          Details of the interchanges and the block structure of D
00085 *>          as determined by ZSYTRF.
00086 *> \endverbatim
00087 *>
00088 *> \param[in,out] B
00089 *> \verbatim
00090 *>          B is COMPLEX*16 array, dimension (LDB,NRHS)
00091 *>          On entry, the right hand side matrix B.
00092 *>          On exit, the solution matrix X.
00093 *> \endverbatim
00094 *>
00095 *> \param[in] LDB
00096 *> \verbatim
00097 *>          LDB is INTEGER
00098 *>          The leading dimension of the array B.  LDB >= max(1,N).
00099 *> \endverbatim
00100 *>
00101 *> \param[out] INFO
00102 *> \verbatim
00103 *>          INFO is INTEGER
00104 *>          = 0:  successful exit
00105 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
00106 *> \endverbatim
00107 *
00108 *  Authors:
00109 *  ========
00110 *
00111 *> \author Univ. of Tennessee 
00112 *> \author Univ. of California Berkeley 
00113 *> \author Univ. of Colorado Denver 
00114 *> \author NAG Ltd. 
00115 *
00116 *> \date November 2011
00117 *
00118 *> \ingroup complex16SYcomputational
00119 *
00120 *  =====================================================================
00121       SUBROUTINE ZSYTRS( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, INFO )
00122 *
00123 *  -- LAPACK computational routine (version 3.4.0) --
00124 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00125 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00126 *     November 2011
00127 *
00128 *     .. Scalar Arguments ..
00129       CHARACTER          UPLO
00130       INTEGER            INFO, LDA, LDB, N, NRHS
00131 *     ..
00132 *     .. Array Arguments ..
00133       INTEGER            IPIV( * )
00134       COMPLEX*16         A( LDA, * ), B( LDB, * )
00135 *     ..
00136 *
00137 *  =====================================================================
00138 *
00139 *     .. Parameters ..
00140       COMPLEX*16         ONE
00141       PARAMETER          ( ONE = ( 1.0D+0, 0.0D+0 ) )
00142 *     ..
00143 *     .. Local Scalars ..
00144       LOGICAL            UPPER
00145       INTEGER            J, K, KP
00146       COMPLEX*16         AK, AKM1, AKM1K, BK, BKM1, DENOM
00147 *     ..
00148 *     .. External Functions ..
00149       LOGICAL            LSAME
00150       EXTERNAL           LSAME
00151 *     ..
00152 *     .. External Subroutines ..
00153       EXTERNAL           XERBLA, ZGEMV, ZGERU, ZSCAL, ZSWAP
00154 *     ..
00155 *     .. Intrinsic Functions ..
00156       INTRINSIC          MAX
00157 *     ..
00158 *     .. Executable Statements ..
00159 *
00160       INFO = 0
00161       UPPER = LSAME( UPLO, 'U' )
00162       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
00163          INFO = -1
00164       ELSE IF( N.LT.0 ) THEN
00165          INFO = -2
00166       ELSE IF( NRHS.LT.0 ) THEN
00167          INFO = -3
00168       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
00169          INFO = -5
00170       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
00171          INFO = -8
00172       END IF
00173       IF( INFO.NE.0 ) THEN
00174          CALL XERBLA( 'ZSYTRS', -INFO )
00175          RETURN
00176       END IF
00177 *
00178 *     Quick return if possible
00179 *
00180       IF( N.EQ.0 .OR. NRHS.EQ.0 )
00181      $   RETURN
00182 *
00183       IF( UPPER ) THEN
00184 *
00185 *        Solve A*X = B, where A = U*D*U**T.
00186 *
00187 *        First solve U*D*X = B, overwriting B with X.
00188 *
00189 *        K is the main loop index, decreasing from N to 1 in steps of
00190 *        1 or 2, depending on the size of the diagonal blocks.
00191 *
00192          K = N
00193    10    CONTINUE
00194 *
00195 *        If K < 1, exit from loop.
00196 *
00197          IF( K.LT.1 )
00198      $      GO TO 30
00199 *
00200          IF( IPIV( K ).GT.0 ) THEN
00201 *
00202 *           1 x 1 diagonal block
00203 *
00204 *           Interchange rows K and IPIV(K).
00205 *
00206             KP = IPIV( K )
00207             IF( KP.NE.K )
00208      $         CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
00209 *
00210 *           Multiply by inv(U(K)), where U(K) is the transformation
00211 *           stored in column K of A.
00212 *
00213             CALL ZGERU( K-1, NRHS, -ONE, A( 1, K ), 1, B( K, 1 ), LDB,
00214      $                  B( 1, 1 ), LDB )
00215 *
00216 *           Multiply by the inverse of the diagonal block.
00217 *
00218             CALL ZSCAL( NRHS, ONE / A( K, K ), B( K, 1 ), LDB )
00219             K = K - 1
00220          ELSE
00221 *
00222 *           2 x 2 diagonal block
00223 *
00224 *           Interchange rows K-1 and -IPIV(K).
00225 *
00226             KP = -IPIV( K )
00227             IF( KP.NE.K-1 )
00228      $         CALL ZSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ), LDB )
00229 *
00230 *           Multiply by inv(U(K)), where U(K) is the transformation
00231 *           stored in columns K-1 and K of A.
00232 *
00233             CALL ZGERU( K-2, NRHS, -ONE, A( 1, K ), 1, B( K, 1 ), LDB,
00234      $                  B( 1, 1 ), LDB )
00235             CALL ZGERU( K-2, NRHS, -ONE, A( 1, K-1 ), 1, B( K-1, 1 ),
00236      $                  LDB, B( 1, 1 ), LDB )
00237 *
00238 *           Multiply by the inverse of the diagonal block.
00239 *
00240             AKM1K = A( K-1, K )
00241             AKM1 = A( K-1, K-1 ) / AKM1K
00242             AK = A( K, K ) / AKM1K
00243             DENOM = AKM1*AK - ONE
00244             DO 20 J = 1, NRHS
00245                BKM1 = B( K-1, J ) / AKM1K
00246                BK = B( K, J ) / AKM1K
00247                B( K-1, J ) = ( AK*BKM1-BK ) / DENOM
00248                B( K, J ) = ( AKM1*BK-BKM1 ) / DENOM
00249    20       CONTINUE
00250             K = K - 2
00251          END IF
00252 *
00253          GO TO 10
00254    30    CONTINUE
00255 *
00256 *        Next solve U**T *X = B, overwriting B with X.
00257 *
00258 *        K is the main loop index, increasing from 1 to N in steps of
00259 *        1 or 2, depending on the size of the diagonal blocks.
00260 *
00261          K = 1
00262    40    CONTINUE
00263 *
00264 *        If K > N, exit from loop.
00265 *
00266          IF( K.GT.N )
00267      $      GO TO 50
00268 *
00269          IF( IPIV( K ).GT.0 ) THEN
00270 *
00271 *           1 x 1 diagonal block
00272 *
00273 *           Multiply by inv(U**T(K)), where U(K) is the transformation
00274 *           stored in column K of A.
00275 *
00276             CALL ZGEMV( 'Transpose', K-1, NRHS, -ONE, B, LDB, A( 1, K ),
00277      $                  1, ONE, B( K, 1 ), LDB )
00278 *
00279 *           Interchange rows K and IPIV(K).
00280 *
00281             KP = IPIV( K )
00282             IF( KP.NE.K )
00283      $         CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
00284             K = K + 1
00285          ELSE
00286 *
00287 *           2 x 2 diagonal block
00288 *
00289 *           Multiply by inv(U**T(K+1)), where U(K+1) is the transformation
00290 *           stored in columns K and K+1 of A.
00291 *
00292             CALL ZGEMV( 'Transpose', K-1, NRHS, -ONE, B, LDB, A( 1, K ),
00293      $                  1, ONE, B( K, 1 ), LDB )
00294             CALL ZGEMV( 'Transpose', K-1, NRHS, -ONE, B, LDB,
00295      $                  A( 1, K+1 ), 1, ONE, B( K+1, 1 ), LDB )
00296 *
00297 *           Interchange rows K and -IPIV(K).
00298 *
00299             KP = -IPIV( K )
00300             IF( KP.NE.K )
00301      $         CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
00302             K = K + 2
00303          END IF
00304 *
00305          GO TO 40
00306    50    CONTINUE
00307 *
00308       ELSE
00309 *
00310 *        Solve A*X = B, where A = L*D*L**T.
00311 *
00312 *        First solve L*D*X = B, overwriting B with X.
00313 *
00314 *        K is the main loop index, increasing from 1 to N in steps of
00315 *        1 or 2, depending on the size of the diagonal blocks.
00316 *
00317          K = 1
00318    60    CONTINUE
00319 *
00320 *        If K > N, exit from loop.
00321 *
00322          IF( K.GT.N )
00323      $      GO TO 80
00324 *
00325          IF( IPIV( K ).GT.0 ) THEN
00326 *
00327 *           1 x 1 diagonal block
00328 *
00329 *           Interchange rows K and IPIV(K).
00330 *
00331             KP = IPIV( K )
00332             IF( KP.NE.K )
00333      $         CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
00334 *
00335 *           Multiply by inv(L(K)), where L(K) is the transformation
00336 *           stored in column K of A.
00337 *
00338             IF( K.LT.N )
00339      $         CALL ZGERU( N-K, NRHS, -ONE, A( K+1, K ), 1, B( K, 1 ),
00340      $                     LDB, B( K+1, 1 ), LDB )
00341 *
00342 *           Multiply by the inverse of the diagonal block.
00343 *
00344             CALL ZSCAL( NRHS, ONE / A( K, K ), B( K, 1 ), LDB )
00345             K = K + 1
00346          ELSE
00347 *
00348 *           2 x 2 diagonal block
00349 *
00350 *           Interchange rows K+1 and -IPIV(K).
00351 *
00352             KP = -IPIV( K )
00353             IF( KP.NE.K+1 )
00354      $         CALL ZSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ), LDB )
00355 *
00356 *           Multiply by inv(L(K)), where L(K) is the transformation
00357 *           stored in columns K and K+1 of A.
00358 *
00359             IF( K.LT.N-1 ) THEN
00360                CALL ZGERU( N-K-1, NRHS, -ONE, A( K+2, K ), 1, B( K, 1 ),
00361      $                     LDB, B( K+2, 1 ), LDB )
00362                CALL ZGERU( N-K-1, NRHS, -ONE, A( K+2, K+1 ), 1,
00363      $                     B( K+1, 1 ), LDB, B( K+2, 1 ), LDB )
00364             END IF
00365 *
00366 *           Multiply by the inverse of the diagonal block.
00367 *
00368             AKM1K = A( K+1, K )
00369             AKM1 = A( K, K ) / AKM1K
00370             AK = A( K+1, K+1 ) / AKM1K
00371             DENOM = AKM1*AK - ONE
00372             DO 70 J = 1, NRHS
00373                BKM1 = B( K, J ) / AKM1K
00374                BK = B( K+1, J ) / AKM1K
00375                B( K, J ) = ( AK*BKM1-BK ) / DENOM
00376                B( K+1, J ) = ( AKM1*BK-BKM1 ) / DENOM
00377    70       CONTINUE
00378             K = K + 2
00379          END IF
00380 *
00381          GO TO 60
00382    80    CONTINUE
00383 *
00384 *        Next solve L**T *X = B, overwriting B with X.
00385 *
00386 *        K is the main loop index, decreasing from N to 1 in steps of
00387 *        1 or 2, depending on the size of the diagonal blocks.
00388 *
00389          K = N
00390    90    CONTINUE
00391 *
00392 *        If K < 1, exit from loop.
00393 *
00394          IF( K.LT.1 )
00395      $      GO TO 100
00396 *
00397          IF( IPIV( K ).GT.0 ) THEN
00398 *
00399 *           1 x 1 diagonal block
00400 *
00401 *           Multiply by inv(L**T(K)), where L(K) is the transformation
00402 *           stored in column K of A.
00403 *
00404             IF( K.LT.N )
00405      $         CALL ZGEMV( 'Transpose', N-K, NRHS, -ONE, B( K+1, 1 ),
00406      $                     LDB, A( K+1, K ), 1, ONE, B( K, 1 ), LDB )
00407 *
00408 *           Interchange rows K and IPIV(K).
00409 *
00410             KP = IPIV( K )
00411             IF( KP.NE.K )
00412      $         CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
00413             K = K - 1
00414          ELSE
00415 *
00416 *           2 x 2 diagonal block
00417 *
00418 *           Multiply by inv(L**T(K-1)), where L(K-1) is the transformation
00419 *           stored in columns K-1 and K of A.
00420 *
00421             IF( K.LT.N ) THEN
00422                CALL ZGEMV( 'Transpose', N-K, NRHS, -ONE, B( K+1, 1 ),
00423      $                     LDB, A( K+1, K ), 1, ONE, B( K, 1 ), LDB )
00424                CALL ZGEMV( 'Transpose', N-K, NRHS, -ONE, B( K+1, 1 ),
00425      $                     LDB, A( K+1, K-1 ), 1, ONE, B( K-1, 1 ),
00426      $                     LDB )
00427             END IF
00428 *
00429 *           Interchange rows K and -IPIV(K).
00430 *
00431             KP = -IPIV( K )
00432             IF( KP.NE.K )
00433      $         CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
00434             K = K - 2
00435          END IF
00436 *
00437          GO TO 90
00438   100    CONTINUE
00439       END IF
00440 *
00441       RETURN
00442 *
00443 *     End of ZSYTRS
00444 *
00445       END
 All Files Functions