LAPACK 3.3.1 Linear Algebra PACKage

# zppt01.f

Go to the documentation of this file.
```00001       SUBROUTINE ZPPT01( UPLO, N, A, AFAC, 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            N
00010       DOUBLE PRECISION   RESID
00011 *     ..
00012 *     .. Array Arguments ..
00013       DOUBLE PRECISION   RWORK( * )
00014       COMPLEX*16         A( * ), AFAC( * )
00015 *     ..
00016 *
00017 *  Purpose
00018 *  =======
00019 *
00020 *  ZPPT01 reconstructs a Hermitian positive definite packed matrix A
00021 *  from 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
00025 *  L, 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 (N*(N+1)/2)
00040 *          The original Hermitian matrix A, stored as a packed
00041 *          triangular matrix.
00042 *
00043 *  AFAC    (input/output) COMPLEX*16 array, dimension (N*(N+1)/2)
00044 *          On entry, the factor L or U from the L*L' or U'*U
00045 *          factorization of A, stored as a packed triangular matrix.
00046 *          Overwritten with the reconstructed matrix, and then with the
00047 *          difference L*L' - A (or U'*U - A).
00048 *
00049 *  RWORK   (workspace) DOUBLE PRECISION array, dimension (N)
00050 *
00051 *  RESID   (output) DOUBLE PRECISION
00052 *          If UPLO = 'L', norm(L*L' - A) / ( N * norm(A) * EPS )
00053 *          If UPLO = 'U', norm(U'*U - A) / ( N * norm(A) * EPS )
00054 *
00055 *  =====================================================================
00056 *
00057 *     .. Parameters ..
00058       DOUBLE PRECISION   ZERO, ONE
00059       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
00060 *     ..
00061 *     .. Local Scalars ..
00062       INTEGER            I, K, KC
00063       DOUBLE PRECISION   ANORM, EPS, TR
00064       COMPLEX*16         TC
00065 *     ..
00066 *     .. External Functions ..
00067       LOGICAL            LSAME
00068       DOUBLE PRECISION   DLAMCH, ZLANHP
00069       COMPLEX*16         ZDOTC
00070       EXTERNAL           LSAME, DLAMCH, ZLANHP, ZDOTC
00071 *     ..
00072 *     .. External Subroutines ..
00073       EXTERNAL           ZHPR, ZSCAL, ZTPMV
00074 *     ..
00075 *     .. Intrinsic Functions ..
00076       INTRINSIC          DBLE, DIMAG
00077 *     ..
00078 *     .. Executable Statements ..
00079 *
00080 *     Quick exit if N = 0
00081 *
00082       IF( N.LE.0 ) THEN
00083          RESID = ZERO
00084          RETURN
00085       END IF
00086 *
00087 *     Exit with RESID = 1/EPS if ANORM = 0.
00088 *
00089       EPS = DLAMCH( 'Epsilon' )
00090       ANORM = ZLANHP( '1', UPLO, N, A, RWORK )
00091       IF( ANORM.LE.ZERO ) THEN
00092          RESID = ONE / EPS
00093          RETURN
00094       END IF
00095 *
00096 *     Check the imaginary parts of the diagonal elements and return with
00097 *     an error code if any are nonzero.
00098 *
00099       KC = 1
00100       IF( LSAME( UPLO, 'U' ) ) THEN
00101          DO 10 K = 1, N
00102             IF( DIMAG( AFAC( KC ) ).NE.ZERO ) THEN
00103                RESID = ONE / EPS
00104                RETURN
00105             END IF
00106             KC = KC + K + 1
00107    10    CONTINUE
00108       ELSE
00109          DO 20 K = 1, N
00110             IF( DIMAG( AFAC( KC ) ).NE.ZERO ) THEN
00111                RESID = ONE / EPS
00112                RETURN
00113             END IF
00114             KC = KC + N - K + 1
00115    20    CONTINUE
00116       END IF
00117 *
00118 *     Compute the product U'*U, overwriting U.
00119 *
00120       IF( LSAME( UPLO, 'U' ) ) THEN
00121          KC = ( N*( N-1 ) ) / 2 + 1
00122          DO 30 K = N, 1, -1
00123 *
00124 *           Compute the (K,K) element of the result.
00125 *
00126             TR = ZDOTC( K, AFAC( KC ), 1, AFAC( KC ), 1 )
00127             AFAC( KC+K-1 ) = TR
00128 *
00129 *           Compute the rest of column K.
00130 *
00131             IF( K.GT.1 ) THEN
00132                CALL ZTPMV( 'Upper', 'Conjugate', 'Non-unit', K-1, AFAC,
00133      \$                     AFAC( KC ), 1 )
00134                KC = KC - ( K-1 )
00135             END IF
00136    30    CONTINUE
00137 *
00138 *        Compute the difference  L*L' - A
00139 *
00140          KC = 1
00141          DO 50 K = 1, N
00142             DO 40 I = 1, K - 1
00143                AFAC( KC+I-1 ) = AFAC( KC+I-1 ) - A( KC+I-1 )
00144    40       CONTINUE
00145             AFAC( KC+K-1 ) = AFAC( KC+K-1 ) - DBLE( A( KC+K-1 ) )
00146             KC = KC + K
00147    50    CONTINUE
00148 *
00149 *     Compute the product L*L', overwriting L.
00150 *
00151       ELSE
00152          KC = ( N*( N+1 ) ) / 2
00153          DO 60 K = N, 1, -1
00154 *
00155 *           Add a multiple of column K of the factor L to each of
00156 *           columns K+1 through N.
00157 *
00158             IF( K.LT.N )
00159      \$         CALL ZHPR( 'Lower', N-K, ONE, AFAC( KC+1 ), 1,
00160      \$                    AFAC( KC+N-K+1 ) )
00161 *
00162 *           Scale column K by the diagonal element.
00163 *
00164             TC = AFAC( KC )
00165             CALL ZSCAL( N-K+1, TC, AFAC( KC ), 1 )
00166 *
00167             KC = KC - ( N-K+2 )
00168    60    CONTINUE
00169 *
00170 *        Compute the difference  U'*U - A
00171 *
00172          KC = 1
00173          DO 80 K = 1, N
00174             AFAC( KC ) = AFAC( KC ) - DBLE( A( KC ) )
00175             DO 70 I = K + 1, N
00176                AFAC( KC+I-K ) = AFAC( KC+I-K ) - A( KC+I-K )
00177    70       CONTINUE
00178             KC = KC + N - K + 1
00179    80    CONTINUE
00180       END IF
00181 *
00182 *     Compute norm( L*U - A ) / ( N * norm(A) * EPS )
00183 *
00184       RESID = ZLANHP( '1', UPLO, N, AFAC, RWORK )
00185 *
00186       RESID = ( ( RESID / DBLE( N ) ) / ANORM ) / EPS
00187 *
00188       RETURN
00189 *
00190 *     End of ZPPT01
00191 *
00192       END
```