LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
chpt01.f
Go to the documentation of this file.
00001 *> \brief \b CHPT01
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *  Definition:
00009 *  ===========
00010 *
00011 *       SUBROUTINE CHPT01( UPLO, N, A, AFAC, IPIV, C, LDC, RWORK, RESID )
00012 * 
00013 *       .. Scalar Arguments ..
00014 *       CHARACTER          UPLO
00015 *       INTEGER            LDC, N
00016 *       REAL               RESID
00017 *       ..
00018 *       .. Array Arguments ..
00019 *       INTEGER            IPIV( * )
00020 *       REAL               RWORK( * )
00021 *       COMPLEX            A( * ), AFAC( * ), C( LDC, * )
00022 *       ..
00023 *  
00024 *
00025 *> \par Purpose:
00026 *  =============
00027 *>
00028 *> \verbatim
00029 *>
00030 *> CHPT01 reconstructs a Hermitian indefinite packed matrix A from its
00031 *> block L*D*L' or U*D*U' factorization and computes the residual
00032 *>    norm( C - A ) / ( N * norm(A) * EPS ),
00033 *> where C is the reconstructed matrix, EPS is the machine epsilon,
00034 *> L' is the conjugate transpose of L, and U' is the conjugate transpose
00035 *> of U.
00036 *> \endverbatim
00037 *
00038 *  Arguments:
00039 *  ==========
00040 *
00041 *> \param[in] UPLO
00042 *> \verbatim
00043 *>          UPLO is CHARACTER*1
00044 *>          Specifies whether the upper or lower triangular part of the
00045 *>          Hermitian matrix A is stored:
00046 *>          = 'U':  Upper triangular
00047 *>          = 'L':  Lower triangular
00048 *> \endverbatim
00049 *>
00050 *> \param[in] N
00051 *> \verbatim
00052 *>          N is INTEGER
00053 *>          The number of rows and columns of the matrix A.  N >= 0.
00054 *> \endverbatim
00055 *>
00056 *> \param[in] A
00057 *> \verbatim
00058 *>          A is COMPLEX array, dimension (N*(N+1)/2)
00059 *>          The original Hermitian matrix A, stored as a packed
00060 *>          triangular matrix.
00061 *> \endverbatim
00062 *>
00063 *> \param[in] AFAC
00064 *> \verbatim
00065 *>          AFAC is COMPLEX array, dimension (N*(N+1)/2)
00066 *>          The factored form of the matrix A, stored as a packed
00067 *>          triangular matrix.  AFAC contains the block diagonal matrix D
00068 *>          and the multipliers used to obtain the factor L or U from the
00069 *>          block L*D*L' or U*D*U' factorization as computed by CHPTRF.
00070 *> \endverbatim
00071 *>
00072 *> \param[in] IPIV
00073 *> \verbatim
00074 *>          IPIV is INTEGER array, dimension (N)
00075 *>          The pivot indices from CHPTRF.
00076 *> \endverbatim
00077 *>
00078 *> \param[out] C
00079 *> \verbatim
00080 *>          C is COMPLEX array, dimension (LDC,N)
00081 *> \endverbatim
00082 *>
00083 *> \param[in] LDC
00084 *> \verbatim
00085 *>          LDC is INTEGER
00086 *>          The leading dimension of the array C.  LDC >= max(1,N).
00087 *> \endverbatim
00088 *>
00089 *> \param[out] RWORK
00090 *> \verbatim
00091 *>          RWORK is REAL array, dimension (N)
00092 *> \endverbatim
00093 *>
00094 *> \param[out] RESID
00095 *> \verbatim
00096 *>          RESID is REAL
00097 *>          If UPLO = 'L', norm(L*D*L' - A) / ( N * norm(A) * EPS )
00098 *>          If UPLO = 'U', norm(U*D*U' - A) / ( N * norm(A) * EPS )
00099 *> \endverbatim
00100 *
00101 *  Authors:
00102 *  ========
00103 *
00104 *> \author Univ. of Tennessee 
00105 *> \author Univ. of California Berkeley 
00106 *> \author Univ. of Colorado Denver 
00107 *> \author NAG Ltd. 
00108 *
00109 *> \date November 2011
00110 *
00111 *> \ingroup complex_lin
00112 *
00113 *  =====================================================================
00114       SUBROUTINE CHPT01( UPLO, N, A, AFAC, IPIV, C, LDC, RWORK, RESID )
00115 *
00116 *  -- LAPACK test routine (version 3.4.0) --
00117 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00118 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00119 *     November 2011
00120 *
00121 *     .. Scalar Arguments ..
00122       CHARACTER          UPLO
00123       INTEGER            LDC, N
00124       REAL               RESID
00125 *     ..
00126 *     .. Array Arguments ..
00127       INTEGER            IPIV( * )
00128       REAL               RWORK( * )
00129       COMPLEX            A( * ), AFAC( * ), C( LDC, * )
00130 *     ..
00131 *
00132 *  =====================================================================
00133 *
00134 *     .. Parameters ..
00135       REAL               ZERO, ONE
00136       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
00137       COMPLEX            CZERO, CONE
00138       PARAMETER          ( CZERO = ( 0.0E+0, 0.0E+0 ),
00139      $                   CONE = ( 1.0E+0, 0.0E+0 ) )
00140 *     ..
00141 *     .. Local Scalars ..
00142       INTEGER            I, INFO, J, JC
00143       REAL               ANORM, EPS
00144 *     ..
00145 *     .. External Functions ..
00146       LOGICAL            LSAME
00147       REAL               CLANHE, CLANHP, SLAMCH
00148       EXTERNAL           LSAME, CLANHE, CLANHP, SLAMCH
00149 *     ..
00150 *     .. External Subroutines ..
00151       EXTERNAL           CLAVHP, CLASET
00152 *     ..
00153 *     .. Intrinsic Functions ..
00154       INTRINSIC          AIMAG, REAL
00155 *     ..
00156 *     .. Executable Statements ..
00157 *
00158 *     Quick exit if N = 0.
00159 *
00160       IF( N.LE.0 ) THEN
00161          RESID = ZERO
00162          RETURN
00163       END IF
00164 *
00165 *     Determine EPS and the norm of A.
00166 *
00167       EPS = SLAMCH( 'Epsilon' )
00168       ANORM = CLANHP( '1', UPLO, N, A, RWORK )
00169 *
00170 *     Check the imaginary parts of the diagonal elements and return with
00171 *     an error code if any are nonzero.
00172 *
00173       JC = 1
00174       IF( LSAME( UPLO, 'U' ) ) THEN
00175          DO 10 J = 1, N
00176             IF( AIMAG( AFAC( JC ) ).NE.ZERO ) THEN
00177                RESID = ONE / EPS
00178                RETURN
00179             END IF
00180             JC = JC + J + 1
00181    10    CONTINUE
00182       ELSE
00183          DO 20 J = 1, N
00184             IF( AIMAG( AFAC( JC ) ).NE.ZERO ) THEN
00185                RESID = ONE / EPS
00186                RETURN
00187             END IF
00188             JC = JC + N - J + 1
00189    20    CONTINUE
00190       END IF
00191 *
00192 *     Initialize C to the identity matrix.
00193 *
00194       CALL CLASET( 'Full', N, N, CZERO, CONE, C, LDC )
00195 *
00196 *     Call CLAVHP to form the product D * U' (or D * L' ).
00197 *
00198       CALL CLAVHP( UPLO, 'Conjugate', 'Non-unit', N, N, AFAC, IPIV, C,
00199      $             LDC, INFO )
00200 *
00201 *     Call CLAVHP again to multiply by U ( or L ).
00202 *
00203       CALL CLAVHP( UPLO, 'No transpose', 'Unit', N, N, AFAC, IPIV, C,
00204      $             LDC, INFO )
00205 *
00206 *     Compute the difference  C - A .
00207 *
00208       IF( LSAME( UPLO, 'U' ) ) THEN
00209          JC = 0
00210          DO 40 J = 1, N
00211             DO 30 I = 1, J - 1
00212                C( I, J ) = C( I, J ) - A( JC+I )
00213    30       CONTINUE
00214             C( J, J ) = C( J, J ) - REAL( A( JC+J ) )
00215             JC = JC + J
00216    40    CONTINUE
00217       ELSE
00218          JC = 1
00219          DO 60 J = 1, N
00220             C( J, J ) = C( J, J ) - REAL( A( JC ) )
00221             DO 50 I = J + 1, N
00222                C( I, J ) = C( I, J ) - A( JC+I-J )
00223    50       CONTINUE
00224             JC = JC + N - J + 1
00225    60    CONTINUE
00226       END IF
00227 *
00228 *     Compute norm( C - A ) / ( N * norm(A) * EPS )
00229 *
00230       RESID = CLANHE( '1', UPLO, N, C, LDC, RWORK )
00231 *
00232       IF( ANORM.LE.ZERO ) THEN
00233          IF( RESID.NE.ZERO )
00234      $      RESID = ONE / EPS
00235       ELSE
00236          RESID = ( ( RESID / REAL( N ) ) / ANORM ) / EPS
00237       END IF
00238 *
00239       RETURN
00240 *
00241 *     End of CHPT01
00242 *
00243       END
 All Files Functions