Annotation of rpl/lapack/lapack/zhetrs_aa.f, revision 1.6

1.1       bertrand    1: *> \brief \b ZHETRS_AA
                      2: *
                      3: *  =========== DOCUMENTATION ===========
                      4: *
                      5: * Online html documentation available at
                      6: *            http://www.netlib.org/lapack/explore-html/
                      7: *
                      8: *> \htmlonly
                      9: *> Download ZHETRS_AA + dependencies
                     10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhetrs_aa.f">
                     11: *> [TGZ]</a>
                     12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhetrs_aa.f">
                     13: *> [ZIP]</a>
                     14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhetrs_aa.f">
                     15: *> [TXT]</a>
                     16: *> \endhtmlonly
                     17: *
                     18: *  Definition:
                     19: *  ===========
                     20: *
                     21: *       SUBROUTINE ZHETRS_AA( UPLO, N, NRHS, A, LDA, IPIV, B, LDB,
                     22: *                             WORK, LWORK, INFO )
                     23: *
                     24: *       .. Scalar Arguments ..
                     25: *       CHARACTER          UPLO
                     26: *       INTEGER            N, NRHS, LDA, LDB, LWORK, INFO
                     27: *       ..
                     28: *       .. Array Arguments ..
                     29: *       INTEGER            IPIV( * )
                     30: *       COMPLEX*16         A( LDA, * ), B( LDB, * ), WORK( * )
                     31: *       ..
                     32: *
                     33: *
                     34: *
                     35: *> \par Purpose:
                     36: *  =============
                     37: *>
                     38: *> \verbatim
                     39: *>
                     40: *> ZHETRS_AA solves a system of linear equations A*X = B with a complex
1.5       bertrand   41: *> hermitian matrix A using the factorization A = U**H*T*U or
                     42: *> A = L*T*L**H computed by ZHETRF_AA.
1.1       bertrand   43: *> \endverbatim
                     44: *
                     45: *  Arguments:
                     46: *  ==========
                     47: *
                     48: *> \param[in] UPLO
                     49: *> \verbatim
                     50: *>          UPLO is CHARACTER*1
                     51: *>          Specifies whether the details of the factorization are stored
                     52: *>          as an upper or lower triangular matrix.
1.5       bertrand   53: *>          = 'U':  Upper triangular, form is A = U**H*T*U;
1.1       bertrand   54: *>          = 'L':  Lower triangular, form is A = L*T*L**H.
                     55: *> \endverbatim
                     56: *>
                     57: *> \param[in] N
                     58: *> \verbatim
                     59: *>          N is INTEGER
                     60: *>          The order of the matrix A.  N >= 0.
                     61: *> \endverbatim
                     62: *>
                     63: *> \param[in] NRHS
                     64: *> \verbatim
                     65: *>          NRHS is INTEGER
                     66: *>          The number of right hand sides, i.e., the number of columns
                     67: *>          of the matrix B.  NRHS >= 0.
                     68: *> \endverbatim
                     69: *>
1.3       bertrand   70: *> \param[in] A
1.1       bertrand   71: *> \verbatim
                     72: *>          A is COMPLEX*16 array, dimension (LDA,N)
                     73: *>          Details of factors computed by ZHETRF_AA.
                     74: *> \endverbatim
                     75: *>
                     76: *> \param[in] LDA
                     77: *> \verbatim
                     78: *>          LDA is INTEGER
                     79: *>          The leading dimension of the array A.  LDA >= max(1,N).
                     80: *> \endverbatim
                     81: *>
                     82: *> \param[in] IPIV
                     83: *> \verbatim
                     84: *>          IPIV is INTEGER array, dimension (N)
                     85: *>          Details of the interchanges as computed by ZHETRF_AA.
                     86: *> \endverbatim
                     87: *>
                     88: *> \param[in,out] B
                     89: *> \verbatim
                     90: *>          B is COMPLEX*16 array, dimension (LDB,NRHS)
                     91: *>          On entry, the right hand side matrix B.
                     92: *>          On exit, the solution matrix X.
                     93: *> \endverbatim
                     94: *>
                     95: *> \param[in] LDB
                     96: *> \verbatim
                     97: *>          LDB is INTEGER
                     98: *>          The leading dimension of the array B.  LDB >= max(1,N).
                     99: *> \endverbatim
                    100: *>
1.5       bertrand  101: *> \param[out] WORK
1.1       bertrand  102: *> \verbatim
1.5       bertrand  103: *>          WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
1.1       bertrand  104: *> \endverbatim
                    105: *>
                    106: *> \param[in] LWORK
                    107: *> \verbatim
1.5       bertrand  108: *>          LWORK is INTEGER
                    109: *>          The dimension of the array WORK. LWORK >= max(1,3*N-2).
                    110: *> \endverbatim
