Annotation of rpl/lapack/lapack/zhptrs.f, revision 1.5

1.1       bertrand    1:       SUBROUTINE ZHPTRS( UPLO, N, NRHS, AP, IPIV, B, LDB, INFO )
                      2: *
                      3: *  -- LAPACK routine (version 3.2) --
                      4: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
                      5: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
                      6: *     November 2006
                      7: *
                      8: *     .. Scalar Arguments ..
                      9:       CHARACTER          UPLO
                     10:       INTEGER            INFO, LDB, N, NRHS
                     11: *     ..
                     12: *     .. Array Arguments ..
                     13:       INTEGER            IPIV( * )
                     14:       COMPLEX*16         AP( * ), B( LDB, * )
                     15: *     ..
                     16: *
                     17: *  Purpose
                     18: *  =======
                     19: *
                     20: *  ZHPTRS solves a system of linear equations A*X = B with a complex
                     21: *  Hermitian matrix A stored in packed format using the factorization
                     22: *  A = U*D*U**H or A = L*D*L**H computed by ZHPTRF.
                     23: *
                     24: *  Arguments
                     25: *  =========
                     26: *
                     27: *  UPLO    (input) CHARACTER*1
                     28: *          Specifies whether the details of the factorization are stored
                     29: *          as an upper or lower triangular matrix.
                     30: *          = 'U':  Upper triangular, form is A = U*D*U**H;
                     31: *          = 'L':  Lower triangular, form is A = L*D*L**H.
                     32: *
                     33: *  N       (input) INTEGER
                     34: *          The order of the matrix A.  N >= 0.
                     35: *
                     36: *  NRHS    (input) INTEGER
                     37: *          The number of right hand sides, i.e., the number of columns
                     38: *          of the matrix B.  NRHS >= 0.
                     39: *
                     40: *  AP      (input) COMPLEX*16 array, dimension (N*(N+1)/2)
                     41: *          The block diagonal matrix D and the multipliers used to
                     42: *          obtain the factor U or L as computed by ZHPTRF, stored as a
                     43: *          packed triangular matrix.
                     44: *
                     45: *  IPIV    (input) INTEGER array, dimension (N)
                     46: *          Details of the interchanges and the block structure of D
                     47: *          as determined by ZHPTRF.
                     48: *
                     49: *  B       (input/output) COMPLEX*16 array, dimension (LDB,NRHS)
                     50: *          On entry, the right hand side matrix B.
                     51: *          On exit, the solution matrix X.
                     52: *
                     53: *  LDB     (input) INTEGER
                     54: *          The leading dimension of the array B.  LDB >= max(1,N).
                     55: *
                     56: *  INFO    (output) INTEGER
                     57: *          = 0:  successful exit
                     58: *          < 0: if INFO = -i, the i-th argument had an illegal value
                     59: *
                     60: *  =====================================================================
                     61: *
                     62: *     .. Parameters ..
                     63:       COMPLEX*16         ONE
                     64:       PARAMETER          ( ONE = ( 1.0D+0, 0.0D+0 ) )
                     65: *     ..
                     66: *     .. Local Scalars ..
                     67:       LOGICAL            UPPER
                     68:       INTEGER            J, K, KC, KP
                     69:       DOUBLE PRECISION   S
                     70:       COMPLEX*16         AK, AKM1, AKM1K, BK, BKM1, DENOM
                     71: *     ..
                     72: *     .. External Functions ..
                     73:       LOGICAL            LSAME
                     74:       EXTERNAL           LSAME
                     75: *     ..
                     76: *     .. External Subroutines ..
                     77:       EXTERNAL           XERBLA, ZDSCAL, ZGEMV, ZGERU, ZLACGV, ZSWAP
                     78: *     ..
                     79: *     .. Intrinsic Functions ..
                     80:       INTRINSIC          DBLE, DCONJG, MAX
                     81: *     ..
                     82: *     .. Executable Statements ..
                     83: *
                     84:       INFO = 0
                     85:       UPPER = LSAME( UPLO, 'U' )
                     86:       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
                     87:          INFO = -1
                     88:       ELSE IF( N.LT.0 ) THEN
                     89:          INFO = -2
                     90:       ELSE IF( NRHS.LT.0 ) THEN
                     91:          INFO = -3
                     92:       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
                     93:          INFO = -7
                     94:       END IF
                     95:       IF( INFO.NE.0 ) THEN
                     96:          CALL XERBLA( 'ZHPTRS', -INFO )
                     97:          RETURN
                     98:       END IF
                     99: *
                    100: *     Quick return if possible
                    101: *
                    102:       IF( N.EQ.0 .OR. NRHS.EQ.0 )
                    103:      $   RETURN
                    104: *
                    105:       IF( UPPER ) THEN
                    106: *
                    107: *        Solve A*X = B, where A = U*D*U'.
                    108: *
                    109: *        First solve U*D*X = B, overwriting B with X.
                    110: *
                    111: *        K is the main loop index, decreasing from N to 1 in steps of
                    112: *        1 or 2, depending on the size of the diagonal blocks.
                    113: *
                    114:          K = N
                    115:          KC = N*( N+1 ) / 2 + 1
                    116:    10    CONTINUE
                    117: *
                    118: *        If K < 1, exit from loop.
                    119: *
                    120:          IF( K.LT.1 )
                    121:      $      GO TO 30
                    122: *
                    123:          KC = KC - K
                    124:          IF( IPIV( K ).GT.0 ) THEN
                    125: *
                    126: *           1 x 1 diagonal block
                    127: *
                    128: *           Interchange rows K and IPIV(K).
                    129: *
                    130:             KP = IPIV( K )
                    131:             IF( KP.NE.K )
                    132:      $         CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
                    133: *
                    134: *           Multiply by inv(U(K)), where U(K) is the transformation
                    135: *           stored in column K of A.
                    136: *
                    137:             CALL ZGERU( K-1, NRHS, -ONE, AP( KC ), 1, B( K, 1 ), LDB,
                    138:      $                  B( 1, 1 ), LDB )
                    139: *
                    140: *           Multiply by the inverse of the diagonal block.
                    141: *
                    142:             S = DBLE( ONE ) / DBLE( AP( KC+K-1 ) )
                    143:             CALL ZDSCAL( NRHS, S, B( K, 1 ), LDB )
                    144:             K = K - 1
                    145:          ELSE
                    146: *
                    147: *           2 x 2 diagonal block
                    148: *
                    149: *           Interchange rows K-1 and -IPIV(K).
                    150: *
                    151:             KP = -IPIV( K )
                    152:             IF( KP.NE.K-1 )
                    153:      $         CALL ZSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ), LDB )
                    154: *
                    155: *           Multiply by inv(U(K)), where U(K) is the transformation
                    156: *           stored in columns K-1 and K of A.
                    157: *
                    158:             CALL ZGERU( K-2, NRHS, -ONE, AP( KC ), 1, B( K, 1 ), LDB,
                    159:      $                  B( 1, 1 ), LDB )
                    160:             CALL ZGERU( K-2, NRHS, -ONE, AP( KC-( K-1 ) ), 1,
                    161:      $                  B( K-1, 1 ), LDB, B( 1, 1 ), LDB )
                    162: *
                    163: *           Multiply by the inverse of the diagonal block.
                    164: *
                    165:             AKM1K = AP( KC+K-2 )
                    166:             AKM1 = AP( KC-1 ) / AKM1K
                    167:             AK = AP( KC+K-1 ) / DCONJG( AKM1K )
                    168:             DENOM = AKM1*AK - ONE
                    169:             DO 20 J = 1, NRHS
                    170:                BKM1 = B( K-1, J ) / AKM1K
                    171:                BK = B( K, J ) / DCONJG( AKM1K )
                    172:                B( K-1, J ) = ( AK*BKM1-BK ) / DENOM
                    173:                B( K, J ) = ( AKM1*BK-BKM1 ) / DENOM
                    174:    20       CONTINUE
                    175:             KC = KC - K + 1
                    176:             K = K - 2
                    177:          END IF
                    178: *
                    179:          GO TO 10
                    180:    30    CONTINUE
                    181: *
                    182: *        Next solve U'*X = B, overwriting B with X.
                    183: *
                    184: *        K is the main loop index, increasing from 1 to N in steps of
                    185: *        1 or 2, depending on the size of the diagonal blocks.
                    186: *
                    187:          K = 1
                    188:          KC = 1
                    189:    40    CONTINUE
                    190: *
                    191: *        If K > N, exit from loop.
                    192: *
                    193:          IF( K.GT.N )
                    194:      $      GO TO 50
                    195: *
                    196:          IF( IPIV( K ).GT.0 ) THEN
                    197: *
                    198: *           1 x 1 diagonal block
                    199: *
                    200: *           Multiply by inv(U'(K)), where U(K) is the transformation
                    201: *           stored in column K of A.
                    202: *
                    203:             IF( K.GT.1 ) THEN
                    204:                CALL ZLACGV( NRHS, B( K, 1 ), LDB )
                    205:                CALL ZGEMV( 'Conjugate transpose', K-1, NRHS, -ONE, B,
                    206:      $                     LDB, AP( KC ), 1, ONE, B( K, 1 ), LDB )
                    207:                CALL ZLACGV( NRHS, B( K, 1 ), LDB )
                    208:             END IF
                    209: *
                    210: *           Interchange rows K and IPIV(K).
                    211: *
                    212:             KP = IPIV( K )
                    213:             IF( KP.NE.K )
                    214:      $         CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
                    215:             KC = KC + K
                    216:             K = K + 1
                    217:          ELSE
                    218: *
                    219: *           2 x 2 diagonal block
                    220: *
                    221: *           Multiply by inv(U'(K+1)), where U(K+1) is the transformation
                    222: *           stored in columns K and K+1 of A.
                    223: *
                    224:             IF( K.GT.1 ) THEN
                    225:                CALL ZLACGV( NRHS, B( K, 1 ), LDB )
                    226:                CALL ZGEMV( 'Conjugate transpose', K-1, NRHS, -ONE, B,
                    227:      $                     LDB, AP( KC ), 1, ONE, B( K, 1 ), LDB )
                    228:                CALL ZLACGV( NRHS, B( K, 1 ), LDB )
                    229: *
                    230:                CALL ZLACGV( NRHS, B( K+1, 1 ), LDB )
                    231:                CALL ZGEMV( 'Conjugate transpose', K-1, NRHS, -ONE, B,
                    232:      $                     LDB, AP( KC+K ), 1, ONE, B( K+1, 1 ), LDB )
                    233:                CALL ZLACGV( NRHS, B( K+1, 1 ), LDB )
                    234:             END IF
                    235: *
                    236: *           Interchange rows K and -IPIV(K).
                    237: *
                    238:             KP = -IPIV( K )
                    239:             IF( KP.NE.K )
                    240:      $         CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
                    241:             KC = KC + 2*K + 1
                    242:             K = K + 2
                    243:          END IF
                    244: *
                    245:          GO TO 40
                    246:    50    CONTINUE
                    247: *
                    248:       ELSE
                    249: *
                    250: *        Solve A*X = B, where A = L*D*L'.
                    251: *
                    252: *        First solve L*D*X = B, overwriting B with X.
                    253: *
                    254: *        K is the main loop index, increasing from 1 to N in steps of
                    255: *        1 or 2, depending on the size of the diagonal blocks.
                    256: *
                    257:          K = 1
                    258:          KC = 1
                    259:    60    CONTINUE
                    260: *
                    261: *        If K > N, exit from loop.
                    262: *
                    263:          IF( K.GT.N )
                    264:      $      GO TO 80
                    265: *
                    266:          IF( IPIV( K ).GT.0 ) THEN
                    267: *
                    268: *           1 x 1 diagonal block
                    269: *
                    270: *           Interchange rows K and IPIV(K).
                    271: *
                    272:             KP = IPIV( K )
                    273:             IF( KP.NE.K )
                    274:      $         CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
                    275: *
                    276: *           Multiply by inv(L(K)), where L(K) is the transformation
                    277: *           stored in column K of A.
                    278: *
                    279:             IF( K.LT.N )
                    280:      $         CALL ZGERU( N-K, NRHS, -ONE, AP( KC+1 ), 1, B( K, 1 ),
                    281:      $                     LDB, B( K+1, 1 ), LDB )
                    282: *
                    283: *           Multiply by the inverse of the diagonal block.
                    284: *
                    285:             S = DBLE( ONE ) / DBLE( AP( KC ) )
                    286:             CALL ZDSCAL( NRHS, S, B( K, 1 ), LDB )
                    287:             KC = KC + N - K + 1
                    288:             K = K + 1
                    289:          ELSE
                    290: *
                    291: *           2 x 2 diagonal block
                    292: *
                    293: *           Interchange rows K+1 and -IPIV(K).
                    294: *
                    295:             KP = -IPIV( K )
                    296:             IF( KP.NE.K+1 )
                    297:      $         CALL ZSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ), LDB )
                    298: *
                    299: *           Multiply by inv(L(K)), where L(K) is the transformation
                    300: *           stored in columns K and K+1 of A.
                    301: *
                    302:             IF( K.LT.N-1 ) THEN
                    303:                CALL ZGERU( N-K-1, NRHS, -ONE, AP( KC+2 ), 1, B( K, 1 ),
                    304:      $                     LDB, B( K+2, 1 ), LDB )
                    305:                CALL ZGERU( N-K-1, NRHS, -ONE, AP( KC+N-K+2 ), 1,
                    306:      $                     B( K+1, 1 ), LDB, B( K+2, 1 ), LDB )
                    307:             END IF
                    308: *
                    309: *           Multiply by the inverse of the diagonal block.
                    310: *
                    311:             AKM1K = AP( KC+1 )
                    312:             AKM1 = AP( KC ) / DCONJG( AKM1K )
                    313:             AK = AP( KC+N-K+1 ) / AKM1K
                    314:             DENOM = AKM1*AK - ONE
                    315:             DO 70 J = 1, NRHS
                    316:                BKM1 = B( K, J ) / DCONJG( AKM1K )
                    317:                BK = B( K+1, J ) / AKM1K
                    318:                B( K, J ) = ( AK*BKM1-BK ) / DENOM
                    319:                B( K+1, J ) = ( AKM1*BK-BKM1 ) / DENOM
                    320:    70       CONTINUE
                    321:             KC = KC + 2*( N-K ) + 1
                    322:             K = K + 2
                    323:          END IF
                    324: *
                    325:          GO TO 60
                    326:    80    CONTINUE
                    327: *
                    328: *        Next solve L'*X = B, overwriting B with X.
                    329: *
                    330: *        K is the main loop index, decreasing from N to 1 in steps of
                    331: *        1 or 2, depending on the size of the diagonal blocks.
                    332: *
                    333:          K = N
                    334:          KC = N*( N+1 ) / 2 + 1
                    335:    90    CONTINUE
                    336: *
                    337: *        If K < 1, exit from loop.
                    338: *
                    339:          IF( K.LT.1 )
                    340:      $      GO TO 100
                    341: *
                    342:          KC = KC - ( N-K+1 )
                    343:          IF( IPIV( K ).GT.0 ) THEN
                    344: *
                    345: *           1 x 1 diagonal block
                    346: *
                    347: *           Multiply by inv(L'(K)), where L(K) is the transformation
                    348: *           stored in column K of A.
                    349: *
                    350:             IF( K.LT.N ) THEN
                    351:                CALL ZLACGV( NRHS, B( K, 1 ), LDB )
                    352:                CALL ZGEMV( 'Conjugate transpose', N-K, NRHS, -ONE,
                    353:      $                     B( K+1, 1 ), LDB, AP( KC+1 ), 1, ONE,
                    354:      $                     B( K, 1 ), LDB )
                    355:                CALL ZLACGV( NRHS, B( K, 1 ), LDB )
                    356:             END IF
                    357: *
                    358: *           Interchange rows K and IPIV(K).
                    359: *
                    360:             KP = IPIV( K )
                    361:             IF( KP.NE.K )
                    362:      $         CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
                    363:             K = K - 1
                    364:          ELSE
                    365: *
                    366: *           2 x 2 diagonal block
                    367: *
                    368: *           Multiply by inv(L'(K-1)), where L(K-1) is the transformation
                    369: *           stored in columns K-1 and K of A.
                    370: *
                    371:             IF( K.LT.N ) THEN
                    372:                CALL ZLACGV( NRHS, B( K, 1 ), LDB )
                    373:                CALL ZGEMV( 'Conjugate transpose', N-K, NRHS, -ONE,
                    374:      $                     B( K+1, 1 ), LDB, AP( KC+1 ), 1, ONE,
                    375:      $                     B( K, 1 ), LDB )
                    376:                CALL ZLACGV( NRHS, B( K, 1 ), LDB )
                    377: *
                    378:                CALL ZLACGV( NRHS, B( K-1, 1 ), LDB )
                    379:                CALL ZGEMV( 'Conjugate transpose', N-K, NRHS, -ONE,
                    380:      $                     B( K+1, 1 ), LDB, AP( KC-( N-K ) ), 1, ONE,
                    381:      $                     B( K-1, 1 ), LDB )
                    382:                CALL ZLACGV( NRHS, B( K-1, 1 ), LDB )
                    383:             END IF
                    384: *
                    385: *           Interchange rows K and -IPIV(K).
                    386: *
                    387:             KP = -IPIV( K )
                    388:             IF( KP.NE.K )
                    389:      $         CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
                    390:             KC = KC - ( N-K+2 )
                    391:             K = K - 2
                    392:          END IF
                    393: *
                    394:          GO TO 90
                    395:   100    CONTINUE
                    396:       END IF
                    397: *
                    398:       RETURN
                    399: *
                    400: *     End of ZHPTRS
                    401: *
                    402:       END

CVSweb interface <joel.bertrand@systella.fr>