Return to dsytrs_3.f CVS log

Up to [local] / rpl / lapack / lapack

Mise à jour de Lapack.

    1: *> \brief \b DSYTRS_3
    2: *
    3: *  =========== DOCUMENTATION ===========
    4: *
    5: * Online html documentation available at
    6: *            http://www.netlib.org/lapack/explore-html/
    7: *
    8: *> \htmlonly
    9: *> Download DSYTRS_3 + dependencies
   10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsytrs_3.f">
   11: *> [TGZ]</a>
   12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsytrs_3.f">
   13: *> [ZIP]</a>
   14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsytrs_3.f">
   15: *> [TXT]</a>
   16: *> \endhtmlonly
   17: *
   18: *  Definition:
   19: *  ===========
   20: *
   21: *       SUBROUTINE DSYTRS_3( UPLO, N, NRHS, A, LDA, E, IPIV, B, LDB,
   22: *                            INFO )
   23: *
   24: *       .. Scalar Arguments ..
   25: *       CHARACTER          UPLO
   26: *       INTEGER            INFO, LDA, LDB, N, NRHS
   27: *       ..
   28: *       .. Array Arguments ..
   29: *       INTEGER            IPIV( * )
   30: *       DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), E( * )
   31: *       ..
   32: *
   33: *
   34: *> \par Purpose:
   35: *  =============
   36: *>
   37: *> \verbatim
   38: *> DSYTRS_3 solves a system of linear equations A * X = B with a real
   39: *> symmetric matrix A using the factorization computed
   40: *> by DSYTRF_RK or DSYTRF_BK:
   41: *>
   42: *>    A = P*U*D*(U**T)*(P**T) or A = P*L*D*(L**T)*(P**T),
   43: *>
   44: *> where U (or L) is unit upper (or lower) triangular matrix,
   45: *> U**T (or L**T) is the transpose of U (or L), P is a permutation
   46: *> matrix, P**T is the transpose of P, and D is symmetric and block
   47: *> diagonal with 1-by-1 and 2-by-2 diagonal blocks.
   48: *>
   49: *> This algorithm is using Level 3 BLAS.
   50: *> \endverbatim
   51: *
   52: *  Arguments:
   53: *  ==========
   54: *
   55: *> \param[in] UPLO
   56: *> \verbatim
   57: *>          UPLO is CHARACTER*1
   58: *>          Specifies whether the details of the factorization are
   59: *>          stored as an upper or lower triangular matrix:
   60: *>          = 'U':  Upper triangular, form is A = P*U*D*(U**T)*(P**T);
   61: *>          = 'L':  Lower triangular, form is A = P*L*D*(L**T)*(P**T).
   62: *> \endverbatim
   63: *>
   64: *> \param[in] N
   65: *> \verbatim
   66: *>          N is INTEGER
   67: *>          The order of the matrix A.  N >= 0.
   68: *> \endverbatim
   69: *>
   70: *> \param[in] NRHS
   71: *> \verbatim
   72: *>          NRHS is INTEGER
   73: *>          The number of right hand sides, i.e., the number of columns
   74: *>          of the matrix B.  NRHS >= 0.
   75: *> \endverbatim
   76: *>
   77: *> \param[in] A
   78: *> \verbatim
   79: *>          A is DOUBLE PRECISION array, dimension (LDA,N)
   80: *>          Diagonal of the block diagonal matrix D and factors U or L
   81: *>          as computed by DSYTRF_RK and DSYTRF_BK:
   82: *>            a) ONLY diagonal elements of the symmetric block diagonal
   83: *>               matrix D on the diagonal of A, i.e. D(k,k) = A(k,k);
   84: *>               (superdiagonal (or subdiagonal) elements of D
   85: *>                should be provided on entry in array E), and
   86: *>            b) If UPLO = 'U': factor U in the superdiagonal part of A.
   87: *>               If UPLO = 'L': factor L in the subdiagonal part of A.
   88: *> \endverbatim
   89: *>
   90: *> \param[in] LDA
   91: *> \verbatim
   92: *>          LDA is INTEGER
   93: *>          The leading dimension of the array A.  LDA >= max(1,N).
   94: *> \endverbatim
   95: *>
   96: *> \param[in] E
   97: *> \verbatim
   98: *>          E is DOUBLE PRECISION array, dimension (N)
   99: *>          On entry, contains the superdiagonal (or subdiagonal)
  100: *>          elements of the symmetric block diagonal matrix D
  101: *>          with 1-by-1 or 2-by-2 diagonal blocks, where
  102: *>          If UPLO = 'U': E(i) = D(i-1,i),i=2:N, E(1) not referenced;
  103: *>          If UPLO = 'L': E(i) = D(i+1,i),i=1:N-1, E(N) not referenced.
  104: *>
  105: *>          NOTE: For 1-by-1 diagonal block D(k), where
  106: *>          1 <= k <= N, the element E(k) is not referenced in both
  107: *>          UPLO = 'U' or UPLO = 'L' cases.
  108: *> \endverbatim
  109: *>
  110: *> \param[in] IPIV
  111: *> \verbatim
  112: *>          IPIV is INTEGER array, dimension (N)
  113: *>          Details of the interchanges and the block structure of D
  114: *>          as determined by DSYTRF_RK or DSYTRF_BK.
  115: *> \endverbatim
  116: *>
  117: *> \param[in,out] B
  118: *> \verbatim
  119: *>          B is DOUBLE PRECISION array, dimension (LDB,NRHS)
  120: *>          On entry, the right hand side matrix B.
  121: *>          On exit, the solution matrix X.
  122: *> \endverbatim
  123: *>
  124: *> \param[in] LDB
  125: *> \verbatim
  126: *>          LDB is INTEGER
  127: *>          The leading dimension of the array B.  LDB >= max(1,N).
  128: *> \endverbatim
  129: *>
  130: *> \param[out] INFO
  131: *> \verbatim
  132: *>          INFO is INTEGER
  133: *>          = 0:  successful exit
  134: *>          < 0:  if INFO = -i, the i-th argument had an illegal value
  135: *> \endverbatim
  136: *
  137: *  Authors:
  138: *  ========
  139: *
  140: *> \author Univ. of Tennessee
  141: *> \author Univ. of California Berkeley
  142: *> \author Univ. of Colorado Denver
  143: *> \author NAG Ltd.
  144: *
  145: *> \date June 2017
  146: *
  147: *> \ingroup doubleSYcomputational
  148: *
  149: *> \par Contributors:
  150: *  ==================
  151: *>
  152: *> \verbatim
  153: *>
  154: *>  June 2017,  Igor Kozachenko,
  155: *>                  Computer Science Division,
  156: *>                  University of California, Berkeley
  157: *>
  158: *>  September 2007, Sven Hammarling, Nicholas J. Higham, Craig Lucas,
  159: *>                  School of Mathematics,
  160: *>                  University of Manchester
  161: *>
  162: *> \endverbatim
  163: *
  164: *  =====================================================================
  165:       SUBROUTINE DSYTRS_3( UPLO, N, NRHS, A, LDA, E, IPIV, B, LDB,
  166:      $                     INFO )
  167: *
  168: *  -- LAPACK computational routine (version 3.7.1) --
  169: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  170: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  171: *     June 2017
  172: *
  173: *     .. Scalar Arguments ..
  174:       CHARACTER          UPLO
  175:       INTEGER            INFO, LDA, LDB, N, NRHS
  176: *     ..
  177: *     .. Array Arguments ..
  178:       INTEGER            IPIV( * )
  179:       DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), E( * )
  180: *     ..
  181: *
  182: *  =====================================================================
  183: *
  184: *     .. Parameters ..
  185:       DOUBLE PRECISION   ONE
  186:       PARAMETER          ( ONE = 1.0D+0 )
  187: *     ..
  188: *     .. Local Scalars ..
  189:       LOGICAL            UPPER
  190:       INTEGER            I, J, K, KP
  191:       DOUBLE PRECISION   AK, AKM1, AKM1K, BK, BKM1, DENOM
  192: *     ..
  193: *     .. External Functions ..
  194:       LOGICAL            LSAME
  195:       EXTERNAL           LSAME
  196: *     ..
  197: *     .. External Subroutines ..
  198:       EXTERNAL           DSCAL, DSWAP, DTRSM, XERBLA
  199: *     ..
  200: *     .. Intrinsic Functions ..
  201:       INTRINSIC          ABS, MAX
  202: *     ..
  203: *     .. Executable Statements ..
  204: *
  205:       INFO = 0
  206:       UPPER = LSAME( UPLO, 'U' )
  207:       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
  208:          INFO = -1
  209:       ELSE IF( N.LT.0 ) THEN
  210:          INFO = -2
  211:       ELSE IF( NRHS.LT.0 ) THEN
  212:          INFO = -3
  213:       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
  214:          INFO = -5
  215:       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
  216:          INFO = -9
  217:       END IF
  218:       IF( INFO.NE.0 ) THEN
  219:          CALL XERBLA( 'DSYTRS_3', -INFO )
  220:          RETURN
  221:       END IF
  222: *
  223: *     Quick return if possible
  224: *
  225:       IF( N.EQ.0 .OR. NRHS.EQ.0 )
  226:      $   RETURN
  227: *
  228:       IF( UPPER ) THEN
  229: *
  230: *        Begin Upper
  231: *
  232: *        Solve A*X = B, where A = U*D*U**T.
  233: *
  234: *        P**T * B
  235: *
  236: *        Interchange rows K and IPIV(K) of matrix B in the same order
  237: *        that the formation order of IPIV(I) vector for Upper case.
  238: *
  239: *        (We can do the simple loop over IPIV with decrement -1,
  240: *        since the ABS value of IPIV( I ) represents the row index
  241: *        of the interchange with row i in both 1x1 and 2x2 pivot cases)
  242: *
  243:          DO K = N, 1, -1
  244:             KP = ABS( IPIV( K ) )
  245:             IF( KP.NE.K ) THEN
  246:                CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
  247:             END IF
  248:          END DO
  249: *
  250: *        Compute (U \P**T * B) -> B    [ (U \P**T * B) ]
  251: *
  252:          CALL DTRSM( 'L', 'U', 'N', 'U', N, NRHS, ONE, A, LDA, B, LDB )
  253: *
  254: *        Compute D \ B -> B   [ D \ (U \P**T * B) ]
  255: *
  256:          I = N
  257:          DO WHILE ( I.GE.1 )
  258:             IF( IPIV( I ).GT.0 ) THEN
  259:                CALL DSCAL( NRHS, ONE / A( I, I ), B( I, 1 ), LDB )
  260:             ELSE IF ( I.GT.1 ) THEN
  261:                AKM1K = E( I )
  262:                AKM1 = A( I-1, I-1 ) / AKM1K
  263:                AK = A( I, I ) / AKM1K
  264:                DENOM = AKM1*AK - ONE
  265:                DO J = 1, NRHS
  266:                   BKM1 = B( I-1, J ) / AKM1K
  267:                   BK = B( I, J ) / AKM1K
  268:                   B( I-1, J ) = ( AK*BKM1-BK ) / DENOM
  269:                   B( I, J ) = ( AKM1*BK-BKM1 ) / DENOM
  270:                END DO
  271:                I = I - 1
  272:             END IF
  273:             I = I - 1
  274:          END DO
  275: *
  276: *        Compute (U**T \ B) -> B   [ U**T \ (D \ (U \P**T * B) ) ]
  277: *
  278:          CALL DTRSM( 'L', 'U', 'T', 'U', N, NRHS, ONE, A, LDA, B, LDB )
  279: *
  280: *        P * B  [ P * (U**T \ (D \ (U \P**T * B) )) ]
  281: *
  282: *        Interchange rows K and IPIV(K) of matrix B in reverse order
  283: *        from the formation order of IPIV(I) vector for Upper case.
  284: *
  285: *        (We can do the simple loop over IPIV with increment 1,
  286: *        since the ABS value of IPIV(I) represents the row index
  287: *        of the interchange with row i in both 1x1 and 2x2 pivot cases)
  288: *
  289:          DO K = 1, N
  290:             KP = ABS( IPIV( K ) )
  291:             IF( KP.NE.K ) THEN
  292:                CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
  293:             END IF
  294:          END DO
  295: *
  296:       ELSE
  297: *
  298: *        Begin Lower
  299: *
  300: *        Solve A*X = B, where A = L*D*L**T.
  301: *
  302: *        P**T * B
  303: *        Interchange rows K and IPIV(K) of matrix B in the same order
  304: *        that the formation order of IPIV(I) vector for Lower case.
  305: *
  306: *        (We can do the simple loop over IPIV with increment 1,
  307: *        since the ABS value of IPIV(I) represents the row index
  308: *        of the interchange with row i in both 1x1 and 2x2 pivot cases)
  309: *
  310:          DO K = 1, N
  311:             KP = ABS( IPIV( K ) )
  312:             IF( KP.NE.K ) THEN
  313:                CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
  314:             END IF
  315:          END DO
  316: *
  317: *        Compute (L \P**T * B) -> B    [ (L \P**T * B) ]
  318: *
  319:          CALL DTRSM( 'L', 'L', 'N', 'U', N, NRHS, ONE, A, LDA, B, LDB )
  320: *
  321: *        Compute D \ B -> B   [ D \ (L \P**T * B) ]
  322: *
  323:          I = 1
  324:          DO WHILE ( I.LE.N )
  325:             IF( IPIV( I ).GT.0 ) THEN
  326:                CALL DSCAL( NRHS, ONE / A( I, I ), B( I, 1 ), LDB )
  327:             ELSE IF( I.LT.N ) THEN
  328:                AKM1K = E( I )
  329:                AKM1 = A( I, I ) / AKM1K
  330:                AK = A( I+1, I+1 ) / AKM1K
  331:                DENOM = AKM1*AK - ONE
  332:                DO  J = 1, NRHS
  333:                   BKM1 = B( I, J ) / AKM1K
  334:                   BK = B( I+1, J ) / AKM1K
  335:                   B( I, J ) = ( AK*BKM1-BK ) / DENOM
  336:                   B( I+1, J ) = ( AKM1*BK-BKM1 ) / DENOM
  337:                END DO
  338:                I = I + 1
  339:             END IF
  340:             I = I + 1
  341:          END DO
  342: *
  343: *        Compute (L**T \ B) -> B   [ L**T \ (D \ (L \P**T * B) ) ]
  344: *
  345:          CALL DTRSM('L', 'L', 'T', 'U', N, NRHS, ONE, A, LDA, B, LDB )
  346: *
  347: *        P * B  [ P * (L**T \ (D \ (L \P**T * B) )) ]
  348: *
  349: *        Interchange rows K and IPIV(K) of matrix B in reverse order
  350: *        from the formation order of IPIV(I) vector for Lower case.
  351: *
  352: *        (We can do the simple loop over IPIV with decrement -1,
  353: *        since the ABS value of IPIV(I) represents the row index
  354: *        of the interchange with row i in both 1x1 and 2x2 pivot cases)
  355: *
  356:          DO K = N, 1, -1
  357:             KP = ABS( IPIV( K ) )
  358:             IF( KP.NE.K ) THEN
  359:                CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
  360:             END IF
  361:          END DO
  362: *
  363: *        END Lower
  364: *
  365:       END IF
  366: *
  367:       RETURN
  368: *
  369: *     End of DSYTRS_3
  370: *
  371:       END