LAPACK 3.3.1 Linear Algebra PACKage

# dppt01.f

Go to the documentation of this file.
```00001       SUBROUTINE DPPT01( 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   A( * ), AFAC( * ), RWORK( * )
00014 *     ..
00015 *
00016 *  Purpose
00017 *  =======
00018 *
00019 *  DPPT01 reconstructs a symmetric positive definite packed matrix A
00020 *  from its L*L' or U'*U factorization and computes the residual
00021 *     norm( L*L' - A ) / ( N * norm(A) * EPS ) or
00022 *     norm( U'*U - A ) / ( N * norm(A) * EPS ),
00023 *  where EPS is the machine epsilon.
00024 *
00025 *  Arguments
00026 *  ==========
00027 *
00028 *  UPLO    (input) CHARACTER*1
00029 *          Specifies whether the upper or lower triangular part of the
00030 *          symmetric matrix A is stored:
00031 *          = 'U':  Upper triangular
00032 *          = 'L':  Lower triangular
00033 *
00034 *  N       (input) INTEGER
00035 *          The number of rows and columns of the matrix A.  N >= 0.
00036 *
00037 *  A       (input) DOUBLE PRECISION array, dimension (N*(N+1)/2)
00038 *          The original symmetric matrix A, stored as a packed
00039 *          triangular matrix.
00040 *
00041 *  AFAC    (input/output) DOUBLE PRECISION array, dimension (N*(N+1)/2)
00042 *          On entry, the factor L or U from the L*L' or U'*U
00043 *          factorization of A, stored as a packed triangular matrix.
00044 *          Overwritten with the reconstructed matrix, and then with the
00045 *          difference L*L' - A (or U'*U - A).
00046 *
00047 *  RWORK   (workspace) DOUBLE PRECISION array, dimension (N)
00048 *
00049 *  RESID   (output) DOUBLE PRECISION
00050 *          If UPLO = 'L', norm(L*L' - A) / ( N * norm(A) * EPS )
00051 *          If UPLO = 'U', norm(U'*U - A) / ( N * norm(A) * EPS )
00052 *
00053 *  =====================================================================
00054 *
00055 *     .. Parameters ..
00056       DOUBLE PRECISION   ZERO, ONE
00057       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
00058 *     ..
00059 *     .. Local Scalars ..
00060       INTEGER            I, K, KC, NPP
00061       DOUBLE PRECISION   ANORM, EPS, T
00062 *     ..
00063 *     .. External Functions ..
00064       LOGICAL            LSAME
00065       DOUBLE PRECISION   DDOT, DLAMCH, DLANSP
00066       EXTERNAL           LSAME, DDOT, DLAMCH, DLANSP
00067 *     ..
00068 *     .. External Subroutines ..
00069       EXTERNAL           DSCAL, DSPR, DTPMV
00070 *     ..
00071 *     .. Intrinsic Functions ..
00072       INTRINSIC          DBLE
00073 *     ..
00074 *     .. Executable Statements ..
00075 *
00076 *     Quick exit if N = 0
00077 *
00078       IF( N.LE.0 ) THEN
00079          RESID = ZERO
00080          RETURN
00081       END IF
00082 *
00083 *     Exit with RESID = 1/EPS if ANORM = 0.
00084 *
00085       EPS = DLAMCH( 'Epsilon' )
00086       ANORM = DLANSP( '1', UPLO, N, A, RWORK )
00087       IF( ANORM.LE.ZERO ) THEN
00088          RESID = ONE / EPS
00089          RETURN
00090       END IF
00091 *
00092 *     Compute the product U'*U, overwriting U.
00093 *
00094       IF( LSAME( UPLO, 'U' ) ) THEN
00095          KC = ( N*( N-1 ) ) / 2 + 1
00096          DO 10 K = N, 1, -1
00097 *
00098 *           Compute the (K,K) element of the result.
00099 *
00100             T = DDOT( K, AFAC( KC ), 1, AFAC( KC ), 1 )
00101             AFAC( KC+K-1 ) = T
00102 *
00103 *           Compute the rest of column K.
00104 *
00105             IF( K.GT.1 ) THEN
00106                CALL DTPMV( 'Upper', 'Transpose', 'Non-unit', K-1, AFAC,
00107      \$                     AFAC( KC ), 1 )
00108                KC = KC - ( K-1 )
00109             END IF
00110    10    CONTINUE
00111 *
00112 *     Compute the product L*L', overwriting L.
00113 *
00114       ELSE
00115          KC = ( N*( N+1 ) ) / 2
00116          DO 20 K = N, 1, -1
00117 *
00118 *           Add a multiple of column K of the factor L to each of
00119 *           columns K+1 through N.
00120 *
00121             IF( K.LT.N )
00122      \$         CALL DSPR( 'Lower', N-K, ONE, AFAC( KC+1 ), 1,
00123      \$                    AFAC( KC+N-K+1 ) )
00124 *
00125 *           Scale column K by the diagonal element.
00126 *
00127             T = AFAC( KC )
00128             CALL DSCAL( N-K+1, T, AFAC( KC ), 1 )
00129 *
00130             KC = KC - ( N-K+2 )
00131    20    CONTINUE
00132       END IF
00133 *
00134 *     Compute the difference  L*L' - A (or U'*U - A).
00135 *
00136       NPP = N*( N+1 ) / 2
00137       DO 30 I = 1, NPP
00138          AFAC( I ) = AFAC( I ) - A( I )
00139    30 CONTINUE
00140 *
00141 *     Compute norm( L*U - A ) / ( N * norm(A) * EPS )
00142 *
00143       RESID = DLANSP( '1', UPLO, N, AFAC, RWORK )
00144 *
00145       RESID = ( ( RESID / DBLE( N ) ) / ANORM ) / EPS
00146 *
00147       RETURN
00148 *
00149 *     End of DPPT01
00150 *
00151       END
```