LAPACK 3.3.0

dpot01.f

Go to the documentation of this file.
00001       SUBROUTINE DPOT01( 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   A( LDA, * ), AFAC( LDAFAC, * ), RWORK( * )
00014 *     ..
00015 *
00016 *  Purpose
00017 *  =======
00018 *
00019 *  DPOT01 reconstructs a symmetric positive definite matrix  A  from
00020 *  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 (LDA,N)
00038 *          The original symmetric matrix A.
00039 *
00040 *  LDA     (input) INTEGER
00041 *          The leading dimension of the array A.  LDA >= max(1,N)
00042 *
00043 *  AFAC    (input/output) DOUBLE PRECISION array, dimension (LDAFAC,N)
00044 *          On entry, the factor L or U from the L*L' or U'*U
00045 *          factorization of A.
00046 *          Overwritten with the reconstructed matrix, and then with the
00047 *          difference L*L' - A (or U'*U - A).
00048 *
00049 *  LDAFAC  (input) INTEGER
00050 *          The leading dimension of the array AFAC.  LDAFAC >= max(1,N).
00051 *
00052 *  RWORK   (workspace) DOUBLE PRECISION array, dimension (N)
00053 *
00054 *  RESID   (output) DOUBLE PRECISION
00055 *          If UPLO = 'L', norm(L*L' - A) / ( N * norm(A) * EPS )
00056 *          If UPLO = 'U', norm(U'*U - A) / ( N * norm(A) * EPS )
00057 *
00058 *  =====================================================================
00059 *
00060 *     .. Parameters ..
00061       DOUBLE PRECISION   ZERO, ONE
00062       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
00063 *     ..
00064 *     .. Local Scalars ..
00065       INTEGER            I, J, K
00066       DOUBLE PRECISION   ANORM, EPS, T
00067 *     ..
00068 *     .. External Functions ..
00069       LOGICAL            LSAME
00070       DOUBLE PRECISION   DDOT, DLAMCH, DLANSY
00071       EXTERNAL           LSAME, DDOT, DLAMCH, DLANSY
00072 *     ..
00073 *     .. External Subroutines ..
00074       EXTERNAL           DSCAL, DSYR, DTRMV
00075 *     ..
00076 *     .. Intrinsic Functions ..
00077       INTRINSIC          DBLE
00078 *     ..
00079 *     .. Executable Statements ..
00080 *
00081 *     Quick exit if N = 0.
00082 *
00083       IF( N.LE.0 ) THEN
00084          RESID = ZERO
00085          RETURN
00086       END IF
00087 *
00088 *     Exit with RESID = 1/EPS if ANORM = 0.
00089 *
00090       EPS = DLAMCH( 'Epsilon' )
00091       ANORM = DLANSY( '1', UPLO, N, A, LDA, RWORK )
00092       IF( ANORM.LE.ZERO ) THEN
00093          RESID = ONE / EPS
00094          RETURN
00095       END IF
00096 *
00097 *     Compute the product U'*U, overwriting U.
00098 *
00099       IF( LSAME( UPLO, 'U' ) ) THEN
00100          DO 10 K = N, 1, -1
00101 *
00102 *           Compute the (K,K) element of the result.
00103 *
00104             T = DDOT( K, AFAC( 1, K ), 1, AFAC( 1, K ), 1 )
00105             AFAC( K, K ) = T
00106 *
00107 *           Compute the rest of column K.
00108 *
00109             CALL DTRMV( 'Upper', 'Transpose', 'Non-unit', K-1, AFAC,
00110      $                  LDAFAC, AFAC( 1, K ), 1 )
00111 *
00112    10    CONTINUE
00113 *
00114 *     Compute the product L*L', overwriting L.
00115 *
00116       ELSE
00117          DO 20 K = N, 1, -1
00118 *
00119 *           Add a multiple of column K of the factor L to each of
00120 *           columns K+1 through N.
00121 *
00122             IF( K+1.LE.N )
00123      $         CALL DSYR( 'Lower', N-K, ONE, AFAC( K+1, K ), 1,
00124      $                    AFAC( K+1, K+1 ), LDAFAC )
00125 *
00126 *           Scale column K by the diagonal element.
00127 *
00128             T = AFAC( K, K )
00129             CALL DSCAL( N-K+1, T, AFAC( K, K ), 1 )
00130 *
00131    20    CONTINUE
00132       END IF
00133 *
00134 *     Compute the difference  L*L' - A (or U'*U - A).
00135 *
00136       IF( LSAME( UPLO, 'U' ) ) THEN
00137          DO 40 J = 1, N
00138             DO 30 I = 1, J
00139                AFAC( I, J ) = AFAC( I, J ) - A( I, J )
00140    30       CONTINUE
00141    40    CONTINUE
00142       ELSE
00143          DO 60 J = 1, N
00144             DO 50 I = J, N
00145                AFAC( I, J ) = AFAC( I, J ) - A( I, J )
00146    50       CONTINUE
00147    60    CONTINUE
00148       END IF
00149 *
00150 *     Compute norm( L*U - A ) / ( N * norm(A) * EPS )
00151 *
00152       RESID = DLANSY( '1', UPLO, N, AFAC, LDAFAC, RWORK )
00153 *
00154       RESID = ( ( RESID / DBLE( N ) ) / ANORM ) / EPS
00155 *
00156       RETURN
00157 *
00158 *     End of DPOT01
00159 *
00160       END
 All Files Functions