LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
zppt01.f
Go to the documentation of this file.
00001 *> \brief \b ZPPT01
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 ZPPT01( UPLO, N, A, AFAC, RWORK, RESID )
00012 * 
00013 *       .. Scalar Arguments ..
00014 *       CHARACTER          UPLO
00015 *       INTEGER            N
00016 *       DOUBLE PRECISION   RESID
00017 *       ..
00018 *       .. Array Arguments ..
00019 *       DOUBLE PRECISION   RWORK( * )
00020 *       COMPLEX*16         A( * ), AFAC( * )
00021 *       ..
00022 *  
00023 *
00024 *> \par Purpose:
00025 *  =============
00026 *>
00027 *> \verbatim
00028 *>
00029 *> ZPPT01 reconstructs a Hermitian positive definite packed matrix A
00030 *> from its L*L' or U'*U factorization and computes the residual
00031 *>    norm( L*L' - A ) / ( N * norm(A) * EPS ) or
00032 *>    norm( U'*U - A ) / ( N * norm(A) * EPS ),
00033 *> where EPS is the machine epsilon, L' is the conjugate transpose of
00034 *> L, and U' is the conjugate transpose of U.
00035 *> \endverbatim
00036 *
00037 *  Arguments:
00038 *  ==========
00039 *
00040 *> \param[in] UPLO
00041 *> \verbatim
00042 *>          UPLO is CHARACTER*1
00043 *>          Specifies whether the upper or lower triangular part of the
00044 *>          Hermitian matrix A is stored:
00045 *>          = 'U':  Upper triangular
00046 *>          = 'L':  Lower triangular
00047 *> \endverbatim
00048 *>
00049 *> \param[in] N
00050 *> \verbatim
00051 *>          N is INTEGER
00052 *>          The number of rows and columns of the matrix A.  N >= 0.
00053 *> \endverbatim
00054 *>
00055 *> \param[in] A
00056 *> \verbatim
00057 *>          A is COMPLEX*16 array, dimension (N*(N+1)/2)
00058 *>          The original Hermitian matrix A, stored as a packed
00059 *>          triangular matrix.
00060 *> \endverbatim
00061 *>
00062 *> \param[in,out] AFAC
00063 *> \verbatim
00064 *>          AFAC is COMPLEX*16 array, dimension (N*(N+1)/2)
00065 *>          On entry, the factor L or U from the L*L' or U'*U
00066 *>          factorization of A, stored as a packed triangular matrix.
00067 *>          Overwritten with the reconstructed matrix, and then with the
00068 *>          difference L*L' - A (or U'*U - A).
00069 *> \endverbatim
00070 *>
00071 *> \param[out] RWORK
00072 *> \verbatim
00073 *>          RWORK is DOUBLE PRECISION array, dimension (N)
00074 *> \endverbatim
00075 *>
00076 *> \param[out] RESID
00077 *> \verbatim
00078 *>          RESID is DOUBLE PRECISION
00079 *>          If UPLO = 'L', norm(L*L' - A) / ( N * norm(A) * EPS )
00080 *>          If UPLO = 'U', norm(U'*U - A) / ( N * norm(A) * EPS )
00081 *> \endverbatim
00082 *
00083 *  Authors:
00084 *  ========
00085 *
00086 *> \author Univ. of Tennessee 
00087 *> \author Univ. of California Berkeley 
00088 *> \author Univ. of Colorado Denver 
00089 *> \author NAG Ltd. 
00090 *
00091 *> \date November 2011
00092 *
00093 *> \ingroup complex16_lin
00094 *
00095 *  =====================================================================
00096       SUBROUTINE ZPPT01( UPLO, N, A, AFAC, RWORK, RESID )
00097 *
00098 *  -- LAPACK test routine (version 3.4.0) --
00099 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00100 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00101 *     November 2011
00102 *
00103 *     .. Scalar Arguments ..
00104       CHARACTER          UPLO
00105       INTEGER            N
00106       DOUBLE PRECISION   RESID
00107 *     ..
00108 *     .. Array Arguments ..
00109       DOUBLE PRECISION   RWORK( * )
00110       COMPLEX*16         A( * ), AFAC( * )
00111 *     ..
00112 *
00113 *  =====================================================================
00114 *
00115 *     .. Parameters ..
00116       DOUBLE PRECISION   ZERO, ONE
00117       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
00118 *     ..
00119 *     .. Local Scalars ..
00120       INTEGER            I, K, KC
00121       DOUBLE PRECISION   ANORM, EPS, TR
00122       COMPLEX*16         TC
00123 *     ..
00124 *     .. External Functions ..
00125       LOGICAL            LSAME
00126       DOUBLE PRECISION   DLAMCH, ZLANHP
00127       COMPLEX*16         ZDOTC
00128       EXTERNAL           LSAME, DLAMCH, ZLANHP, ZDOTC
00129 *     ..
00130 *     .. External Subroutines ..
00131       EXTERNAL           ZHPR, ZSCAL, ZTPMV
00132 *     ..
00133 *     .. Intrinsic Functions ..
00134       INTRINSIC          DBLE, DIMAG
00135 *     ..
00136 *     .. Executable Statements ..
00137 *
00138 *     Quick exit if N = 0
00139 *
00140       IF( N.LE.0 ) THEN
00141          RESID = ZERO
00142          RETURN
00143       END IF
00144 *
00145 *     Exit with RESID = 1/EPS if ANORM = 0.
00146 *
00147       EPS = DLAMCH( 'Epsilon' )
00148       ANORM = ZLANHP( '1', UPLO, N, A, RWORK )
00149       IF( ANORM.LE.ZERO ) THEN
00150          RESID = ONE / EPS
00151          RETURN
00152       END IF
00153 *
00154 *     Check the imaginary parts of the diagonal elements and return with
00155 *     an error code if any are nonzero.
00156 *
00157       KC = 1
00158       IF( LSAME( UPLO, 'U' ) ) THEN
00159          DO 10 K = 1, N
00160             IF( DIMAG( AFAC( KC ) ).NE.ZERO ) THEN
00161                RESID = ONE / EPS
00162                RETURN
00163             END IF
00164             KC = KC + K + 1
00165    10    CONTINUE
00166       ELSE
00167          DO 20 K = 1, N
00168             IF( DIMAG( AFAC( KC ) ).NE.ZERO ) THEN
00169                RESID = ONE / EPS
00170                RETURN
00171             END IF
00172             KC = KC + N - K + 1
00173    20    CONTINUE
00174       END IF
00175 *
00176 *     Compute the product U'*U, overwriting U.
00177 *
00178       IF( LSAME( UPLO, 'U' ) ) THEN
00179          KC = ( N*( N-1 ) ) / 2 + 1
00180          DO 30 K = N, 1, -1
00181 *
00182 *           Compute the (K,K) element of the result.
00183 *
00184             TR = ZDOTC( K, AFAC( KC ), 1, AFAC( KC ), 1 )
00185             AFAC( KC+K-1 ) = TR
00186 *
00187 *           Compute the rest of column K.
00188 *
00189             IF( K.GT.1 ) THEN
00190                CALL ZTPMV( 'Upper', 'Conjugate', 'Non-unit', K-1, AFAC,
00191      $                     AFAC( KC ), 1 )
00192                KC = KC - ( K-1 )
00193             END IF
00194    30    CONTINUE
00195 *
00196 *        Compute the difference  L*L' - A
00197 *
00198          KC = 1
00199          DO 50 K = 1, N
00200             DO 40 I = 1, K - 1
00201                AFAC( KC+I-1 ) = AFAC( KC+I-1 ) - A( KC+I-1 )
00202    40       CONTINUE
00203             AFAC( KC+K-1 ) = AFAC( KC+K-1 ) - DBLE( A( KC+K-1 ) )
00204             KC = KC + K
00205    50    CONTINUE
00206 *
00207 *     Compute the product L*L', overwriting L.
00208 *
00209       ELSE
00210          KC = ( N*( N+1 ) ) / 2
00211          DO 60 K = N, 1, -1
00212 *
00213 *           Add a multiple of column K of the factor L to each of
00214 *           columns K+1 through N.
00215 *
00216             IF( K.LT.N )
00217      $         CALL ZHPR( 'Lower', N-K, ONE, AFAC( KC+1 ), 1,
00218      $                    AFAC( KC+N-K+1 ) )
00219 *
00220 *           Scale column K by the diagonal element.
00221 *
00222             TC = AFAC( KC )
00223             CALL ZSCAL( N-K+1, TC, AFAC( KC ), 1 )
00224 *
00225             KC = KC - ( N-K+2 )
00226    60    CONTINUE
00227 *
00228 *        Compute the difference  U'*U - A
00229 *
00230          KC = 1
00231          DO 80 K = 1, N
00232             AFAC( KC ) = AFAC( KC ) - DBLE( A( KC ) )
00233             DO 70 I = K + 1, N
00234                AFAC( KC+I-K ) = AFAC( KC+I-K ) - A( KC+I-K )
00235    70       CONTINUE
00236             KC = KC + N - K + 1
00237    80    CONTINUE
00238       END IF
00239 *
00240 *     Compute norm( L*U - A ) / ( N * norm(A) * EPS )
00241 *
00242       RESID = ZLANHP( '1', UPLO, N, AFAC, RWORK )
00243 *
00244       RESID = ( ( RESID / DBLE( N ) ) / ANORM ) / EPS
00245 *
00246       RETURN
00247 *
00248 *     End of ZPPT01
00249 *
00250       END
 All Files Functions