File:  [local] / rpl / lapack / blas / ztrsv.f
Revision 1.7: download - view: text, annotated - select for diffs - revision graph
Fri Jul 22 07:38:03 2011 UTC (12 years, 10 months ago) by bertrand
Branches: MAIN
CVS tags: rpl-4_1_3, rpl-4_1_2, rpl-4_1_1, HEAD
En route vers la 4.4.1.

    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: *
   15: *     A*x = b,   or   A**T*x = b,   or   A**H*x = b,
   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: *
   42: *              TRANS = 'T' or 't'   A**T*x = b.
   43: *
   44: *              TRANS = 'C' or 'c'   A**H*x = b.
   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: *
  227: *        Form  x := inv( A**T )*x  or  x := inv( A**H )*x.
  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>