LAPACK: SRC/chegs2.f Source File

Go to the documentation of this file.
00001       SUBROUTINE CHEGS2( ITYPE, UPLO, N, A, LDA, B, LDB, INFO )
00002 *
00003 *  -- LAPACK routine (version 3.2) --
00004 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00005 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00006 *     November 2006
00007 *
00008 *     .. Scalar Arguments ..
00009       CHARACTER          UPLO
00010       INTEGER            INFO, ITYPE, LDA, LDB, N
00011 *     ..
00012 *     .. Array Arguments ..
00013       COMPLEX            A( LDA, * ), B( LDB, * )
00014 *     ..
00015 *
00016 *  Purpose
00017 *  =======
00018 *
00019 *  CHEGS2 reduces a complex Hermitian-definite generalized
00020 *  eigenproblem to standard form.
00021 *
00022 *  If ITYPE = 1, the problem is A*x = lambda*B*x,
00023 *  and A is overwritten by inv(U')*A*inv(U) or inv(L)*A*inv(L')
00024 *
00025 *  If ITYPE = 2 or 3, the problem is A*B*x = lambda*x or
00026 *  B*A*x = lambda*x, and A is overwritten by U*A*U` or L'*A*L.
00027 *
00028 *  B must have been previously factorized as U'*U or L*L' by CPOTRF.
00029 *
00030 *  Arguments
00031 *  =========
00032 *
00033 *  ITYPE   (input) INTEGER
00034 *          = 1: compute inv(U')*A*inv(U) or inv(L)*A*inv(L');
00035 *          = 2 or 3: compute U*A*U' or L'*A*L.
00036 *
00037 *  UPLO    (input) CHARACTER*1
00038 *          Specifies whether the upper or lower triangular part of the
00039 *          Hermitian matrix A is stored, and how B has been factorized.
00040 *          = 'U':  Upper triangular
00041 *          = 'L':  Lower triangular
00042 *
00043 *  N       (input) INTEGER
00044 *          The order of the matrices A and B.  N >= 0.
00045 *
00046 *  A       (input/output) COMPLEX array, dimension (LDA,N)
00047 *          On entry, the Hermitian matrix A.  If UPLO = 'U', the leading
00048 *          n by n upper triangular part of A contains the upper
00049 *          triangular part of the matrix A, and the strictly lower
00050 *          triangular part of A is not referenced.  If UPLO = 'L', the
00051 *          leading n by n lower triangular part of A contains the lower
00052 *          triangular part of the matrix A, and the strictly upper
00053 *          triangular part of A is not referenced.
00054 *
00055 *          On exit, if INFO = 0, the transformed matrix, stored in the
00056 *          same format as A.
00057 *
00058 *  LDA     (input) INTEGER
00059 *          The leading dimension of the array A.  LDA >= max(1,N).
00060 *
00061 *  B       (input) COMPLEX array, dimension (LDB,N)
00062 *          The triangular factor from the Cholesky factorization of B,
00063 *          as returned by CPOTRF.
00064 *
00065 *  LDB     (input) INTEGER
00066 *          The leading dimension of the array B.  LDB >= max(1,N).
00067 *
00068 *  INFO    (output) INTEGER
00069 *          = 0:  successful exit.
00070 *          < 0:  if INFO = -i, the i-th argument had an illegal value.
00071 *
00072 *  =====================================================================
00073 *
00074 *     .. Parameters ..
00075       REAL               ONE, HALF
00076       PARAMETER          ( ONE = 1.0E+0, HALF = 0.5E+0 )
00077       COMPLEX            CONE
00078       PARAMETER          ( CONE = ( 1.0E+0, 0.0E+0 ) )
00079 *     ..
00080 *     .. Local Scalars ..
00081       LOGICAL            UPPER
00082       INTEGER            K
00083       REAL               AKK, BKK
00084       COMPLEX            CT
00085 *     ..
00086 *     .. External Subroutines ..
00087       EXTERNAL           CAXPY, CHER2, CLACGV, CSSCAL, CTRMV, CTRSV,
00088      $                   XERBLA
00089 *     ..
00090 *     .. Intrinsic Functions ..
00091       INTRINSIC          MAX
00092 *     ..
00093 *     .. External Functions ..
00094       LOGICAL            LSAME
00095       EXTERNAL           LSAME
00096 *     ..
00097 *     .. Executable Statements ..
00098 *
00099 *     Test the input parameters.
00100 *
00101       INFO = 0
00102       UPPER = LSAME( UPLO, 'U' )
00103       IF( ITYPE.LT.1 .OR. ITYPE.GT.3 ) THEN
00104          INFO = -1
00105       ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
00106          INFO = -2
00107       ELSE IF( N.LT.0 ) THEN
00108          INFO = -3
00109       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
00110          INFO = -5
00111       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
00112          INFO = -7
00113       END IF
00114       IF( INFO.NE.0 ) THEN
00115          CALL XERBLA( 'CHEGS2', -INFO )
00116          RETURN
00117       END IF
00118 *
00119       IF( ITYPE.EQ.1 ) THEN
00120          IF( UPPER ) THEN
00121 *
00122 *           Compute inv(U')*A*inv(U)
00123 *
00124             DO 10 K = 1, N
00125 *
00126 *              Update the upper triangle of A(k:n,k:n)
00127 *
00128                AKK = A( K, K )
00129                BKK = B( K, K )
00130                AKK = AKK / BKK**2
00131                A( K, K ) = AKK
00132                IF( K.LT.N ) THEN
00133                   CALL CSSCAL( N-K, ONE / BKK, A( K, K+1 ), LDA )
00134                   CT = -HALF*AKK
00135                   CALL CLACGV( N-K, A( K, K+1 ), LDA )
00136                   CALL CLACGV( N-K, B( K, K+1 ), LDB )
00137                   CALL CAXPY( N-K, CT, B( K, K+1 ), LDB, A( K, K+1 ),
00138      $                        LDA )
00139                   CALL CHER2( UPLO, N-K, -CONE, A( K, K+1 ), LDA,
00140      $                        B( K, K+1 ), LDB, A( K+1, K+1 ), LDA )
00141                   CALL CAXPY( N-K, CT, B( K, K+1 ), LDB, A( K, K+1 ),
00142      $                        LDA )
00143                   CALL CLACGV( N-K, B( K, K+1 ), LDB )
00144                   CALL CTRSV( UPLO, 'Conjugate transpose', 'Non-unit',
00145      $                        N-K, B( K+1, K+1 ), LDB, A( K, K+1 ),
00146      $                        LDA )
00147                   CALL CLACGV( N-K, A( K, K+1 ), LDA )
00148                END IF
00149    10       CONTINUE
00150          ELSE
00151 *
00152 *           Compute inv(L)*A*inv(L')
00153 *
00154             DO 20 K = 1, N
00155 *
00156 *              Update the lower triangle of A(k:n,k:n)
00157 *
00158                AKK = A( K, K )
00159                BKK = B( K, K )
00160                AKK = AKK / BKK**2
00161                A( K, K ) = AKK
00162                IF( K.LT.N ) THEN
00163                   CALL CSSCAL( N-K, ONE / BKK, A( K+1, K ), 1 )
00164                   CT = -HALF*AKK
00165                   CALL CAXPY( N-K, CT, B( K+1, K ), 1, A( K+1, K ), 1 )
00166                   CALL CHER2( UPLO, N-K, -CONE, A( K+1, K ), 1,
00167      $                        B( K+1, K ), 1, A( K+1, K+1 ), LDA )
00168                   CALL CAXPY( N-K, CT, B( K+1, K ), 1, A( K+1, K ), 1 )
00169                   CALL CTRSV( UPLO, 'No transpose', 'Non-unit', N-K,
00170      $                        B( K+1, K+1 ), LDB, A( K+1, K ), 1 )
00171                END IF
00172    20       CONTINUE
00173          END IF
00174       ELSE
00175          IF( UPPER ) THEN
00176 *
00177 *           Compute U*A*U'
00178 *
00179             DO 30 K = 1, N
00180 *
00181 *              Update the upper triangle of A(1:k,1:k)
00182 *
00183                AKK = A( K, K )
00184                BKK = B( K, K )
00185                CALL CTRMV( UPLO, 'No transpose', 'Non-unit', K-1, B,
00186      $                     LDB, A( 1, K ), 1 )
00187                CT = HALF*AKK
00188                CALL CAXPY( K-1, CT, B( 1, K ), 1, A( 1, K ), 1 )
00189                CALL CHER2( UPLO, K-1, CONE, A( 1, K ), 1, B( 1, K ), 1,
00190      $                     A, LDA )
00191                CALL CAXPY( K-1, CT, B( 1, K ), 1, A( 1, K ), 1 )
00192                CALL CSSCAL( K-1, BKK, A( 1, K ), 1 )
00193                A( K, K ) = AKK*BKK**2
00194    30       CONTINUE
00195          ELSE
00196 *
00197 *           Compute L'*A*L
00198 *
00199             DO 40 K = 1, N
00200 *
00201 *              Update the lower triangle of A(1:k,1:k)
00202 *
00203                AKK = A( K, K )
00204                BKK = B( K, K )
00205                CALL CLACGV( K-1, A( K, 1 ), LDA )
00206                CALL CTRMV( UPLO, 'Conjugate transpose', 'Non-unit', K-1,
00207      $                     B, LDB, A( K, 1 ), LDA )
00208                CT = HALF*AKK
00209                CALL CLACGV( K-1, B( K, 1 ), LDB )
00210                CALL CAXPY( K-1, CT, B( K, 1 ), LDB, A( K, 1 ), LDA )
00211                CALL CHER2( UPLO, K-1, CONE, A( K, 1 ), LDA, B( K, 1 ),
00212      $                     LDB, A, LDA )
00213                CALL CAXPY( K-1, CT, B( K, 1 ), LDB, A( K, 1 ), LDA )
00214                CALL CLACGV( K-1, B( K, 1 ), LDB )
00215                CALL CSSCAL( K-1, BKK, A( K, 1 ), LDA )
00216                CALL CLACGV( K-1, A( K, 1 ), LDA )
00217                A( K, K ) = AKK*BKK**2
00218    40       CONTINUE
00219          END IF
00220       END IF
00221       RETURN
00222 *
00223 *     End of CHEGS2
00224 *
00225       END