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