Return to ztbsv.f CVS log

Up to [local] / rpl / lapack / blas

Cohérence.

    1: *> \brief \b ZTBSV
    2: *
    3: *  =========== DOCUMENTATION ===========
    4: *
    5: * Online html documentation available at 
    6: *            http://www.netlib.org/lapack/explore-html/ 
    7: *
    8: *  Definition:
    9: *  ===========
   10: *
   11: *       SUBROUTINE ZTBSV(UPLO,TRANS,DIAG,N,K,A,LDA,X,INCX)
   12: * 
   13: *       .. Scalar Arguments ..
   14: *       INTEGER INCX,K,LDA,N
   15: *       CHARACTER DIAG,TRANS,UPLO
   16: *       ..
   17: *       .. Array Arguments ..
   18: *       COMPLEX*16 A(LDA,*),X(*)
   19: *       ..
   20: *  
   21: *
   22: *> \par Purpose:
   23: *  =============
   24: *>
   25: *> \verbatim
   26: *>
   27: *> ZTBSV  solves one of the systems of equations
   28: *>
   29: *>    A*x = b,   or   A**T*x = b,   or   A**H*x = b,
   30: *>
   31: *> where b and x are n element vectors and A is an n by n unit, or
   32: *> non-unit, upper or lower triangular band matrix, with ( k + 1 )
   33: *> diagonals.
   34: *>
   35: *> No test for singularity or near-singularity is included in this
   36: *> routine. Such tests must be performed before calling this routine.
   37: *> \endverbatim
   38: *
   39: *  Arguments:
   40: *  ==========
   41: *
   42: *> \param[in] UPLO
   43: *> \verbatim
   44: *>          UPLO is CHARACTER*1
   45: *>           On entry, UPLO specifies whether the matrix is an upper or
   46: *>           lower triangular matrix as follows:
   47: *>
   48: *>              UPLO = 'U' or 'u'   A is an upper triangular matrix.
   49: *>
   50: *>              UPLO = 'L' or 'l'   A is a lower triangular matrix.
   51: *> \endverbatim
   52: *>
   53: *> \param[in] TRANS
   54: *> \verbatim
   55: *>          TRANS is CHARACTER*1
   56: *>           On entry, TRANS specifies the equations to be solved as
   57: *>           follows:
   58: *>
   59: *>              TRANS = 'N' or 'n'   A*x = b.
   60: *>
   61: *>              TRANS = 'T' or 't'   A**T*x = b.
   62: *>
   63: *>              TRANS = 'C' or 'c'   A**H*x = b.
   64: *> \endverbatim
   65: *>
   66: *> \param[in] DIAG
   67: *> \verbatim
   68: *>          DIAG is CHARACTER*1
   69: *>           On entry, DIAG specifies whether or not A is unit
   70: *>           triangular as follows:
   71: *>
   72: *>              DIAG = 'U' or 'u'   A is assumed to be unit triangular.
   73: *>
   74: *>              DIAG = 'N' or 'n'   A is not assumed to be unit
   75: *>                                  triangular.
   76: *> \endverbatim
   77: *>
   78: *> \param[in] N
   79: *> \verbatim
   80: *>          N is INTEGER
   81: *>           On entry, N specifies the order of the matrix A.
   82: *>           N must be at least zero.
   83: *> \endverbatim
   84: *>
   85: *> \param[in] K
   86: *> \verbatim
   87: *>          K is INTEGER
   88: *>           On entry with UPLO = 'U' or 'u', K specifies the number of
   89: *>           super-diagonals of the matrix A.
   90: *>           On entry with UPLO = 'L' or 'l', K specifies the number of
   91: *>           sub-diagonals of the matrix A.
   92: *>           K must satisfy  0 .le. K.
   93: *> \endverbatim
   94: *>
   95: *> \param[in] A
   96: *> \verbatim
   97: *>          A is COMPLEX*16 array of DIMENSION ( LDA, n ).
   98: *>           Before entry with UPLO = 'U' or 'u', the leading ( k + 1 )
   99: *>           by n part of the array A must contain the upper triangular
  100: *>           band part of the matrix of coefficients, supplied column by
  101: *>           column, with the leading diagonal of the matrix in row
  102: *>           ( k + 1 ) of the array, the first super-diagonal starting at
  103: *>           position 2 in row k, and so on. The top left k by k triangle
  104: *>           of the array A is not referenced.
  105: *>           The following program segment will transfer an upper
  106: *>           triangular band matrix from conventional full matrix storage
  107: *>           to band storage:
  108: *>
  109: *>                 DO 20, J = 1, N
  110: *>                    M = K + 1 - J
  111: *>                    DO 10, I = MAX( 1, J - K ), J
  112: *>                       A( M + I, J ) = matrix( I, J )
  113: *>              10    CONTINUE
  114: *>              20 CONTINUE
  115: *>
  116: *>           Before entry with UPLO = 'L' or 'l', the leading ( k + 1 )
  117: *>           by n part of the array A must contain the lower triangular
  118: *>           band part of the matrix of coefficients, supplied column by
  119: *>           column, with the leading diagonal of the matrix in row 1 of
  120: *>           the array, the first sub-diagonal starting at position 1 in
  121: *>           row 2, and so on. The bottom right k by k triangle of the
  122: *>           array A is not referenced.
  123: *>           The following program segment will transfer a lower
  124: *>           triangular band matrix from conventional full matrix storage
  125: *>           to band storage:
  126: *>
  127: *>                 DO 20, J = 1, N
  128: *>                    M = 1 - J
  129: *>                    DO 10, I = J, MIN( N, J + K )
  130: *>                       A( M + I, J ) = matrix( I, J )
  131: *>              10    CONTINUE
  132: *>              20 CONTINUE
  133: *>
  134: *>           Note that when DIAG = 'U' or 'u' the elements of the array A
  135: *>           corresponding to the diagonal elements of the matrix are not
  136: *>           referenced, but are assumed to be unity.
  137: *> \endverbatim
  138: *>
  139: *> \param[in] LDA
  140: *> \verbatim
  141: *>          LDA is INTEGER
  142: *>           On entry, LDA specifies the first dimension of A as declared
  143: *>           in the calling (sub) program. LDA must be at least
  144: *>           ( k + 1 ).
  145: *> \endverbatim
  146: *>
  147: *> \param[in,out] X
  148: *> \verbatim
  149: *>          X is COMPLEX*16 array of dimension at least
  150: *>           ( 1 + ( n - 1 )*abs( INCX ) ).
  151: *>           Before entry, the incremented array X must contain the n
  152: *>           element right-hand side vector b. On exit, X is overwritten
  153: *>           with the solution vector x.
  154: *> \endverbatim
  155: *>
  156: *> \param[in] INCX
  157: *> \verbatim
  158: *>          INCX is INTEGER
  159: *>           On entry, INCX specifies the increment for the elements of
  160: *>           X. INCX must not be zero.
  161: *> \endverbatim
  162: *
  163: *  Authors:
  164: *  ========
  165: *
  166: *> \author Univ. of Tennessee 
  167: *> \author Univ. of California Berkeley 
  168: *> \author Univ. of Colorado Denver 
  169: *> \author NAG Ltd. 
  170: *
  171: *> \date November 2011
  172: *
  173: *> \ingroup complex16_blas_level2
  174: *
  175: *> \par Further Details:
  176: *  =====================
  177: *>
  178: *> \verbatim
  179: *>
  180: *>  Level 2 Blas routine.
  181: *>
  182: *>  -- Written on 22-October-1986.
  183: *>     Jack Dongarra, Argonne National Lab.
  184: *>     Jeremy Du Croz, Nag Central Office.
  185: *>     Sven Hammarling, Nag Central Office.
  186: *>     Richard Hanson, Sandia National Labs.
  187: *> \endverbatim
  188: *>
  189: *  =====================================================================
  190:       SUBROUTINE ZTBSV(UPLO,TRANS,DIAG,N,K,A,LDA,X,INCX)
  191: *
  192: *  -- Reference BLAS level2 routine (version 3.4.0) --
  193: *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
  194: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  195: *     November 2011
  196: *
  197: *     .. Scalar Arguments ..
  198:       INTEGER INCX,K,LDA,N
  199:       CHARACTER DIAG,TRANS,UPLO
  200: *     ..
  201: *     .. Array Arguments ..
  202:       COMPLEX*16 A(LDA,*),X(*)
  203: *     ..
  204: *
  205: *  =====================================================================
  206: *
  207: *     .. Parameters ..
  208:       COMPLEX*16 ZERO
  209:       PARAMETER (ZERO= (0.0D+0,0.0D+0))
  210: *     ..
  211: *     .. Local Scalars ..
  212:       COMPLEX*16 TEMP
  213:       INTEGER I,INFO,IX,J,JX,KPLUS1,KX,L
  214:       LOGICAL NOCONJ,NOUNIT
  215: *     ..
  216: *     .. External Functions ..
  217:       LOGICAL LSAME
  218:       EXTERNAL LSAME
  219: *     ..
  220: *     .. External Subroutines ..
  221:       EXTERNAL XERBLA
  222: *     ..
  223: *     .. Intrinsic Functions ..
  224:       INTRINSIC DCONJG,MAX,MIN
  225: *     ..
  226: *
  227: *     Test the input parameters.
  228: *
  229:       INFO = 0
  230:       IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN
  231:           INFO = 1
  232:       ELSE IF (.NOT.LSAME(TRANS,'N') .AND. .NOT.LSAME(TRANS,'T') .AND.
  233:      +         .NOT.LSAME(TRANS,'C')) THEN
  234:           INFO = 2
  235:       ELSE IF (.NOT.LSAME(DIAG,'U') .AND. .NOT.LSAME(DIAG,'N')) THEN
  236:           INFO = 3
  237:       ELSE IF (N.LT.0) THEN
  238:           INFO = 4
  239:       ELSE IF (K.LT.0) THEN
  240:           INFO = 5
  241:       ELSE IF (LDA.LT. (K+1)) THEN
  242:           INFO = 7
  243:       ELSE IF (INCX.EQ.0) THEN
  244:           INFO = 9
  245:       END IF
  246:       IF (INFO.NE.0) THEN
  247:           CALL XERBLA('ZTBSV ',INFO)
  248:           RETURN
  249:       END IF
  250: *
  251: *     Quick return if possible.
  252: *
  253:       IF (N.EQ.0) RETURN
  254: *
  255:       NOCONJ = LSAME(TRANS,'T')
  256:       NOUNIT = LSAME(DIAG,'N')
  257: *
  258: *     Set up the start point in X if the increment is not unity. This
  259: *     will be  ( N - 1 )*INCX  too small for descending loops.
  260: *
  261:       IF (INCX.LE.0) THEN
  262:           KX = 1 - (N-1)*INCX
  263:       ELSE IF (INCX.NE.1) THEN
  264:           KX = 1
  265:       END IF
  266: *
  267: *     Start the operations. In this version the elements of A are
  268: *     accessed by sequentially with one pass through A.
  269: *
  270:       IF (LSAME(TRANS,'N')) THEN
  271: *
  272: *        Form  x := inv( A )*x.
  273: *
  274:           IF (LSAME(UPLO,'U')) THEN
  275:               KPLUS1 = K + 1
  276:               IF (INCX.EQ.1) THEN
  277:                   DO 20 J = N,1,-1
  278:                       IF (X(J).NE.ZERO) THEN
  279:                           L = KPLUS1 - J
  280:                           IF (NOUNIT) X(J) = X(J)/A(KPLUS1,J)
  281:                           TEMP = X(J)
  282:                           DO 10 I = J - 1,MAX(1,J-K),-1
  283:                               X(I) = X(I) - TEMP*A(L+I,J)
  284:    10                     CONTINUE
  285:                       END IF
  286:    20             CONTINUE
  287:               ELSE
  288:                   KX = KX + (N-1)*INCX
  289:                   JX = KX
  290:                   DO 40 J = N,1,-1
  291:                       KX = KX - INCX
  292:                       IF (X(JX).NE.ZERO) THEN
  293:                           IX = KX
  294:                           L = KPLUS1 - J
  295:                           IF (NOUNIT) X(JX) = X(JX)/A(KPLUS1,J)
  296:                           TEMP = X(JX)
  297:                           DO 30 I = J - 1,MAX(1,J-K),-1
  298:                               X(IX) = X(IX) - TEMP*A(L+I,J)
  299:                               IX = IX - INCX
  300:    30                     CONTINUE
  301:                       END IF
  302:                       JX = JX - INCX
  303:    40             CONTINUE
  304:               END IF
  305:           ELSE
  306:               IF (INCX.EQ.1) THEN
  307:                   DO 60 J = 1,N
  308:                       IF (X(J).NE.ZERO) THEN
  309:                           L = 1 - J
  310:                           IF (NOUNIT) X(J) = X(J)/A(1,J)
  311:                           TEMP = X(J)
  312:                           DO 50 I = J + 1,MIN(N,J+K)
  313:                               X(I) = X(I) - TEMP*A(L+I,J)
  314:    50                     CONTINUE
  315:                       END IF
  316:    60             CONTINUE
  317:               ELSE
  318:                   JX = KX
  319:                   DO 80 J = 1,N
  320:                       KX = KX + INCX
  321:                       IF (X(JX).NE.ZERO) THEN
  322:                           IX = KX
  323:                           L = 1 - J
  324:                           IF (NOUNIT) X(JX) = X(JX)/A(1,J)
  325:                           TEMP = X(JX)
  326:                           DO 70 I = J + 1,MIN(N,J+K)
  327:                               X(IX) = X(IX) - TEMP*A(L+I,J)
  328:                               IX = IX + INCX
  329:    70                     CONTINUE
  330:                       END IF
  331:                       JX = JX + INCX
  332:    80             CONTINUE
  333:               END IF
  334:           END IF
  335:       ELSE
  336: *
  337: *        Form  x := inv( A**T )*x  or  x := inv( A**H )*x.
  338: *
  339:           IF (LSAME(UPLO,'U')) THEN
  340:               KPLUS1 = K + 1
  341:               IF (INCX.EQ.1) THEN
  342:                   DO 110 J = 1,N
  343:                       TEMP = X(J)
  344:                       L = KPLUS1 - J
  345:                       IF (NOCONJ) THEN
  346:                           DO 90 I = MAX(1,J-K),J - 1
  347:                               TEMP = TEMP - A(L+I,J)*X(I)
  348:    90                     CONTINUE
  349:                           IF (NOUNIT) TEMP = TEMP/A(KPLUS1,J)
  350:                       ELSE
  351:                           DO 100 I = MAX(1,J-K),J - 1
  352:                               TEMP = TEMP - DCONJG(A(L+I,J))*X(I)
  353:   100                     CONTINUE
  354:                           IF (NOUNIT) TEMP = TEMP/DCONJG(A(KPLUS1,J))
  355:                       END IF
  356:                       X(J) = TEMP
  357:   110             CONTINUE
  358:               ELSE
  359:                   JX = KX
  360:                   DO 140 J = 1,N
  361:                       TEMP = X(JX)
  362:                       IX = KX
  363:                       L = KPLUS1 - J
  364:                       IF (NOCONJ) THEN
  365:                           DO 120 I = MAX(1,J-K),J - 1
  366:                               TEMP = TEMP - A(L+I,J)*X(IX)
  367:                               IX = IX + INCX
  368:   120                     CONTINUE
  369:                           IF (NOUNIT) TEMP = TEMP/A(KPLUS1,J)
  370:                       ELSE
  371:                           DO 130 I = MAX(1,J-K),J - 1
  372:                               TEMP = TEMP - DCONJG(A(L+I,J))*X(IX)
  373:                               IX = IX + INCX
  374:   130                     CONTINUE
  375:                           IF (NOUNIT) TEMP = TEMP/DCONJG(A(KPLUS1,J))
  376:                       END IF
  377:                       X(JX) = TEMP
  378:                       JX = JX + INCX
  379:                       IF (J.GT.K) KX = KX + INCX
  380:   140             CONTINUE
  381:               END IF
  382:           ELSE
  383:               IF (INCX.EQ.1) THEN
  384:                   DO 170 J = N,1,-1
  385:                       TEMP = X(J)
  386:                       L = 1 - J
  387:                       IF (NOCONJ) THEN
  388:                           DO 150 I = MIN(N,J+K),J + 1,-1
  389:                               TEMP = TEMP - A(L+I,J)*X(I)
  390:   150                     CONTINUE
  391:                           IF (NOUNIT) TEMP = TEMP/A(1,J)
  392:                       ELSE
  393:                           DO 160 I = MIN(N,J+K),J + 1,-1
  394:                               TEMP = TEMP - DCONJG(A(L+I,J))*X(I)
  395:   160                     CONTINUE
  396:                           IF (NOUNIT) TEMP = TEMP/DCONJG(A(1,J))
  397:                       END IF
  398:                       X(J) = TEMP
  399:   170             CONTINUE
  400:               ELSE
  401:                   KX = KX + (N-1)*INCX
  402:                   JX = KX
  403:                   DO 200 J = N,1,-1
  404:                       TEMP = X(JX)
  405:                       IX = KX
  406:                       L = 1 - J
  407:                       IF (NOCONJ) THEN
  408:                           DO 180 I = MIN(N,J+K),J + 1,-1
  409:                               TEMP = TEMP - A(L+I,J)*X(IX)
  410:                               IX = IX - INCX
  411:   180                     CONTINUE
  412:                           IF (NOUNIT) TEMP = TEMP/A(1,J)
  413:                       ELSE
  414:                           DO 190 I = MIN(N,J+K),J + 1,-1
  415:                               TEMP = TEMP - DCONJG(A(L+I,J))*X(IX)
  416:                               IX = IX - INCX
  417:   190                     CONTINUE
  418:                           IF (NOUNIT) TEMP = TEMP/DCONJG(A(1,J))
  419:                       END IF
  420:                       X(JX) = TEMP
  421:                       JX = JX - INCX
  422:                       IF ((N-J).GE.K) KX = KX - INCX
  423:   200             CONTINUE
  424:               END IF
  425:           END IF
  426:       END IF
  427: *
  428:       RETURN
  429: *
  430: *     End of ZTBSV .
  431: *
  432:       END