Return to zsysv_aa_2stage.f CVS log

Up to [local] / rpl / lapack / lapack

Première mise à jour de lapack et blas.

    1: *> \brief <b> ZSYSV_AA_2STAGE 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 ZSYSV_AA_2STAGE + dependencies
   10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zsysv_aa_2stage.f">
   11: *> [TGZ]</a>
   12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zsysv_aa_2stage.f">
   13: *> [ZIP]</a>
   14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zsysv_aa_2stage.f">
   15: *> [TXT]</a>
   16: *> \endhtmlonly
   17: *
   18: *  Definition:
   19: *  ===========
   20: *
   21: *      SUBROUTINE ZSYSV_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: *> \par Purpose:
   35: *  =============
   36: *>
   37: *> \verbatim
   38: *>
   39: *> ZSYSV_AA_2STAGE computes the solution to a complex system of
   40: *> linear equations
   41: *>    A * X = B,
   42: *> where A is an N-by-N symmetric matrix and X and B are N-by-NRHS
   43: *> matrices.
   44: *>
   45: *> Aasen's 2-stage algorithm is used to factor A as
   46: *>    A = U**T * T * U,  if UPLO = 'U', or
   47: *>    A = L * T * L**T,  if UPLO = 'L',
   48: *> where U (or L) is a product of permutation and unit upper (lower)
   49: *> triangular matrices, and T is symmetric and band. The matrix T is
   50: *> then LU-factored with partial pivoting. The factored form of A
   51: *> is then used to solve the system of equations A * X = B.
   52: *>
   53: *> This is the blocked version of the algorithm, calling Level 3 BLAS.
   54: *> \endverbatim
   55: *
   56: *  Arguments:
   57: *  ==========
   58: *
   59: *> \param[in] UPLO
   60: *> \verbatim
   61: *>          UPLO is CHARACTER*1
   62: *>          = 'U':  Upper triangle of A is stored;
   63: *>          = 'L':  Lower triangle of A is stored.
   64: *> \endverbatim
   65: *>
   66: *> \param[in] N
   67: *> \verbatim
   68: *>          N is INTEGER
   69: *>          The order of the matrix A.  N >= 0.
   70: *> \endverbatim
   71: *>
   72: *> \param[in] NRHS
   73: *> \verbatim
   74: *>          NRHS is INTEGER
   75: *>          The number of right hand sides, i.e., the number of columns
   76: *>          of the matrix B.  NRHS >= 0.
   77: *> \endverbatim
   78: *>
   79: *> \param[in,out] A
   80: *> \verbatim
   81: *>          A is COMPLEX*16 array, dimension (LDA,N)
   82: *>          On entry, the symmetric matrix A.  If UPLO = 'U', the leading
   83: *>          N-by-N upper triangular part of A contains the upper
   84: *>          triangular part of the matrix A, and the strictly lower
   85: *>          triangular part of A is not referenced.  If UPLO = 'L', the
   86: *>          leading N-by-N lower triangular part of A contains the lower
   87: *>          triangular part of the matrix A, and the strictly upper
   88: *>          triangular part of A is not referenced.
   89: *>
   90: *>          On exit, L is stored below (or above) the subdiaonal blocks,
   91: *>          when UPLO  is 'L' (or 'U').
   92: *> \endverbatim
   93: *>
   94: *> \param[in] LDA
   95: *> \verbatim
   96: *>          LDA is INTEGER
   97: *>          The leading dimension of the array A.  LDA >= max(1,N).
   98: *> \endverbatim
   99: *>
  100: *> \param[out] TB
  101: *> \verbatim
  102: *>          TB is COMPLEX*16 array, dimension (LTB)
  103: *>          On exit, details of the LU factorization of the band matrix.
  104: *> \endverbatim
  105: *>
  106: *> \param[in] LTB
  107: *> \verbatim
  108: *>          LTB is INTEGER
  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: *>          IPIV2 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: *>          LWORK is INTEGER
  155: *>          The size of WORK. LWORK >= N, internally used to select NB
  156: *>          such that LWORK >= N*NB.
  157: *>
  158: *>          If LWORK = -1, then a workspace query is assumed; the
  159: *>          routine only calculates the optimal size of the WORK array,
  160: *>          returns this value as the first entry of the WORK array, and
  161: *>          no error message related to LWORK is issued by XERBLA.
  162: *> \endverbatim
  163: *>
  164: *> \param[out] INFO
  165: *> \verbatim
  166: *>          INFO is INTEGER
  167: *>          = 0:  successful exit
  168: *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
  169: *>          > 0:  if INFO = i, band LU factorization failed on i-th column
  170: *> \endverbatim
  171: *
  172: *  Authors:
  173: *  ========
  174: *
  175: *> \author Univ. of Tennessee
  176: *> \author Univ. of California Berkeley
  177: *> \author Univ. of Colorado Denver
  178: *> \author NAG Ltd.
  179: *
  180: *> \ingroup complex16SYcomputational
  181: *
  182: *  =====================================================================
  183:       SUBROUTINE ZSYSV_AA_2STAGE( UPLO, N, NRHS, A, LDA, TB, LTB,
  184:      $                            IPIV, IPIV2, B, LDB, WORK, LWORK,
  185:      $                            INFO )
  186: *
  187: *  -- LAPACK computational routine --
  188: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  189: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  190: *
  191:       IMPLICIT NONE
  192: *
  193: *     .. Scalar Arguments ..
  194:       CHARACTER          UPLO
  195:       INTEGER            N, NRHS, LDA, LDB, LTB, LWORK, INFO
  196: *     ..
  197: *     .. Array Arguments ..
  198:       INTEGER            IPIV( * ), IPIV2( * )
  199:       COMPLEX*16         A( LDA, * ), B( LDB, * ), TB( * ), WORK( * )
  200: *     ..
  201: *
  202: *  =====================================================================
  203: *
  204: *     .. Local Scalars ..
  205:       LOGICAL            UPPER, TQUERY, WQUERY
  206:       INTEGER            LWKOPT
  207: *     ..
  208: *     .. External Functions ..
  209:       LOGICAL            LSAME
  210:       EXTERNAL           LSAME
  211: *     ..
  212: *     .. External Subroutines ..
  213:       EXTERNAL           XERBLA, ZSYTRF_AA_2STAGE,
  214:      $                   ZSYTRS_AA_2STAGE
  215: *     ..
  216: *     .. Intrinsic Functions ..
  217:       INTRINSIC          MAX
  218: *     ..
  219: *     .. Executable Statements ..
  220: *
  221: *     Test the input parameters.
  222: *
  223:       INFO = 0
  224:       UPPER = LSAME( UPLO, 'U' )
  225:       WQUERY = ( LWORK.EQ.-1 )
  226:       TQUERY = ( LTB.EQ.-1 )
  227:       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
  228:          INFO = -1
  229:       ELSE IF( N.LT.0 ) THEN
  230:          INFO = -2
  231:       ELSE IF( NRHS.LT.0 ) THEN
  232:          INFO = -3
  233:       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
  234:          INFO = -5
  235:       ELSE IF( LTB.LT.( 4*N ) .AND. .NOT.TQUERY ) THEN
  236:          INFO = -7
  237:       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
  238:          INFO = -11
  239:       ELSE IF( LWORK.LT.N .AND. .NOT.WQUERY ) THEN
  240:          INFO = -13
  241:       END IF
  242: *
  243:       IF( INFO.EQ.0 ) THEN
  244:          CALL ZSYTRF_AA_2STAGE( UPLO, N, A, LDA, TB, -1, IPIV,
  245:      $                          IPIV2, WORK, -1, INFO )
  246:          LWKOPT = INT( WORK(1) )
  247:       END IF
  248: *
  249:       IF( INFO.NE.0 ) THEN
  250:          CALL XERBLA( 'ZSYSV_AA_2STAGE', -INFO )
  251:          RETURN
  252:       ELSE IF( WQUERY .OR. TQUERY ) THEN
  253:          RETURN
  254:       END IF
  255: *
  256: *
  257: *     Compute the factorization A = U**T*T*U or A = L*T*L**T.
  258: *
  259:       CALL ZSYTRF_AA_2STAGE( UPLO, N, A, LDA, TB, LTB, IPIV, IPIV2,
  260:      $                       WORK, LWORK, INFO )
  261:       IF( INFO.EQ.0 ) THEN
  262: *
  263: *        Solve the system A*X = B, overwriting B with X.
  264: *
  265:          CALL ZSYTRS_AA_2STAGE( UPLO, N, NRHS, A, LDA, TB, LTB, IPIV,
  266:      $                          IPIV2, B, LDB, INFO )
  267: *
  268:       END IF
  269: *
  270:       WORK( 1 ) = LWKOPT
  271: *
  272: *     End of ZSYSV_AA_2STAGE
  273: *
  274:       END