Return to zhesv_aa_2stage.f CVS log

Up to [local] / rpl / lapack / lapack

Première mise à jour de lapack et blas.

    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**H * T * U,  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: *>          LTB is INTEGER
  110: *>          The size of the array TB. LTB >= 4*N, internally
  111: *>          used to select NB such that LTB >= (3*NB+1)*N.
  112: *>
  113: *>          If LTB = -1, then a workspace query is assumed; the
  114: *>          routine only calculates the optimal size of LTB, 
  115: *>          returns this value as the first entry of TB, and
  116: *>          no error message related to LTB is issued by XERBLA.
  117: *> \endverbatim
  118: *>
  119: *> \param[out] IPIV
  120: *> \verbatim
  121: *>          IPIV is INTEGER array, dimension (N)
  122: *>          On exit, it contains the details of the interchanges, i.e.,
  123: *>          the row and column k of A were interchanged with the
  124: *>          row and column IPIV(k).
  125: *> \endverbatim
  126: *>
  127: *> \param[out] IPIV2
  128: *> \verbatim
  129: *>          IPIV2 is INTEGER array, dimension (N)
  130: *>          On exit, it contains the details of the interchanges, i.e.,
  131: *>          the row and column k of T were interchanged with the
  132: *>          row and column IPIV(k).
  133: *> \endverbatim
  134: *>
  135: *> \param[in,out] B
  136: *> \verbatim
  137: *>          B is COMPLEX*16 array, dimension (LDB,NRHS)
  138: *>          On entry, the right hand side matrix B.
  139: *>          On exit, the solution matrix X.
  140: *> \endverbatim
  141: *>
  142: *> \param[in] LDB
  143: *> \verbatim
  144: *>          LDB is INTEGER
  145: *>          The leading dimension of the array B.  LDB >= max(1,N).
  146: *> \endverbatim
  147: *>
  148: *> \param[out] WORK
  149: *> \verbatim
  150: *>          WORK is COMPLEX*16 workspace of size LWORK
  151: *> \endverbatim
  152: *>
  153: *> \param[in] LWORK
  154: *> \verbatim
  155: *>          LWORK is INTEGER
  156: *>          The size of WORK. LWORK >= N, internally used to select NB
  157: *>          such that LWORK >= N*NB.
  158: *>
  159: *>          If LWORK = -1, then a workspace query is assumed; the
  160: *>          routine only calculates the optimal size of the WORK array,
  161: *>          returns this value as the first entry of the WORK array, and
  162: *>          no error message related to LWORK is issued by XERBLA.
  163: *> \endverbatim
  164: *>
  165: *> \param[out] INFO
  166: *> \verbatim
  167: *>          INFO is INTEGER
  168: *>          = 0:  successful exit
  169: *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
  170: *>          > 0:  if INFO = i, band LU factorization failed on i-th column
  171: *> \endverbatim
  172: *
  173: *  Authors:
  174: *  ========
  175: *
  176: *> \author Univ. of Tennessee
  177: *> \author Univ. of California Berkeley
  178: *> \author Univ. of Colorado Denver
  179: *> \author NAG Ltd.
  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 --
  189: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  190: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  191: *
  192:       IMPLICIT NONE
  193: *
  194: *     .. Scalar Arguments ..
  195:       CHARACTER          UPLO
  196:       INTEGER            N, NRHS, LDA, LDB, LTB, LWORK, INFO
  197: *     ..
  198: *     .. Array Arguments ..
  199:       INTEGER            IPIV( * ), IPIV2( * )
  200:       COMPLEX*16         A( LDA, * ), B( LDB, * ), TB( * ), WORK( * )
  201: *     ..
  202: *
  203: *  =====================================================================
  204: *     .. Parameters ..
  205:       COMPLEX*16         ZERO, ONE
  206:       PARAMETER          ( ZERO = ( 0.0D+0, 0.0D+0 ),
  207:      $                     ONE  = ( 1.0D+0, 0.0D+0 ) )
  208: *
  209: *     .. Local Scalars ..
  210:       LOGICAL            UPPER, TQUERY, WQUERY
  211:       INTEGER            LWKOPT
  212: *     ..
  213: *     .. External Functions ..
  214:       LOGICAL            LSAME
  215:       INTEGER            ILAENV
  216:       EXTERNAL           LSAME, ILAENV
  217: *     ..
  218: *     .. External Subroutines ..
  219:       EXTERNAL           XERBLA, ZHETRF_AA_2STAGE, ZHETRS_AA_2STAGE
  220: *     ..
  221: *     .. Intrinsic Functions ..
  222:       INTRINSIC          MAX
  223: *     ..
  224: *     .. Executable Statements ..
  225: *
  226: *     Test the input parameters.
  227: *
  228:       INFO = 0
  229:       UPPER = LSAME( UPLO, 'U' )
  230:       WQUERY = ( LWORK.EQ.-1 )
  231:       TQUERY = ( LTB.EQ.-1 )
  232:       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
  233:          INFO = -1
  234:       ELSE IF( N.LT.0 ) THEN
  235:          INFO = -2
  236:       ELSE IF( NRHS.LT.0 ) THEN
  237:          INFO = -3
  238:       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
  239:          INFO = -5
  240:       ELSE IF( LTB.LT.( 4*N ) .AND. .NOT.TQUERY ) THEN
  241:          INFO = -7
  242:       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
  243:          INFO = -11
  244:       ELSE IF( LWORK.LT.N .AND. .NOT.WQUERY ) THEN
  245:          INFO = -13
  246:       END IF
  247: *
  248:       IF( INFO.EQ.0 ) THEN
  249:          CALL ZHETRF_AA_2STAGE( UPLO, N, A, LDA, TB, -1, IPIV,
  250:      $                          IPIV2, WORK, -1, INFO )
  251:          LWKOPT = INT( WORK(1) )
  252:       END IF
  253: *
  254:       IF( INFO.NE.0 ) THEN
  255:          CALL XERBLA( 'ZHESV_AA_2STAGE', -INFO )
  256:          RETURN
  257:       ELSE IF( WQUERY .OR. TQUERY ) THEN
  258:          RETURN
  259:       END IF
  260: *
  261: *     Compute the factorization A = U**H*T*U or A = L*T*L**H.
  262: *
  263:       CALL ZHETRF_AA_2STAGE( UPLO, N, A, LDA, TB, LTB, IPIV, IPIV2,
  264:      $                       WORK, LWORK, INFO )
  265:       IF( INFO.EQ.0 ) THEN
  266: *
  267: *        Solve the system A*X = B, overwriting B with X.
  268: *
  269:          CALL ZHETRS_AA_2STAGE( UPLO, N, NRHS, A, LDA, TB, LTB, IPIV,
  270:      $                          IPIV2, B, LDB, INFO )
  271: *
  272:       END IF
  273: *
  274:       WORK( 1 ) = LWKOPT
  275: *
  276:       RETURN
  277: *
  278: *     End of ZHESV_AA_2STAGE
  279: *
  280:       END