1.1       bertrand  111: *>
                    112: *> \param[out] INFO
                    113: *> \verbatim
                    114: *>          INFO is INTEGER
                    115: *>          = 0:  successful exit
                    116: *>          < 0:  if INFO = -i, the i-th argument had an illegal value
                    117: *> \endverbatim
                    118: *
                    119: *  Authors:
                    120: *  ========
                    121: *
                    122: *> \author Univ. of Tennessee
                    123: *> \author Univ. of California Berkeley
                    124: *> \author Univ. of Colorado Denver
                    125: *> \author NAG Ltd.
                    126: *
                    127: *> \ingroup complex16HEcomputational
                    128: *
                    129: *  =====================================================================
                    130:       SUBROUTINE ZHETRS_AA( UPLO, N, NRHS, A, LDA, IPIV, B, LDB,
                    131:      $                      WORK, LWORK, INFO )
                    132: *
1.6     ! bertrand  133: *  -- LAPACK computational routine --
1.1       bertrand  134: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
                    135: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
                    136: *
                    137:       IMPLICIT NONE
                    138: *
                    139: *     .. Scalar Arguments ..
                    140:       CHARACTER          UPLO
                    141:       INTEGER            N, NRHS, LDA, LDB, LWORK, INFO
                    142: *     ..
                    143: *     .. Array Arguments ..
                    144:       INTEGER            IPIV( * )
                    145:       COMPLEX*16         A( LDA, * ), B( LDB, * ), WORK( * )
                    146: *     ..
                    147: *
                    148: *  =====================================================================
                    149: *
                    150:       COMPLEX*16         ONE
                    151:       PARAMETER          ( ONE = 1.0D+0 )
                    152: *     ..
                    153: *     .. Local Scalars ..
                    154:       LOGICAL            LQUERY, UPPER
                    155:       INTEGER            K, KP, LWKOPT
                    156: *     ..
                    157: *     .. External Functions ..
                    158:       LOGICAL            LSAME
                    159:       EXTERNAL           LSAME
                    160: *     ..
                    161: *     .. External Subroutines ..
1.3       bertrand  162:       EXTERNAL           ZGTSV, ZSWAP, ZTRSM, ZLACGV, ZLACPY, XERBLA
1.1       bertrand  163: *     ..
                    164: *     .. Intrinsic Functions ..
                    165:       INTRINSIC          MAX
                    166: *     ..
                    167: *     .. Executable Statements ..
                    168: *
                    169:       INFO = 0
                    170:       UPPER = LSAME( UPLO, 'U' )
                    171:       LQUERY = ( LWORK.EQ.-1 )
                    172:       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
                    173:          INFO = -1
                    174:       ELSE IF( N.LT.0 ) THEN
                    175:          INFO = -2
                    176:       ELSE IF( NRHS.LT.0 ) THEN
                    177:          INFO = -3
                    178:       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
                    179:          INFO = -5
                    180:       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
                    181:          INFO = -8
                    182:       ELSE IF( LWORK.LT.MAX( 1, 3*N-2 ) .AND. .NOT.LQUERY ) THEN
                    183:          INFO = -10
                    184:       END IF
                    185:       IF( INFO.NE.0 ) THEN
                    186:          CALL XERBLA( 'ZHETRS_AA', -INFO )
                    187:          RETURN
                    188:       ELSE IF( LQUERY ) THEN
                    189:          LWKOPT = (3*N-2)
                    190:          WORK( 1 ) = LWKOPT
                    191:          RETURN
                    192:       END IF
                    193: *
                    194: *     Quick return if possible
                    195: *
                    196:       IF( N.EQ.0 .OR. NRHS.EQ.0 )
                    197:      $   RETURN
                    198: *
                    199:       IF( UPPER ) THEN
                    200: *
1.5       bertrand  201: *        Solve A*X = B, where A = U**H*T*U.
                    202: *
                    203: *        1) Forward substitution with U**H
                    204: *
                    205:          IF( N.GT.1 ) THEN
                    206: *
                    207: *           Pivot, P**T * B -> B
1.1       bertrand  208: *
1.5       bertrand  209:             DO K = 1, N
                    210:                KP = IPIV( K )
                    211:                IF( KP.NE.K )
                    212:      $            CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
                    213:             END DO
1.1       bertrand  214: *
1.5       bertrand  215: *           Compute U**H \ B -> B    [ (U**H \P**T * B) ]
1.1       bertrand  216: *
1.5       bertrand  217:             CALL ZTRSM( 'L', 'U', 'C', 'U', N-1, NRHS, ONE, A( 1, 2 ),
                    218:      $                  LDA, B( 2, 1 ), LDB )
                    219:          END IF
