Annotation of rpl/lapack/blas/ztrsv.f, revision 1.7

1.1       bertrand    1:       SUBROUTINE ZTRSV(UPLO,TRANS,DIAG,N,A,LDA,X,INCX)
                      2: *     .. Scalar Arguments ..
                      3:       INTEGER INCX,LDA,N
                      4:       CHARACTER DIAG,TRANS,UPLO
                      5: *     ..
                      6: *     .. Array Arguments ..
                      7:       DOUBLE COMPLEX A(LDA,*),X(*)
                      8: *     ..
                      9: *
                     10: *  Purpose
                     11: *  =======
                     12: *
                     13: *  ZTRSV  solves one of the systems of equations
                     14: *
1.7     ! bertrand   15: *     A*x = b,   or   A**T*x = b,   or   A**H*x = b,
1.1       bertrand   16: *
                     17: *  where b and x are n element vectors and A is an n by n unit, or
                     18: *  non-unit, upper or lower triangular matrix.
                     19: *
                     20: *  No test for singularity or near-singularity is included in this
                     21: *  routine. Such tests must be performed before calling this routine.
                     22: *
                     23: *  Arguments
                     24: *  ==========
                     25: *
                     26: *  UPLO   - CHARACTER*1.
                     27: *           On entry, UPLO specifies whether the matrix is an upper or
                     28: *           lower triangular matrix as follows:
                     29: *
                     30: *              UPLO = 'U' or 'u'   A is an upper triangular matrix.
                     31: *
                     32: *              UPLO = 'L' or 'l'   A is a lower triangular matrix.
                     33: *
                     34: *           Unchanged on exit.
                     35: *
                     36: *  TRANS  - CHARACTER*1.
                     37: *           On entry, TRANS specifies the equations to be solved as
                     38: *           follows:
                     39: *
                     40: *              TRANS = 'N' or 'n'   A*x = b.
                     41: *
