LAPACK 3.3.0

zpot01.f

Go to the documentation of this file.
00001       SUBROUTINE ZPOT01( UPLO, N, A, LDA, AFAC, LDAFAC, RWORK, RESID )
00002 *
00003 *  -- LAPACK test routine (version 3.1) --
00004 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
00005 *     November 2006
00006 *
00007 *     .. Scalar Arguments ..
00008       CHARACTER          UPLO
00009       INTEGER            LDA, LDAFAC, N
00010       DOUBLE PRECISION   RESID
00011 *     ..
00012 *     .. Array Arguments ..
00013       DOUBLE PRECISION   RWORK( * )
00014       COMPLEX*16         A( LDA, * ), AFAC( LDAFAC, * )
00015 *     ..
00016 *
00017 *  Purpose
00018 *  =======
00019 *
00020 *  ZPOT01 reconstructs a Hermitian positive definite matrix  A  from
00021 *  its L*L' or U'*U factorization and computes the residual
00022 *     norm( L*L' - A ) / ( N * norm(A) * EPS ) or
00023 *     norm( U'*U - A ) / ( N * norm(A) * EPS ),
00024 *  where EPS is the machine epsilon, L' is the conjugate transpose of L,
00025 *  and U' is the conjugate transpose of U.
00026 *
00027 *  Arguments
00028 *  ==========
00029 *
00030 *  UPLO    (input) CHARACTER*1
00031 *          Specifies whether the upper or lower triangular part of the
00032 *          Hermitian matrix A is stored:
00033 *          = 'U':  Upper triangular
00034 *          = 'L':  Lower triangular
00035 *
00036 *  N       (input) INTEGER
00037 *          The number of rows and columns of the matrix A.  N >= 0.
00038 *
00039 *  A       (input) COMPLEX*16 array, dimension (LDA,N)
00040 *          The original Hermitian matrix A.
00041 *
00042 *  LDA     (input) INTEGER
00043 *          The leading dimension of the array A.  LDA >= max(1,N)
00044 *
00045 *  AFAC    (input/output) COMPLEX*16 array, dimension (LDAFAC,N)
00046 *          On entry, the factor L or U from the L*L' or U'*U
00047 *          factorization of A.
00048 *          Overwritten with the reconstructed matrix, and then with the
00049 *          difference L*L' - A (or U'*U - A).
00050 *
00051 *  LDAFAC  (input) INTEGER
00052 *          The leading dimension of the array AFAC.  LDAFAC >= max(1,N).
00053 *
00054 *  RWORK   (workspace) DOUBLE PRECISION array, dimension (N)
00055 *
00056 *  RESID   (output) DOUBLE PRECISION
00057 *          If UPLO = 'L', norm(L*L' - A) / ( N * norm(A) * EPS )
00058 *          If UPLO = 'U', norm(U'*U - A) / ( N * norm(A) * EPS )
00059 *
00060 *  =====================================================================
00061 *
00062 *     .. Parameters ..
00063       DOUBLE PRECISION   ZERO, ONE
00064       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
00065 *     ..
00066 *     .. Local Scalars ..
00067       INTEGER            I, J, K
00068       DOUBLE PRECISION   ANORM, EPS, TR
00069       COMPLEX*16         TC
00070 *     ..
00071 *     .. External Functions ..
00072       LOGICAL            LSAME
00073       DOUBLE PRECISION   DLAMCH, ZLANHE
00074       COMPLEX*16         ZDOTC
00075       EXTERNAL           LSAME, DLAMCH, ZLANHE, ZDOTC
00076 *     ..
00077 *     .. External Subroutines ..
00078       EXTERNAL           ZHER, ZSCAL, ZTRMV
00079 *     ..
00080 *     .. Intrinsic Functions ..
00081       INTRINSIC          DBLE, DIMAG
00082 *     ..
00083 *     .. Executable Statements ..
00084 *
00085 *     Quick exit if N = 0.
00086 *
00087       IF( N.LE.0 ) THEN
00088          RESID = ZERO
00089          RETURN
00090       END IF
00091 *
00092 *     Exit with RESID = 1/EPS if ANORM = 0.
00093 *
00094       EPS = DLAMCH( 'Epsilon' )
00095       ANORM = ZLANHE( '1', UPLO, N, A, LDA, RWORK )
00096       IF( ANORM.LE.ZERO ) THEN
00097          RESID = ONE / EPS
00098          RETURN
00099       END IF
00100 *
00101 *     Check the imaginary parts of the diagonal elements and return with
00102 *     an error code if any are nonzero.
00103 *
00104       DO 10 J = 1, N
00105          IF( DIMAG( AFAC( J, J ) ).NE.ZERO ) THEN
00106             RESID = ONE / EPS
00107             RETURN
00108          END IF
00109    10 CONTINUE
00110 *
00111 *     Compute the product U'*U, overwriting U.
00112 *
00113       IF( LSAME( UPLO, 'U' ) ) THEN
00114          DO 20 K = N, 1, -1
00115 *
00116 *           Compute the (K,K) element of the result.
00117 *
00118             TR = ZDOTC( K, AFAC( 1, K ), 1, AFAC( 1, K ), 1 )
00119             AFAC( K, K ) = TR
00120 *
00121 *           Compute the rest of column K.
00122 *
00123             CALL ZTRMV( 'Upper', 'Conjugate', 'Non-unit', K-1, AFAC,
00124      $                  LDAFAC, AFAC( 1, K ), 1 )
00125 *
00126    20    CONTINUE
00127 *
00128 *     Compute the product L*L', overwriting L.
00129 *
00130       ELSE
00131          DO 30 K = N, 1, -1
00132 *
00133 *           Add a multiple of column K of the factor L to each of
00134 *           columns K+1 through N.
00135 *
00136             IF( K+1.LE.N )
00137      $         CALL ZHER( 'Lower', N-K, ONE, AFAC( K+1, K ), 1,
00138      $                    AFAC( K+1, K+1 ), LDAFAC )
00139 *
00140 *           Scale column K by the diagonal element.
00141 *
00142             TC = AFAC( K, K )
00143             CALL ZSCAL( N-K+1, TC, AFAC( K, K ), 1 )
00144 *
00145    30    CONTINUE
00146       END IF
00147 *
00148 *     Compute the difference  L*L' - A (or U'*U - A).
00149 *
00150       IF( LSAME( UPLO, 'U' ) ) THEN
00151          DO 50 J = 1, N
00152             DO 40 I = 1, J - 1
00153                AFAC( I, J ) = AFAC( I, J ) - A( I, J )
00154    40       CONTINUE
00155             AFAC( J, J ) = AFAC( J, J ) - DBLE( A( J, J ) )
00156    50    CONTINUE
00157       ELSE
00158          DO 70 J = 1, N
00159             AFAC( J, J ) = AFAC( J, J ) - DBLE( A( J, J ) )
00160             DO 60 I = J + 1, N
00161                AFAC( I, J ) = AFAC( I, J ) - A( I, J )
00162    60       CONTINUE
00163    70    CONTINUE
00164       END IF
00165 *
00166 *     Compute norm( L*U - A ) / ( N * norm(A) * EPS )
00167 *
00168       RESID = ZLANHE( '1', UPLO, N, AFAC, LDAFAC, RWORK )
00169 *
00170       RESID = ( ( RESID / DBLE( N ) ) / ANORM ) / EPS
00171 *
00172       RETURN
00173 *
00174 *     End of ZPOT01
00175 *
00176       END
 All Files Functions