1.1       bertrand  220: *
1.5       bertrand  221: *        2) Solve with triangular matrix T
1.1       bertrand  222: *
1.5       bertrand  223: *        Compute T \ B -> B   [ T \ (U**H \P**T * B) ]
1.1       bertrand  224: *
1.5       bertrand  225:          CALL ZLACPY( 'F', 1, N, A(1, 1), LDA+1, WORK(N), 1 )
1.1       bertrand  226:          IF( N.GT.1 ) THEN
                    227:              CALL ZLACPY( 'F', 1, N-1, A( 1, 2 ), LDA+1, WORK( 2*N ), 1)
1.5       bertrand  228:              CALL ZLACPY( 'F', 1, N-1, A( 1, 2 ), LDA+1, WORK( 1 ), 1 )
1.1       bertrand  229:              CALL ZLACGV( N-1, WORK( 1 ), 1 )
                    230:          END IF
1.5       bertrand  231:          CALL ZGTSV( N, NRHS, WORK(1), WORK(N), WORK(2*N), B, LDB,
                    232:      $               INFO )
                    233: *
                    234: *        3) Backward substitution with U
                    235: *
                    236:          IF( N.GT.1 ) THEN
1.1       bertrand  237: *
1.5       bertrand  238: *           Compute U \ B -> B   [ U \ (T \ (U**H \P**T * B) ) ]
1.1       bertrand  239: *
1.5       bertrand  240:             CALL ZTRSM( 'L', 'U', 'N', 'U', N-1, NRHS, ONE, A( 1, 2 ),
                    241:      $                  LDA, B(2, 1), LDB)
1.1       bertrand  242: *
1.5       bertrand  243: *           Pivot, P * B  [ P * (U**H \ (T \ (U \P**T * B) )) ]
1.1       bertrand  244: *
1.5       bertrand  245:             DO K = N, 1, -1
                    246:                KP = IPIV( K )
                    247:                IF( KP.NE.K )
                    248:      $            CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
                    249:             END DO
                    250:          END IF
1.1       bertrand  251: *
                    252:       ELSE
                    253: *
1.5       bertrand  254: *        Solve A*X = B, where A = L*T*L**H.
                    255: *
                    256: *        1) Forward substitution with L
1.1       bertrand  257: *
1.5       bertrand  258:          IF( N.GT.1 ) THEN
                    259: *
                    260: *           Pivot, P**T * B -> B
1.1       bertrand  261: *
1.5       bertrand  262:             DO K = 1, N
                    263:                KP = IPIV( K )
                    264:                IF( KP.NE.K )
                    265:      $            CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
                    266:             END DO
1.1       bertrand  267: *
1.5       bertrand  268: *           Compute L \ B -> B    [ (L \P**T * B) ]
1.1       bertrand  269: *
1.5       bertrand  270:             CALL ZTRSM( 'L', 'L', 'N', 'U', N-1, NRHS, ONE, A( 2, 1 ),
                    271:      $                  LDA, B(2, 1), LDB)
                    272:          END IF
                    273: *
                    274: *        2) Solve with triangular matrix T
1.1       bertrand  275: *
                    276: *        Compute T \ B -> B   [ T \ (L \P**T * B) ]
                    277: *
                    278:          CALL ZLACPY( 'F', 1, N, A(1, 1), LDA+1, WORK(N), 1)
                    279:          IF( N.GT.1 ) THEN
                    280:              CALL ZLACPY( 'F', 1, N-1, A( 2, 1 ), LDA+1, WORK( 1 ), 1)
                    281:              CALL ZLACPY( 'F', 1, N-1, A( 2, 1 ), LDA+1, WORK( 2*N ), 1)
                    282:              CALL ZLACGV( N-1, WORK( 2*N ), 1 )
                    283:          END IF
                    284:          CALL ZGTSV(N, NRHS, WORK(1), WORK(N), WORK(2*N), B, LDB,
                    285:      $              INFO)
                    286: *
1.5       bertrand  287: *        3) Backward substitution with L**H
                    288: *
                    289:          IF( N.GT.1 ) THEN
1.1       bertrand  290: *
1.5       bertrand  291: *           Compute L**H \ B -> B   [ L**H \ (T \ (L \P**T * B) ) ]
1.1       bertrand  292: *
1.5       bertrand  293:             CALL ZTRSM( 'L', 'L', 'C', 'U', N-1, NRHS, ONE, A( 2, 1 ),
                    294:      $                  LDA, B( 2, 1 ), LDB)
1.1       bertrand  295: *
1.5       bertrand  296: *           Pivot, P * B  [ P * (L**H \ (T \ (L \P**T * B) )) ]
                    297: *
                    298:             DO K = N, 1, -1
                    299:                KP = IPIV( K )
                    300:                IF( KP.NE.K )
                    301:      $            CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
                    302:             END DO
                    303:          END IF
1.1       bertrand  304: *
                    305:       END IF
                    306: *
                    307:       RETURN
                    308: *
                    309: *     End of ZHETRS_AA
                    310: *
                    311:       END

CVSweb interface <joel.bertrand@systella.fr>