1.7     ! bertrand   42: *              TRANS = 'T' or 't'   A**T*x = b.
1.1       bertrand   43: *
1.7     ! bertrand   44: *              TRANS = 'C' or 'c'   A**H*x = b.
1.1       bertrand   45: *
                     46: *           Unchanged on exit.
                     47: *
                     48: *  DIAG   - CHARACTER*1.
                     49: *           On entry, DIAG specifies whether or not A is unit
                     50: *           triangular as follows:
                     51: *
                     52: *              DIAG = 'U' or 'u'   A is assumed to be unit triangular.
                     53: *
                     54: *              DIAG = 'N' or 'n'   A is not assumed to be unit
                     55: *                                  triangular.
                     56: *
                     57: *           Unchanged on exit.
                     58: *
                     59: *  N      - INTEGER.
                     60: *           On entry, N specifies the order of the matrix A.
                     61: *           N must be at least zero.
                     62: *           Unchanged on exit.
                     63: *
                     64: *  A      - COMPLEX*16       array of DIMENSION ( LDA, n ).
                     65: *           Before entry with  UPLO = 'U' or 'u', the leading n by n
                     66: *           upper triangular part of the array A must contain the upper
                     67: *           triangular matrix and the strictly lower triangular part of
                     68: *           A is not referenced.
                     69: *           Before entry with UPLO = 'L' or 'l', the leading n by n
                     70: *           lower triangular part of the array A must contain the lower
                     71: *           triangular matrix and the strictly upper triangular part of
                     72: *           A is not referenced.
                     73: *           Note that when  DIAG = 'U' or 'u', the diagonal elements of
                     74: *           A are not referenced either, but are assumed to be unity.
                     75: *           Unchanged on exit.
                     76: *
                     77: *  LDA    - INTEGER.
                     78: *           On entry, LDA specifies the first dimension of A as declared
                     79: *           in the calling (sub) program. LDA must be at least
                     80: *           max( 1, n ).
                     81: *           Unchanged on exit.
                     82: *
                     83: *  X      - COMPLEX*16       array of dimension at least
                     84: *           ( 1 + ( n - 1 )*abs( INCX ) ).
                     85: *           Before entry, the incremented array X must contain the n
                     86: *           element right-hand side vector b. On exit, X is overwritten
                     87: *           with the solution vector x.
                     88: *
                     89: *  INCX   - INTEGER.
                     90: *           On entry, INCX specifies the increment for the elements of
                     91: *           X. INCX must not be zero.
                     92: *           Unchanged on exit.
                     93: *
                     94: *  Further Details
                     95: *  ===============
                     96: *
                     97: *  Level 2 Blas routine.
                     98: *
                     99: *  -- Written on 22-October-1986.
                    100: *     Jack Dongarra, Argonne National Lab.
                    101: *     Jeremy Du Croz, Nag Central Office.
                    102: *     Sven Hammarling, Nag Central Office.
                    103: *     Richard Hanson, Sandia National Labs.
                    104: *
                    105: *  =====================================================================
                    106: *
                    107: *     .. Parameters ..
                    108:       DOUBLE COMPLEX ZERO
                    109:       PARAMETER (ZERO= (0.0D+0,0.0D+0))
                    110: *     ..
                    111: *     .. Local Scalars ..
                    112:       DOUBLE COMPLEX TEMP
                    113:       INTEGER I,INFO,IX,J,JX,KX
                    114:       LOGICAL NOCONJ,NOUNIT
                    115: *     ..
                    116: *     .. External Functions ..
                    117:       LOGICAL LSAME
                    118:       EXTERNAL LSAME
                    119: *     ..
                    120: *     .. External Subroutines ..
                    121:       EXTERNAL XERBLA
                    122: *     ..
                    123: *     .. Intrinsic Functions ..
                    124:       INTRINSIC DCONJG,MAX
                    125: *     ..
                    126: *
                    127: *     Test the input parameters.
                    128: *
                    129:       INFO = 0
                    130:       IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN
                    131:           INFO = 1
                    132:       ELSE IF (.NOT.LSAME(TRANS,'N') .AND. .NOT.LSAME(TRANS,'T') .AND.
                    133:      +         .NOT.LSAME(TRANS,'C')) THEN
                    134:           INFO = 2
                    135:       ELSE IF (.NOT.LSAME(DIAG,'U') .AND. .NOT.LSAME(DIAG,'N')) THEN
                    136:           INFO = 3
                    137:       ELSE IF (N.LT.0) THEN
                    138:           INFO = 4
                    139:       ELSE IF (LDA.LT.MAX(1,N)) THEN
                    140:           INFO = 6
                    141:       ELSE IF (INCX.EQ.0) THEN
                    142:           INFO = 8
                    143:       END IF
                    144:       IF (INFO.NE.0) THEN
                    145:           CALL XERBLA('ZTRSV ',INFO)
                    146:           RETURN
                    147:       END IF
                    148: *
                    149: *     Quick return if possible.
                    150: *
                    151:       IF (N.EQ.0) RETURN
                    152: *
                    153:       NOCONJ = LSAME(TRANS,'T')
                    154:       NOUNIT = LSAME(DIAG,'N')
                    155: *
                    156: *     Set up the start point in X if the increment is not unity. This
                    157: *     will be  ( N - 1 )*INCX  too small for descending loops.
                    158: *
                    159:       IF (INCX.LE.0) THEN
                    160:           KX = 1 - (N-1)*INCX
                    161:       ELSE IF (INCX.NE.1) THEN
                    162:           KX = 1
                    163:       END IF
                    164: *
                    165: *     Start the operations. In this version the elements of A are
                    166: *     accessed sequentially with one pass through A.
                    167: *
                    168:       IF (LSAME(TRANS,'N')) THEN
                    169: *
                    170: *        Form  x := inv( A )*x.
                    171: *
                    172:           IF (LSAME(UPLO,'U')) THEN
                    173:               IF (INCX.EQ.1) THEN
                    174:                   DO 20 J = N,1,-1
                    175:                       IF (X(J).NE.ZERO) THEN
                    176:                           IF (NOUNIT) X(J) = X(J)/A(J,J)
                    177:                           TEMP = X(J)
                    178:                           DO 10 I = J - 1,1,-1
                    179:                               X(I) = X(I) - TEMP*A(I,J)
                    180:    10                     CONTINUE
                    181:                       END IF
                    182:    20             CONTINUE
                    183:               ELSE
                    184:                   JX = KX + (N-1)*INCX
                    185:                   DO 40 J = N,1,-1
                    186:                       IF (X(JX).NE.ZERO) THEN
                    187:                           IF (NOUNIT) X(JX) = X(JX)/A(J,J)
                    188:                           TEMP = X(JX)
                    189:                           IX = JX
                    190:                           DO 30 I = J - 1,1,-1
                    191:                               IX = IX - INCX
                    192:                               X(IX) = X(IX) - TEMP*A(I,J)
                    193:    30                     CONTINUE
                    194:                       END IF
                    195:                       JX = JX - INCX
                    196:    40             CONTINUE
                    197:               END IF
                    198:           ELSE
                    199:               IF (INCX.EQ.1) THEN
                    200:                   DO 60 J = 1,N
                    201:                       IF (X(J).NE.ZERO) THEN
                    202:                           IF (NOUNIT) X(J) = X(J)/A(J,J)
                    203:                           TEMP = X(J)
                    204:                           DO 50 I = J + 1,N
                    205:                               X(I) = X(I) - TEMP*A(I,J)
                    206:    50                     CONTINUE
                    207:                       END IF
                    208:    60             CONTINUE
                    209:               ELSE
                    210:                   JX = KX
                    211:                   DO 80 J = 1,N
                    212:                       IF (X(JX).NE.ZERO) THEN
                    213:                           IF (NOUNIT) X(JX) = X(JX)/A(J,J)
                    214:                           TEMP = X(JX)
                    215:                           IX = JX
                    216:                           DO 70 I = J + 1,N
                    217:                               IX = IX + INCX
                    218:                               X(IX) = X(IX) - TEMP*A(I,J)
                    219:    70                     CONTINUE
                    220:                       END IF
                    221:                       JX = JX + INCX
                    222:    80             CONTINUE
                    223:               END IF
                    224:           END IF
                    225:       ELSE
                    226: *
