Return to zhesv_aa_2stage.f CVS log

Up to [local] / rpl / lapack / lapack

Ajout d'un fichier oublié.

    1: *> \brief <b> ZHESV_AA_2STAGE computes the solution to system of linear equations A * X = B for HE 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_AA_2STAGE + dependencies
   10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhesv_aa_2stage.f">
   11: *> [TGZ]</a>
   12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhesv_aa_2stage.f">
   13: *> [ZIP]</a>
   14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhesv_aa_2stage.f">
   15: *> [TXT]</a>
   16: *> \endhtmlonly
   17: *
   18: *  Definition:
   19: *  ===========
   20: *
   21: *      SUBROUTINE ZHESV_AA_2STAGE( UPLO, N, NRHS, A, LDA, TB, LTB,
   22: *                                  IPIV, IPIV2, B, LDB, WORK, LWORK,
   23: *                                  INFO )
   24: *
   25: *       .. Scalar Arguments ..
   26: *       CHARACTER          UPLO
   27: *       INTEGER            N, NRHS, LDA, LTB, LDB, LWORK, INFO
   28: *       ..
   29: *       .. Array Arguments ..
   30: *       INTEGER            IPIV( * ), IPIV2( * )
   31: *       COMPLEX*16         A( LDA, * ), TB( * ), B( LDB, *), WORK( * )
   32: *       ..
   33: *
   34: *
   35: *> \par Purpose:
   36: *  =============
   37: *>
   38: *> \verbatim
   39: *>
   40: *> ZHESV_AA_2STAGE computes the solution to a complex system of 
   41: *> linear equations
   42: *>    A * X = B,
   43: *> where A is an N-by-N Hermitian matrix and X and B are N-by-NRHS
   44: *> matrices.
   45: *>
   46: *> Aasen's 2-stage algorithm is used to factor A as
   47: *>    A = U * T * U**H,  if UPLO = 'U', or
   48: *>    A = L * T * L**H,  if UPLO = 'L',
   49: *> where U (or L) is a product of permutation and unit upper (lower)
   50: *> triangular matrices, and T is Hermitian and band. The matrix T is
   51: *> then LU-factored with partial pivoting. The factored form of A
   52: *> is then used to solve the system of equations A * X = B.
   53: *>
   54: *> This is the blocked version of the algorithm, calling Level 3 BLAS.
   55: *> \endverbatim
   56: *
   57: *  Arguments:
   58: *  ==========
   59: *
   60: *> \param[in] UPLO
   61: *> \verbatim
   62: *>          UPLO is CHARACTER*1
   63: *>          = 'U':  Upper triangle of A is stored;
   64: *>          = 'L':  Lower triangle of A is stored.
   65: *> \endverbatim
   66: *>
   67: *> \param[in] N
   68: *> \verbatim
   69: *>          N is INTEGER
   70: *>          The order of the matrix A.  N >= 0.
   71: *> \endverbatim
   72: *>
   73: *> \param[in] NRHS
   74: *> \verbatim
   75: *>          NRHS is INTEGER
   76: *>          The number of right hand sides, i.e., the number of columns
   77: *>          of the matrix B.  NRHS >= 0.
   78: *> \endverbatim
   79: *>
   80: *> \param[in,out] A
   81: *> \verbatim
   82: *>          A is COMPLEX*16 array, dimension (LDA,N)
   83: *>          On entry, the hermitian matrix A.  If UPLO = 'U', the leading
   84: *>          N-by-N upper triangular part of A contains the upper
   85: *>          triangular part of the matrix A, and the strictly lower
   86: *>          triangular part of A is not referenced.  If UPLO = 'L', the
   87: *>          leading N-by-N lower triangular part of A contains the lower
   88: *>          triangular part of the matrix A, and the strictly upper
   89: *>          triangular part of A is not referenced.
   90: *>
   91: *>          On exit, L is stored below (or above) the subdiaonal blocks,
   92: *>          when UPLO  is 'L' (or 'U').
   93: *> \endverbatim
   94: *>
   95: *> \param[in] LDA
   96: *> \verbatim
   97: *>          LDA is INTEGER
   98: *>          The leading dimension of the array A.  LDA >= max(1,N).
   99: *> \endverbatim
  100: *>
  101: *> \param[out] TB
  102: *> \verbatim
  103: *>          TB is COMPLEX*16 array, dimension (LTB)
  104: *>          On exit, details of the LU factorization of the band matrix.
  105: *> \endverbatim
  106: *>
  107: *> \param[in] LTB
  108: *> \verbatim
  109: *>          The size of the array TB. LTB >= 4*N, internally
  110: *>          used to select NB such that LTB >= (3*NB+1)*N.
  111: *>
  112: *>          If LTB = -1, then a workspace query is assumed; the
  113: *>          routine only calculates the optimal size of LTB, 
  114: *>          returns this value as the first entry of TB, and
  115: *>          no error message related to LTB is issued by XERBLA.
  116: *> \endverbatim
  117: *>
  118: *> \param[out] IPIV
  119: *> \verbatim
  120: *>          IPIV is INTEGER array, dimension (N)
  121: *>          On exit, it contains the details of the interchanges, i.e.,
  122: *>          the row and column k of A were interchanged with the
  123: *>          row and column IPIV(k).
  124: *> \endverbatim
  125: *>
  126: *> \param[out] IPIV2
  127: *> \verbatim
  128: *>          IPIV is INTEGER array, dimension (N)
  129: *>          On exit, it contains the details of the interchanges, i.e.,
  130: *>          the row and column k of T were interchanged with the
  131: *>          row and column IPIV(k).
  132: *> \endverbatim
  133: *>
  134: *> \param[in,out] B
  135: *> \verbatim
  136: *>          B is COMPLEX*16 array, dimension (LDB,NRHS)
  137: *>          On entry, the right hand side matrix B.
  138: *>          On exit, the solution matrix X.
  139: *> \endverbatim
  140: *>
  141: *> \param[in] LDB
  142: *> \verbatim
  143: *>          LDB is INTEGER
  144: *>          The leading dimension of the array B.  LDB >= max(1,N).
  145: *> \endverbatim
  146: *>
  147: *> \param[out] WORK
  148: *> \verbatim
  149: *>          WORK is COMPLEX*16 workspace of size LWORK
  150: *> \endverbatim
  151: *>
  152: *> \param[in] LWORK
  153: *> \verbatim
  154: *>          The size of WORK. LWORK >= N, internally used to select NB
  155: *>          such that LWORK >= N*NB.
  156: *>
  157: *>          If LWORK = -1, then a workspace query is assumed; the
  158: *>          routine only calculates the optimal size of the WORK array,
  159: *>          returns this value as the first entry of the WORK array, and
  160: *>          no error message related to LWORK is issued by XERBLA.
  161: *> \endverbatim
  162: *>
  163: *> \param[out] INFO
  164: *> \verbatim
  165: *>          INFO is INTEGER
  166: *>          = 0:  successful exit
  167: *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
  168: *>          > 0:  if INFO = i, band LU factorization failed on i-th column
  169: *> \endverbatim
  170: *
  171: *  Authors:
  172: *  ========
  173: *
  174: *> \author Univ. of Tennessee
  175: *> \author Univ. of California Berkeley
  176: *> \author Univ. of Colorado Denver
  177: *> \author NAG Ltd.
  178: *
  179: *> \date November 2017
  180: *
  181: *> \ingroup complex16HEsolve
  182: *
  183: *  =====================================================================
  184:       SUBROUTINE ZHESV_AA_2STAGE( UPLO, N, NRHS, A, LDA, TB, LTB,
  185:      $                            IPIV, IPIV2, B, LDB, WORK, LWORK,
  186:      $                            INFO )
  187: *
  188: *  -- LAPACK driver routine (version 3.8.0) --
  189: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  190: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  191: *     November 2017
  192: *
  193:       IMPLICIT NONE
  194: *
  195: *     .. Scalar Arguments ..
  196:       CHARACTER          UPLO
  197:       INTEGER            N, NRHS, LDA, LDB, LTB, LWORK, INFO
  198: *     ..
  199: *     .. Array Arguments ..
  200:       INTEGER            IPIV( * ), IPIV2( * )
  201:       COMPLEX*16         A( LDA, * ), B( LDB, * ), TB( * ), WORK( * )
  202: *     ..
  203: *
  204: *  =====================================================================
  205: *     .. Parameters ..
  206:       COMPLEX*16         ZERO, ONE
  207:       PARAMETER          ( ZERO = ( 0.0D+0, 0.0D+0 ),
  208:      $                     ONE  = ( 1.0D+0, 0.0D+0 ) )
  209: *
  210: *     .. Local Scalars ..
  211:       LOGICAL            UPPER, TQUERY, WQUERY
  212:       INTEGER            I, J, K, I1, I2, TD
  213:       INTEGER            LDTB, LWKOPT, NB, KB, NT, IINFO
  214:       COMPLEX            PIV
  215: *     ..
  216: *     .. External Functions ..
  217:       LOGICAL            LSAME
  218:       INTEGER            ILAENV
  219:       EXTERNAL           LSAME, ILAENV
  220: *     ..
  221: *     .. External Subroutines ..
  222:       EXTERNAL           XERBLA, ZHETRF_AA_2STAGE, ZHETRS_AA_2STAGE
  223: *     ..
  224: *     .. Intrinsic Functions ..
  225:       INTRINSIC          MAX
  226: *     ..
  227: *     .. Executable Statements ..
  228: *
  229: *     Test the input parameters.
  230: *
  231:       INFO = 0
  232:       UPPER = LSAME( UPLO, 'U' )
  233:       WQUERY = ( LWORK.EQ.-1 )
  234:       TQUERY = ( LTB.EQ.-1 )
  235:       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
  236:          INFO = -1
  237:       ELSE IF( N.LT.0 ) THEN
  238:          INFO = -2
  239:       ELSE IF( NRHS.LT.0 ) THEN
  240:          INFO = -3
  241:       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
  242:          INFO = -5
  243:       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
  244:          INFO = -11
  245:       END IF
  246: *
  247:       IF( INFO.EQ.0 ) THEN
  248:          CALL ZHETRF_AA_2STAGE( UPLO, N, A, LDA, TB, -1, IPIV,
  249:      $                          IPIV2, WORK, -1, INFO )
  250:          LWKOPT = INT( WORK(1) )
  251:          IF( LTB.LT.INT( TB(1) ) .AND. .NOT.TQUERY ) THEN
  252:             INFO = -7
  253:          ELSE IF( LWORK.LT.LWKOPT .AND. .NOT.WQUERY ) THEN
  254:             INFO = -13
  255:          END IF
  256:       END IF
  257: *
  258:       IF( INFO.NE.0 ) THEN
  259:          CALL XERBLA( 'ZHESV_AA_2STAGE', -INFO )
  260:          RETURN
  261:       ELSE IF( WQUERY .OR. TQUERY ) THEN
  262:          RETURN
  263:       END IF
  264: *
  265: *     Compute the factorization A = U*T*U**H or A = L*T*L**H.
  266: *
  267:       CALL ZHETRF_AA_2STAGE( UPLO, N, A, LDA, TB, LTB, IPIV, IPIV2,
  268:      $                       WORK, LWORK, INFO )
  269:       IF( INFO.EQ.0 ) THEN
  270: *
  271: *        Solve the system A*X = B, overwriting B with X.
  272: *
  273:          CALL ZHETRS_AA_2STAGE( UPLO, N, NRHS, A, LDA, TB, LTB, IPIV,
  274:      $                          IPIV2, B, LDB, INFO )
  275: *
  276:       END IF
  277: *
  278:       WORK( 1 ) = LWKOPT
  279: *
  280:       RETURN
  281: *
  282: *     End of ZHESV_AA_2STAGE
  283: *
  284:       END