LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
zhegst.f
Go to the documentation of this file.
00001 *> \brief \b ZHEGST
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download ZHEGST + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhegst.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhegst.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhegst.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE ZHEGST( ITYPE, UPLO, N, A, LDA, B, LDB, INFO )
00022 * 
00023 *       .. Scalar Arguments ..
00024 *       CHARACTER          UPLO
00025 *       INTEGER            INFO, ITYPE, LDA, LDB, N
00026 *       ..
00027 *       .. Array Arguments ..
00028 *       COMPLEX*16         A( LDA, * ), B( LDB, * )
00029 *       ..
00030 *  
00031 *
00032 *> \par Purpose:
00033 *  =============
00034 *>
00035 *> \verbatim
00036 *>
00037 *> ZHEGST reduces a complex Hermitian-definite generalized
00038 *> eigenproblem to standard form.
00039 *>
00040 *> If ITYPE = 1, the problem is A*x = lambda*B*x,
00041 *> and A is overwritten by inv(U**H)*A*inv(U) or inv(L)*A*inv(L**H)
00042 *>
00043 *> If ITYPE = 2 or 3, the problem is A*B*x = lambda*x or
00044 *> B*A*x = lambda*x, and A is overwritten by U*A*U**H or L**H*A*L.
00045 *>
00046 *> B must have been previously factorized as U**H*U or L*L**H by ZPOTRF.
00047 *> \endverbatim
00048 *
00049 *  Arguments:
00050 *  ==========
00051 *
00052 *> \param[in] ITYPE
00053 *> \verbatim
00054 *>          ITYPE is INTEGER
00055 *>          = 1: compute inv(U**H)*A*inv(U) or inv(L)*A*inv(L**H);
00056 *>          = 2 or 3: compute U*A*U**H or L**H*A*L.
00057 *> \endverbatim
00058 *>
00059 *> \param[in] UPLO
00060 *> \verbatim
00061 *>          UPLO is CHARACTER*1
00062 *>          = 'U':  Upper triangle of A is stored and B is factored as
00063 *>                  U**H*U;
00064 *>          = 'L':  Lower triangle of A is stored and B is factored as
00065 *>                  L*L**H.
00066 *> \endverbatim
00067 *>
00068 *> \param[in] N
00069 *> \verbatim
00070 *>          N is INTEGER
00071 *>          The order of the matrices A and B.  N >= 0.
00072 *> \endverbatim
00073 *>
00074 *> \param[in,out] A
00075 *> \verbatim
00076 *>          A is COMPLEX*16 array, dimension (LDA,N)
00077 *>          On entry, the Hermitian matrix A.  If UPLO = 'U', the leading
00078 *>          N-by-N upper triangular part of A contains the upper
00079 *>          triangular part of the matrix A, and the strictly lower
00080 *>          triangular part of A is not referenced.  If UPLO = 'L', the
00081 *>          leading N-by-N lower triangular part of A contains the lower
00082 *>          triangular part of the matrix A, and the strictly upper
00083 *>          triangular part of A is not referenced.
00084 *>
00085 *>          On exit, if INFO = 0, the transformed matrix, stored in the
00086 *>          same format as A.
00087 *> \endverbatim
00088 *>
00089 *> \param[in] LDA
00090 *> \verbatim
00091 *>          LDA is INTEGER
00092 *>          The leading dimension of the array A.  LDA >= max(1,N).
00093 *> \endverbatim
00094 *>
00095 *> \param[in] B
00096 *> \verbatim
00097 *>          B is COMPLEX*16 array, dimension (LDB,N)
00098 *>          The triangular factor from the Cholesky factorization of B,
00099 *>          as returned by ZPOTRF.
00100 *> \endverbatim
00101 *>
00102 *> \param[in] LDB
00103 *> \verbatim
00104 *>          LDB is INTEGER
00105 *>          The leading dimension of the array B.  LDB >= max(1,N).
00106 *> \endverbatim
00107 *>
00108 *> \param[out] INFO
00109 *> \verbatim
00110 *>          INFO is INTEGER
00111 *>          = 0:  successful exit
00112 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
00113 *> \endverbatim
00114 *
00115 *  Authors:
00116 *  ========
00117 *
00118 *> \author Univ. of Tennessee 
00119 *> \author Univ. of California Berkeley 
00120 *> \author Univ. of Colorado Denver 
00121 *> \author NAG Ltd. 
00122 *
00123 *> \date November 2011
00124 *
00125 *> \ingroup complex16HEcomputational
00126 *
00127 *  =====================================================================
00128       SUBROUTINE ZHEGST( ITYPE, UPLO, N, A, LDA, B, LDB, INFO )
00129 *
00130 *  -- LAPACK computational routine (version 3.4.0) --
00131 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00132 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00133 *     November 2011
00134 *
00135 *     .. Scalar Arguments ..
00136       CHARACTER          UPLO
00137       INTEGER            INFO, ITYPE, LDA, LDB, N
00138 *     ..
00139 *     .. Array Arguments ..
00140       COMPLEX*16         A( LDA, * ), B( LDB, * )
00141 *     ..
00142 *
00143 *  =====================================================================
00144 *
00145 *     .. Parameters ..
00146       DOUBLE PRECISION   ONE
00147       PARAMETER          ( ONE = 1.0D+0 )
00148       COMPLEX*16         CONE, HALF
00149       PARAMETER          ( CONE = ( 1.0D+0, 0.0D+0 ),
00150      $                   HALF = ( 0.5D+0, 0.0D+0 ) )
00151 *     ..
00152 *     .. Local Scalars ..
00153       LOGICAL            UPPER
00154       INTEGER            K, KB, NB
00155 *     ..
00156 *     .. External Subroutines ..
00157       EXTERNAL           XERBLA, ZHEGS2, ZHEMM, ZHER2K, ZTRMM, ZTRSM
00158 *     ..
00159 *     .. Intrinsic Functions ..
00160       INTRINSIC          MAX, MIN
00161 *     ..
00162 *     .. External Functions ..
00163       LOGICAL            LSAME
00164       INTEGER            ILAENV
00165       EXTERNAL           LSAME, ILAENV
00166 *     ..
00167 *     .. Executable Statements ..
00168 *
00169 *     Test the input parameters.
00170 *
00171       INFO = 0
00172       UPPER = LSAME( UPLO, 'U' )
00173       IF( ITYPE.LT.1 .OR. ITYPE.GT.3 ) THEN
00174          INFO = -1
00175       ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
00176          INFO = -2
00177       ELSE IF( N.LT.0 ) THEN
00178          INFO = -3
00179       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
00180          INFO = -5
00181       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
00182          INFO = -7
00183       END IF
00184       IF( INFO.NE.0 ) THEN
00185          CALL XERBLA( 'ZHEGST', -INFO )
00186          RETURN
00187       END IF
00188 *
00189 *     Quick return if possible
00190 *
00191       IF( N.EQ.0 )
00192      $   RETURN
00193 *
00194 *     Determine the block size for this environment.
00195 *
00196       NB = ILAENV( 1, 'ZHEGST', UPLO, N, -1, -1, -1 )
00197 *
00198       IF( NB.LE.1 .OR. NB.GE.N ) THEN
00199 *
00200 *        Use unblocked code
00201 *
00202          CALL ZHEGS2( ITYPE, UPLO, N, A, LDA, B, LDB, INFO )
00203       ELSE
00204 *
00205 *        Use blocked code
00206 *
00207          IF( ITYPE.EQ.1 ) THEN
00208             IF( UPPER ) THEN
00209 *
00210 *              Compute inv(U**H)*A*inv(U)
00211 *
00212                DO 10 K = 1, N, NB
00213                   KB = MIN( N-K+1, NB )
00214 *
00215 *                 Update the upper triangle of A(k:n,k:n)
00216 *
00217                   CALL ZHEGS2( ITYPE, UPLO, KB, A( K, K ), LDA,
00218      $                         B( K, K ), LDB, INFO )
00219                   IF( K+KB.LE.N ) THEN
00220                      CALL ZTRSM( 'Left', UPLO, 'Conjugate transpose',
00221      $                           'Non-unit', KB, N-K-KB+1, CONE,
00222      $                           B( K, K ), LDB, A( K, K+KB ), LDA )
00223                      CALL ZHEMM( 'Left', UPLO, KB, N-K-KB+1, -HALF,
00224      $                           A( K, K ), LDA, B( K, K+KB ), LDB,
00225      $                           CONE, A( K, K+KB ), LDA )
00226                      CALL ZHER2K( UPLO, 'Conjugate transpose', N-K-KB+1,
00227      $                            KB, -CONE, A( K, K+KB ), LDA,
00228      $                            B( K, K+KB ), LDB, ONE,
00229      $                            A( K+KB, K+KB ), LDA )
00230                      CALL ZHEMM( 'Left', UPLO, KB, N-K-KB+1, -HALF,
00231      $                           A( K, K ), LDA, B( K, K+KB ), LDB,
00232      $                           CONE, A( K, K+KB ), LDA )
00233                      CALL ZTRSM( 'Right', UPLO, 'No transpose',
00234      $                           'Non-unit', KB, N-K-KB+1, CONE,
00235      $                           B( K+KB, K+KB ), LDB, A( K, K+KB ),
00236      $                           LDA )
00237                   END IF
00238    10          CONTINUE
00239             ELSE
00240 *
00241 *              Compute inv(L)*A*inv(L**H)
00242 *
00243                DO 20 K = 1, N, NB
00244                   KB = MIN( N-K+1, NB )
00245 *
00246 *                 Update the lower triangle of A(k:n,k:n)
00247 *
00248                   CALL ZHEGS2( ITYPE, UPLO, KB, A( K, K ), LDA,
00249      $                         B( K, K ), LDB, INFO )
00250                   IF( K+KB.LE.N ) THEN
00251                      CALL ZTRSM( 'Right', UPLO, 'Conjugate transpose',
00252      $                           'Non-unit', N-K-KB+1, KB, CONE,
00253      $                           B( K, K ), LDB, A( K+KB, K ), LDA )
00254                      CALL ZHEMM( 'Right', UPLO, N-K-KB+1, KB, -HALF,
00255      $                           A( K, K ), LDA, B( K+KB, K ), LDB,
00256      $                           CONE, A( K+KB, K ), LDA )
00257                      CALL ZHER2K( UPLO, 'No transpose', N-K-KB+1, KB,
00258      $                            -CONE, A( K+KB, K ), LDA,
00259      $                            B( K+KB, K ), LDB, ONE,
00260      $                            A( K+KB, K+KB ), LDA )
00261                      CALL ZHEMM( 'Right', UPLO, N-K-KB+1, KB, -HALF,
00262      $                           A( K, K ), LDA, B( K+KB, K ), LDB,
00263      $                           CONE, A( K+KB, K ), LDA )
00264                      CALL ZTRSM( 'Left', UPLO, 'No transpose',
00265      $                           'Non-unit', N-K-KB+1, KB, CONE,
00266      $                           B( K+KB, K+KB ), LDB, A( K+KB, K ),
00267      $                           LDA )
00268                   END IF
00269    20          CONTINUE
00270             END IF
00271          ELSE
00272             IF( UPPER ) THEN
00273 *
00274 *              Compute U*A*U**H
00275 *
00276                DO 30 K = 1, N, NB
00277                   KB = MIN( N-K+1, NB )
00278 *
00279 *                 Update the upper triangle of A(1:k+kb-1,1:k+kb-1)
00280 *
00281                   CALL ZTRMM( 'Left', UPLO, 'No transpose', 'Non-unit',
00282      $                        K-1, KB, CONE, B, LDB, A( 1, K ), LDA )
00283                   CALL ZHEMM( 'Right', UPLO, K-1, KB, HALF, A( K, K ),
00284      $                        LDA, B( 1, K ), LDB, CONE, A( 1, K ),
00285      $                        LDA )
00286                   CALL ZHER2K( UPLO, 'No transpose', K-1, KB, CONE,
00287      $                         A( 1, K ), LDA, B( 1, K ), LDB, ONE, A,
00288      $                         LDA )
00289                   CALL ZHEMM( 'Right', UPLO, K-1, KB, HALF, A( K, K ),
00290      $                        LDA, B( 1, K ), LDB, CONE, A( 1, K ),
00291      $                        LDA )
00292                   CALL ZTRMM( 'Right', UPLO, 'Conjugate transpose',
00293      $                        'Non-unit', K-1, KB, CONE, B( K, K ), LDB,
00294      $                        A( 1, K ), LDA )
00295                   CALL ZHEGS2( ITYPE, UPLO, KB, A( K, K ), LDA,
00296      $                         B( K, K ), LDB, INFO )
00297    30          CONTINUE
00298             ELSE
00299 *
00300 *              Compute L**H*A*L
00301 *
00302                DO 40 K = 1, N, NB
00303                   KB = MIN( N-K+1, NB )
00304 *
00305 *                 Update the lower triangle of A(1:k+kb-1,1:k+kb-1)
00306 *
00307                   CALL ZTRMM( 'Right', UPLO, 'No transpose', 'Non-unit',
00308      $                        KB, K-1, CONE, B, LDB, A( K, 1 ), LDA )
00309                   CALL ZHEMM( 'Left', UPLO, KB, K-1, HALF, A( K, K ),
00310      $                        LDA, B( K, 1 ), LDB, CONE, A( K, 1 ),
00311      $                        LDA )
00312                   CALL ZHER2K( UPLO, 'Conjugate transpose', K-1, KB,
00313      $                         CONE, A( K, 1 ), LDA, B( K, 1 ), LDB,
00314      $                         ONE, A, LDA )
00315                   CALL ZHEMM( 'Left', UPLO, KB, K-1, HALF, A( K, K ),
00316      $                        LDA, B( K, 1 ), LDB, CONE, A( K, 1 ),
00317      $                        LDA )
00318                   CALL ZTRMM( 'Left', UPLO, 'Conjugate transpose',
00319      $                        'Non-unit', KB, K-1, CONE, B( K, K ), LDB,
00320      $                        A( K, 1 ), LDA )
00321                   CALL ZHEGS2( ITYPE, UPLO, KB, A( K, K ), LDA,
00322      $                         B( K, K ), LDB, INFO )
00323    40          CONTINUE
00324             END IF
00325          END IF
00326       END IF
00327       RETURN
00328 *
00329 *     End of ZHEGST
00330 *
00331       END
 All Files Functions