1.7     ! bertrand  227: *        Form  x := inv( A**T )*x  or  x := inv( A**H )*x.
1.1       bertrand  228: *
                    229:           IF (LSAME(UPLO,'U')) THEN
                    230:               IF (INCX.EQ.1) THEN
                    231:                   DO 110 J = 1,N
                    232:                       TEMP = X(J)
                    233:                       IF (NOCONJ) THEN
                    234:                           DO 90 I = 1,J - 1
                    235:                               TEMP = TEMP - A(I,J)*X(I)
                    236:    90                     CONTINUE
                    237:                           IF (NOUNIT) TEMP = TEMP/A(J,J)
                    238:                       ELSE
                    239:                           DO 100 I = 1,J - 1
                    240:                               TEMP = TEMP - DCONJG(A(I,J))*X(I)
                    241:   100                     CONTINUE
                    242:                           IF (NOUNIT) TEMP = TEMP/DCONJG(A(J,J))
                    243:                       END IF
                    244:                       X(J) = TEMP
                    245:   110             CONTINUE
                    246:               ELSE
                    247:                   JX = KX
                    248:                   DO 140 J = 1,N
                    249:                       IX = KX
                    250:                       TEMP = X(JX)
                    251:                       IF (NOCONJ) THEN
                    252:                           DO 120 I = 1,J - 1
                    253:                               TEMP = TEMP - A(I,J)*X(IX)
                    254:                               IX = IX + INCX
                    255:   120                     CONTINUE
                    256:                           IF (NOUNIT) TEMP = TEMP/A(J,J)
                    257:                       ELSE
                    258:                           DO 130 I = 1,J - 1
                    259:                               TEMP = TEMP - DCONJG(A(I,J))*X(IX)
                    260:                               IX = IX + INCX
                    261:   130                     CONTINUE
                    262:                           IF (NOUNIT) TEMP = TEMP/DCONJG(A(J,J))
                    263:                       END IF
                    264:                       X(JX) = TEMP
                    265:                       JX = JX + INCX
                    266:   140             CONTINUE
                    267:               END IF
                    268:           ELSE
                    269:               IF (INCX.EQ.1) THEN
                    270:                   DO 170 J = N,1,-1
                    271:                       TEMP = X(J)
                    272:                       IF (NOCONJ) THEN
                    273:                           DO 150 I = N,J + 1,-1
                    274:                               TEMP = TEMP - A(I,J)*X(I)
                    275:   150                     CONTINUE
                    276:                           IF (NOUNIT) TEMP = TEMP/A(J,J)
                    277:                       ELSE
                    278:                           DO 160 I = N,J + 1,-1
                    279:                               TEMP = TEMP - DCONJG(A(I,J))*X(I)
                    280:   160                     CONTINUE
                    281:                           IF (NOUNIT) TEMP = TEMP/DCONJG(A(J,J))
                    282:                       END IF
                    283:                       X(J) = TEMP
                    284:   170             CONTINUE
                    285:               ELSE
                    286:                   KX = KX + (N-1)*INCX
                    287:                   JX = KX
                    288:                   DO 200 J = N,1,-1
                    289:                       IX = KX
                    290:                       TEMP = X(JX)
                    291:                       IF (NOCONJ) THEN
                    292:                           DO 180 I = N,J + 1,-1
                    293:                               TEMP = TEMP - A(I,J)*X(IX)
                    294:                               IX = IX - INCX
                    295:   180                     CONTINUE
                    296:                           IF (NOUNIT) TEMP = TEMP/A(J,J)
                    297:                       ELSE
                    298:                           DO 190 I = N,J + 1,-1
                    299:                               TEMP = TEMP - DCONJG(A(I,J))*X(IX)
                    300:                               IX = IX - INCX
                    301:   190                     CONTINUE
                    302:                           IF (NOUNIT) TEMP = TEMP/DCONJG(A(J,J))
                    303:                       END IF
                    304:                       X(JX) = TEMP
                    305:                       JX = JX - INCX
                    306:   200             CONTINUE
                    307:               END IF
                    308:           END IF
                    309:       END IF
                    310: *
                    311:       RETURN
                    312: *
                    313: *     End of ZTRSV .
                    314: *
                    315:       END

CVSweb interface <joel.bertrand@systella.fr>