SRC/zhesv_aasen.f

   1 *> \brief <b> ZHESV_AASEN 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_AASEN + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhesv_aasen.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhesv_aasen.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhesv_aasen.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE ZHESV_AASEN( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, WORK,
  22 *                               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, * ), WORK( * )
  31 *       ..
  32 *
  33 *
  34 *> \par Purpose:
  35 *  =============
  36 *>
  37 *> \verbatim
  38 *>
  39 *> ZHESV_AASEN computes the solution to a complex system of linear equations
  40 *>    A * X = B,
  41 *> where A is an N-by-N Hermitian matrix and X and B are N-by-NRHS
  42 *> matrices.
  43 *>
  44 *> Aasen's algorithm is used to factor A as
  45 *>    A = U * T * U**H,  if UPLO = 'U', or
  46 *>    A = L * T * L**H,  if UPLO = 'L',
  47 *> where U (or L) is a product of permutation and unit upper (lower)
  48 *> triangular matrices, and T is Hermitian and tridiagonal. The factored form
  49 *> of A is then used to solve the system of equations A * X = B.
  50 *> \endverbatim
  51 *
  52 *  Arguments:
  53 *  ==========
  54 *
  55 *> \param[in] UPLO
  56 *> \verbatim
  57 *>          UPLO is CHARACTER*1
  58 *>          = 'U':  Upper triangle of A is stored;
  59 *>          = 'L':  Lower triangle of A is stored.
  60 *> \endverbatim
  61 *>
  62 *> \param[in] N
  63 *> \verbatim
  64 *>          N is INTEGER
  65 *>          The number of linear equations, i.e., the order of the
  66 *>          matrix A.  N >= 0.
  67 *> \endverbatim
  68 *>
  69 *> \param[in] NRHS
  70 *> \verbatim
  71 *>          NRHS is INTEGER
  72 *>          The number of right hand sides, i.e., the number of columns
  73 *>          of the matrix B.  NRHS >= 0.
  74 *> \endverbatim
  75 *>
  76 *> \param[in,out] A
  77 *> \verbatim
  78 *>          A is COMPLEX*16 array, dimension (LDA,N)
  79 *>          On entry, the Hermitian matrix A.  If UPLO = 'U', the leading
  80 *>          N-by-N upper triangular part of A contains the upper
  81 *>          triangular part of the matrix A, and the strictly lower
  82 *>          triangular part of A is not referenced.  If UPLO = 'L', the
  83 *>          leading N-by-N lower triangular part of A contains the lower
  84 *>          triangular part of the matrix A, and the strictly upper
  85 *>          triangular part of A is not referenced.
  86 *>
  87 *>          On exit, if INFO = 0, the tridiagonal matrix T and the
  88 *>          multipliers used to obtain the factor U or L from the
  89 *>          factorization A = U*T*U**H or A = L*T*L**H as computed by
  90 *>          ZHETRF_AASEN.
  91 *> \endverbatim
  92 *>
  93 *> \param[in] LDA
  94 *> \verbatim
  95 *>          LDA is INTEGER
  96 *>          The leading dimension of the array A.  LDA >= max(1,N).
  97 *> \endverbatim
  98 *>
  99 *> \param[out] IPIV
 100 *> \verbatim
 101 *>          IPIV is INTEGER array, dimension (N)
 102 *>          On exit, it contains the details of the interchanges, i.e.,
 103 *>          the row and column k of A were interchanged with the
 104 *>          row and column IPIV(k).
 105 *> \endverbatim
 106 *>
 107 *> \param[in,out] B
 108 *> \verbatim
 109 *>          B is COMPLEX*16 array, dimension (LDB,NRHS)
 110 *>          On entry, the N-by-NRHS right hand side matrix B.
 111 *>          On exit, if INFO = 0, the N-by-NRHS solution matrix X.
 112 *> \endverbatim
 113 *>
 114 *> \param[in] LDB
 115 *> \verbatim
 116 *>          LDB is INTEGER
 117 *>          The leading dimension of the array B.  LDB >= max(1,N).
 118 *> \endverbatim
 119 *>
 120 *> \param[out] WORK
 121 *> \verbatim
 122 *>          WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
 123 *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
 124 *> \endverbatim
 125 *>
 126 *> \param[in] LWORK
 127 *> \verbatim
 128 *>          LWORK is INTEGER
 129 *>          The length of WORK.  LWORK >= 1, and for best performance
 130 *>          LWORK >= max(1,N*NB), where NB is the optimal blocksize for
 131 *>          ZHETRF.
 132 *>          for LWORK < N, TRS will be done with Level BLAS 2
 133 *>          for LWORK >= N, TRS will be done with Level BLAS 3
 134 *>
 135 *>          If LWORK = -1, then a workspace query is assumed; the routine
 136 *>          only calculates the optimal size of the WORK array, returns
 137 *>          this value as the first entry of the WORK array, and no error
 138 *>          message related to LWORK is issued by XERBLA.
 139 *> \endverbatim
 140 *>
 141 *> \param[out] INFO
 142 *> \verbatim
 143 *>          INFO is INTEGER
 144 *>          = 0: successful exit
 145 *>          < 0: if INFO = -i, the i-th argument had an illegal value
 146 *>          > 0: if INFO = i, D(i,i) is exactly zero.  The factorization
 147 *>               has been completed, but the block diagonal matrix D is
 148 *>               exactly singular, so the solution could not be computed.
 149 *> \endverbatim
 150 *
 151 *  Authors:
 152 *  ========
 153 *
 154 *> \author Univ. of Tennessee
 155 *> \author Univ. of California Berkeley
 156 *> \author Univ. of Colorado Denver
 157 *> \author NAG Ltd.
 158 *
 159 *> \date November 2016
 160 *
 161 *> \ingroup complex16HEsolve
 162 *
 163 *  @precisions fortran z -> c
 164 *
 165 *  =====================================================================
 166       SUBROUTINE ZHESV_AASEN( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, WORK,
 167      $                        LWORK, INFO )
 168 *
 169 *  -- LAPACK driver routine (version 3.4.0) --
 170 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 171 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 172 *     November 2016
 173 *
 174 *     .. Scalar Arguments ..
 175       CHARACTER          UPLO
 176       INTEGER            INFO, LDA, LDB, LWORK, N, NRHS
 177 *     ..
 178 *     .. Array Arguments ..
 179       INTEGER            IPIV( * )
 180       COMPLEX*16         A( LDA, * ), B( LDB, * ), WORK( * )
 181 *     ..
 182 *
 183 *  =====================================================================
 184 *
 185 *     .. Local Scalars ..
 186       LOGICAL            LQUERY
 187       INTEGER            LWKOPT, NB
 188 *     ..
 189 *     .. External Functions ..
 190       LOGICAL            LSAME
 191       INTEGER            ILAENV
 192       EXTERNAL           LSAME, ILAENV
 193 *     ..
 194 *     .. External Subroutines ..
 195       EXTERNAL           XERBLA, ZHETRF, ZHETRS, ZHETRS2
 196 *     ..
 197 *     .. Intrinsic Functions ..
 198       INTRINSIC          MAX
 199 *     ..
 200 *     .. Executable Statements ..
 201 *
 202 *     Test the input parameters.
 203 *
 204       INFO = 0
 205       LQUERY = ( LWORK.EQ.-1 )
 206       IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
 207          INFO = -1
 208       ELSE IF( N.LT.0 ) THEN
 209          INFO = -2
 210       ELSE IF( NRHS.LT.0 ) THEN
 211          INFO = -3
 212       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
 213          INFO = -5
 214       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
 215          INFO = -8
 216       ELSE IF( LWORK.LT.MAX(2*N, 3*N-2) .AND. .NOT.LQUERY ) THEN
 217          INFO = -10
 218       END IF
 219 *
 220       IF( INFO.EQ.0 ) THEN
 221          NB = ILAENV( 1, 'ZHETRF_AASEN', UPLO, N, -1, -1, -1 )
 222          LWKOPT = MAX( 3*N-2, (1+NB)*N )
 223          WORK( 1 ) = LWKOPT
 224       END IF
 225 *
 226       IF( INFO.NE.0 ) THEN
 227          CALL XERBLA( 'ZHESV_AASEN ', -INFO )
 228          RETURN
 229       ELSE IF( LQUERY ) THEN
 230          RETURN
 231       END IF
 232 *
 233 *     Compute the factorization A = U*T*U**H or A = L*T*L**H.
 234 *
 235       CALL ZHETRF_AASEN( UPLO, N, A, LDA, IPIV, WORK, LWORK, INFO )
 236       IF( INFO.EQ.0 ) THEN
 237 *
 238 *        Solve the system A*X = B, overwriting B with X.
 239 *
 240          CALL ZHETRS_AASEN( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, WORK,
 241      $                      LWORK, INFO )
 242 *
 243       END IF
 244 *
 245       WORK( 1 ) = LWKOPT
 246 *
 247       RETURN
 248 *
 249 *     End of ZHESV_AASEN
 250 *
 251       END