File:  [local] / rpl / lapack / lapack / zhesv_rk.f
Revision 1.4: download - view: text, annotated - select for diffs - revision graph
Mon Aug 7 08:39:24 2023 UTC (8 months, 3 weeks ago) by bertrand
Branches: MAIN
CVS tags: rpl-4_1_35, rpl-4_1_34, HEAD
Première mise à jour de lapack et blas.

    1: *> \brief <b> ZHESV_RK computes the solution to system of linear equations A * X = B for SY matrices</b>
    2: *
    3: *  =========== DOCUMENTATION ===========
    4: *
    5: * Online html documentation available at
    6: *            http://www.netlib.org/lapack/explore-html/
    7: *
    8: *> \htmlonly
    9: *> Download ZHESV_RK + dependencies
   10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhesv_rk.f">
   11: *> [TGZ]</a>
   12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhesv_rk.f">
   13: *> [ZIP]</a>
   14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhesv_rk.f">
   15: *> [TXT]</a>
   16: *> \endhtmlonly
   17: *
   18: *  Definition:
   19: *  ===========
   20: *
   21: *       SUBROUTINE ZHESV_RK( UPLO, N, NRHS, A, LDA, E, IPIV, B, LDB,
   22: *                            WORK, LWORK, INFO )
   23: *
   24: *       .. Scalar Arguments ..
   25: *       CHARACTER          UPLO
   26: *       INTEGER            INFO, LDA, LDB, LWORK, N, NRHS
   27: *       ..
   28: *       .. Array Arguments ..
   29: *       INTEGER            IPIV( * )
   30: *       COMPLEX*16         A( LDA, * ), B( LDB, * ), E( * ), WORK( * )
   31: *       ..
   32: *
   33: *
   34: *> \par Purpose:
   35: *  =============
   36: *>
   37: *> \verbatim
   38: *> ZHESV_RK computes the solution to a complex system of linear
   39: *> equations A * X = B, where A is an N-by-N Hermitian matrix
   40: *> and X and B are N-by-NRHS matrices.
   41: *>
   42: *> The bounded Bunch-Kaufman (rook) diagonal pivoting method is used
   43: *> to factor A as
   44: *>    A = P*U*D*(U**H)*(P**T),  if UPLO = 'U', or
   45: *>    A = P*L*D*(L**H)*(P**T),  if UPLO = 'L',
   46: *> where U (or L) is unit upper (or lower) triangular matrix,
   47: *> U**H (or L**H) is the conjugate of U (or L), P is a permutation
   48: *> matrix, P**T is the transpose of P, and D is Hermitian and block
   49: *> diagonal with 1-by-1 and 2-by-2 diagonal blocks.
   50: *>
   51: *> ZHETRF_RK is called to compute the factorization of a complex
   52: *> Hermitian matrix.  The factored form of A is then used to solve
   53: *> the system of equations A * X = B by calling BLAS3 routine ZHETRS_3.
   54: *> \endverbatim
   55: *
   56: *  Arguments:
   57: *  ==========
   58: *
   59: *> \param[in] UPLO
   60: *> \verbatim
   61: *>          UPLO is CHARACTER*1
   62: *>          Specifies whether the upper or lower triangular part of the
   63: *>          Hermitian matrix A is stored:
   64: *>          = 'U':  Upper triangle of A is stored;
   65: *>          = 'L':  Lower triangle of A is stored.
   66: *> \endverbatim
   67: *>
   68: *> \param[in] N
   69: *> \verbatim
   70: *>          N is INTEGER
   71: *>          The number of linear equations, i.e., the order of the
   72: *>          matrix A.  N >= 0.
   73: *> \endverbatim
   74: *>
   75: *> \param[in] NRHS
   76: *> \verbatim
   77: *>          NRHS is INTEGER
   78: *>          The number of right hand sides, i.e., the number of columns
   79: *>          of the matrix B.  NRHS >= 0.
   80: *> \endverbatim
   81: *>
   82: *> \param[in,out] A
   83: *> \verbatim
   84: *>          A is COMPLEX*16 array, dimension (LDA,N)
   85: *>          On entry, the Hermitian matrix A.
   86: *>            If UPLO = 'U': the leading N-by-N upper triangular part
   87: *>            of A contains the upper triangular part of the matrix A,
   88: *>            and the strictly lower triangular part of A is not
   89: *>            referenced.
   90: *>
   91: *>            If UPLO = 'L': the leading N-by-N lower triangular part
   92: *>            of A contains the lower triangular part of the matrix A,
   93: *>            and the strictly upper triangular part of A is not
   94: *>            referenced.
   95: *>
   96: *>          On exit, if INFO = 0, diagonal of the block diagonal
   97: *>          matrix D and factors U or L  as computed by ZHETRF_RK:
   98: *>            a) ONLY diagonal elements of the Hermitian block diagonal
   99: *>               matrix D on the diagonal of A, i.e. D(k,k) = A(k,k);
  100: *>               (superdiagonal (or subdiagonal) elements of D
  101: *>                are stored on exit in array E), and
  102: *>            b) If UPLO = 'U': factor U in the superdiagonal part of A.
  103: *>               If UPLO = 'L': factor L in the subdiagonal part of A.
  104: *>
  105: *>          For more info see the description of ZHETRF_RK routine.
  106: *> \endverbatim
  107: *>
  108: *> \param[in] LDA
  109: *> \verbatim
  110: *>          LDA is INTEGER
  111: *>          The leading dimension of the array A.  LDA >= max(1,N).
  112: *> \endverbatim
  113: *>
  114: *> \param[out] E
  115: *> \verbatim
  116: *>          E is COMPLEX*16 array, dimension (N)
  117: *>          On exit, contains the output computed by the factorization
  118: *>          routine ZHETRF_RK, i.e. the superdiagonal (or subdiagonal)
  119: *>          elements of the Hermitian block diagonal matrix D
  120: *>          with 1-by-1 or 2-by-2 diagonal blocks, where
  121: *>          If UPLO = 'U': E(i) = D(i-1,i), i=2:N, E(1) is set to 0;
  122: *>          If UPLO = 'L': E(i) = D(i+1,i), i=1:N-1, E(N) is set to 0.
  123: *>
  124: *>          NOTE: For 1-by-1 diagonal block D(k), where
  125: *>          1 <= k <= N, the element E(k) is set to 0 in both
  126: *>          UPLO = 'U' or UPLO = 'L' cases.
  127: *>
  128: *>          For more info see the description of ZHETRF_RK routine.
  129: *> \endverbatim
  130: *>
  131: *> \param[out] IPIV
  132: *> \verbatim
  133: *>          IPIV is INTEGER array, dimension (N)
  134: *>          Details of the interchanges and the block structure of D,
  135: *>          as determined by ZHETRF_RK.
  136: *>
  137: *>          For more info see the description of ZHETRF_RK routine.
  138: *> \endverbatim
  139: *>
  140: *> \param[in,out] B
  141: *> \verbatim
  142: *>          B is COMPLEX*16 array, dimension (LDB,NRHS)
  143: *>          On entry, the N-by-NRHS right hand side matrix B.
  144: *>          On exit, if INFO = 0, the N-by-NRHS solution matrix X.
  145: *> \endverbatim
  146: *>
  147: *> \param[in] LDB
  148: *> \verbatim
  149: *>          LDB is INTEGER
  150: *>          The leading dimension of the array B.  LDB >= max(1,N).
  151: *> \endverbatim
  152: *>
  153: *> \param[out] WORK
  154: *> \verbatim
  155: *>          WORK is COMPLEX*16 array, dimension ( MAX(1,LWORK) ).
  156: *>          Work array used in the factorization stage.
  157: *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
  158: *> \endverbatim
  159: *>
  160: *> \param[in] LWORK
  161: *> \verbatim
  162: *>          LWORK is INTEGER
  163: *>          The length of WORK.  LWORK >= 1. For best performance
  164: *>          of factorization stage LWORK >= max(1,N*NB), where NB is
  165: *>          the optimal blocksize for ZHETRF_RK.
  166: *>
  167: *>          If LWORK = -1, then a workspace query is assumed;
  168: *>          the routine only calculates the optimal size of the WORK
  169: *>          array for factorization stage, returns this value as
  170: *>          the first entry of the WORK array, and no error message
  171: *>          related to LWORK is issued by XERBLA.
  172: *> \endverbatim
  173: *>
  174: *> \param[out] INFO
  175: *> \verbatim
  176: *>          INFO is INTEGER
  177: *>          = 0: successful exit
  178: *>
  179: *>          < 0: If INFO = -k, the k-th argument had an illegal value
  180: *>
  181: *>          > 0: If INFO = k, the matrix A is singular, because:
  182: *>                 If UPLO = 'U': column k in the upper
  183: *>                 triangular part of A contains all zeros.
  184: *>                 If UPLO = 'L': column k in the lower
  185: *>                 triangular part of A contains all zeros.
  186: *>
  187: *>               Therefore D(k,k) is exactly zero, and superdiagonal
  188: *>               elements of column k of U (or subdiagonal elements of
  189: *>               column k of L ) are all zeros. The factorization has
  190: *>               been completed, but the block diagonal matrix D is
  191: *>               exactly singular, and division by zero will occur if
  192: *>               it is used to solve a system of equations.
  193: *>
  194: *>               NOTE: INFO only stores the first occurrence of
  195: *>               a singularity, any subsequent occurrence of singularity
  196: *>               is not stored in INFO even though the factorization
  197: *>               always completes.
  198: *> \endverbatim
  199: *
  200: *  Authors:
  201: *  ========
  202: *
  203: *> \author Univ. of Tennessee
  204: *> \author Univ. of California Berkeley
  205: *> \author Univ. of Colorado Denver
  206: *> \author NAG Ltd.
  207: *
  208: *> \ingroup complex16HEsolve
  209: *
  210: *> \par Contributors:
  211: *  ==================
  212: *>
  213: *> \verbatim
  214: *>
  215: *>  December 2016,  Igor Kozachenko,
  216: *>                  Computer Science Division,
  217: *>                  University of California, Berkeley
  218: *>
  219: *>  September 2007, Sven Hammarling, Nicholas J. Higham, Craig Lucas,
  220: *>                  School of Mathematics,
  221: *>                  University of Manchester
  222: *>
  223: *> \endverbatim
  224: *
  225: *  =====================================================================
  226:       SUBROUTINE ZHESV_RK( UPLO, N, NRHS, A, LDA, E, IPIV, B, LDB, WORK,
  227:      $                     LWORK, INFO )
  228: *
  229: *  -- LAPACK driver routine --
  230: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  231: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  232: *
  233: *     .. Scalar Arguments ..
  234:       CHARACTER          UPLO
  235:       INTEGER            INFO, LDA, LDB, LWORK, N, NRHS
  236: *     ..
  237: *     .. Array Arguments ..
  238:       INTEGER            IPIV( * )
  239:       COMPLEX*16         A( LDA, * ), B( LDB, * ), E( * ), WORK( * )
  240: *     ..
  241: *
  242: *  =====================================================================
  243: *
  244: *     .. Local Scalars ..
  245:       LOGICAL            LQUERY
  246:       INTEGER            LWKOPT
  247: *     ..
  248: *     .. External Functions ..
  249:       LOGICAL            LSAME
  250:       EXTERNAL           LSAME
  251: *     ..
  252: *     .. External Subroutines ..
  253:       EXTERNAL           XERBLA, ZHETRF_RK, ZHETRS_3
  254: *     ..
  255: *     .. Intrinsic Functions ..
  256:       INTRINSIC          MAX
  257: *     ..
  258: *     .. Executable Statements ..
  259: *
  260: *     Test the input parameters.
  261: *
  262:       INFO = 0
  263:       LQUERY = ( LWORK.EQ.-1 )
  264:       IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
  265:          INFO = -1
  266:       ELSE IF( N.LT.0 ) THEN
  267:          INFO = -2
  268:       ELSE IF( NRHS.LT.0 ) THEN
  269:          INFO = -3
  270:       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
  271:          INFO = -5
  272:       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
  273:          INFO = -9
  274:       ELSE IF( LWORK.LT.1 .AND. .NOT.LQUERY ) THEN
  275:          INFO = -11
  276:       END IF
  277: *
  278:       IF( INFO.EQ.0 ) THEN
  279:          IF( N.EQ.0 ) THEN
  280:             LWKOPT = 1
  281:          ELSE
  282:             CALL ZHETRF_RK( UPLO, N, A, LDA, E, IPIV, WORK, -1, INFO )
  283:             LWKOPT = INT( DBLE( WORK( 1 ) ) )
  284:          END IF
  285:          WORK( 1 ) = LWKOPT
  286:       END IF
  287: *
  288:       IF( INFO.NE.0 ) THEN
  289:          CALL XERBLA( 'ZHESV_RK ', -INFO )
  290:          RETURN
  291:       ELSE IF( LQUERY ) THEN
  292:          RETURN
  293:       END IF
  294: *
  295: *     Compute the factorization A = P*U*D*(U**H)*(P**T) or
  296: *     A = P*U*D*(U**H)*(P**T).
  297: *
  298:       CALL ZHETRF_RK( UPLO, N, A, LDA, E, IPIV, WORK, LWORK, INFO )
  299: *
  300:       IF( INFO.EQ.0 ) THEN
  301: *
  302: *        Solve the system A*X = B with BLAS3 solver, overwriting B with X.
  303: *
  304:          CALL ZHETRS_3( UPLO, N, NRHS, A, LDA, E, IPIV, B, LDB, INFO )
  305: *
  306:       END IF
  307: *
  308:       WORK( 1 ) = LWKOPT
  309: *
  310:       RETURN
  311: *
  312: *     End of ZHESV_RK
  313: *
  314:       END

CVSweb interface <joel.bertrand@systella.fr>