LAPACK: TESTING/LIN/dpst01.f Source File

Go to the documentation of this file.
00001       SUBROUTINE DPST01( UPLO, N, A, LDA, AFAC, LDAFAC, PERM, LDPERM,
00002      $                   PIV, RWORK, RESID, RANK )
00003 *
00004 *  -- LAPACK test routine (version 3.1) --
00005 *     Craig Lucas, University of Manchester / NAG Ltd.
00006 *     October, 2008
00007 *
00008 *     .. Scalar Arguments ..
00009       DOUBLE PRECISION   RESID
00010       INTEGER            LDA, LDAFAC, LDPERM, N, RANK
00011       CHARACTER          UPLO
00012 *     ..
00013 *     .. Array Arguments ..
00014       DOUBLE PRECISION   A( LDA, * ), AFAC( LDAFAC, * ),
00015      $                   PERM( LDPERM, * ), RWORK( * )
00016       INTEGER            PIV( * )
00017 *     ..
00018 *
00019 *  Purpose
00020 *  =======
00021 *
00022 *  DPST01 reconstructs a symmetric positive semidefinite matrix A
00023 *  from its L or U factors and the permutation matrix P and computes
00024 *  the residual
00025 *     norm( P*L*L'*P' - A ) / ( N * norm(A) * EPS ) or
00026 *     norm( P*U'*U*P' - A ) / ( N * norm(A) * EPS ),
00027 *  where EPS is the machine epsilon.
00028 *
00029 *  Arguments
00030 *  ==========
00031 *
00032 *  UPLO    (input) CHARACTER*1
00033 *          Specifies whether the upper or lower triangular part of the
00034 *          symmetric matrix A is stored:
00035 *          = 'U':  Upper triangular
00036 *          = 'L':  Lower triangular
00037 *
00038 *  N       (input) INTEGER
00039 *          The number of rows and columns of the matrix A.  N >= 0.
00040 *
00041 *  A       (input) DOUBLE PRECISION array, dimension (LDA,N)
00042 *          The original symmetric matrix A.
00043 *
00044 *  LDA     (input) INTEGER
00045 *          The leading dimension of the array A.  LDA >= max(1,N)
00046 *
00047 *  AFAC    (input) DOUBLE PRECISION array, dimension (LDAFAC,N)
00048 *          The factor L or U from the L*L' or U'*U
00049 *          factorization of A.
00050 *
00051 *  LDAFAC  (input) INTEGER
00052 *          The leading dimension of the array AFAC.  LDAFAC >= max(1,N).
00053 *
00054 *  PERM    (output) DOUBLE PRECISION array, dimension (LDPERM,N)
00055 *          Overwritten with the reconstructed matrix, and then with the
00056 *          difference P*L*L'*P' - A (or P*U'*U*P' - A)
00057 *
00058 *  LDPERM  (input) INTEGER
00059 *          The leading dimension of the array PERM.
00060 *          LDAPERM >= max(1,N).
00061 *
00062 *  PIV     (input) INTEGER array, dimension (N)
00063 *          PIV is such that the nonzero entries are
00064 *          P( PIV( K ), K ) = 1.
00065 *
00066 *  RWORK   (workspace) DOUBLE PRECISION array, dimension (N)
00067 *
00068 *  RESID   (output) DOUBLE PRECISION
00069 *          If UPLO = 'L', norm(L*L' - A) / ( N * norm(A) * EPS )
00070 *          If UPLO = 'U', norm(U'*U - A) / ( N * norm(A) * EPS )
00071 *
00072 *  =====================================================================
00073 *
00074 *     .. Parameters ..
00075       DOUBLE PRECISION   ZERO, ONE
00076       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
00077 *     ..
00078 *     .. Local Scalars ..
00079       DOUBLE PRECISION   ANORM, EPS, T
00080       INTEGER            I, J, K
00081 *     ..
00082 *     .. External Functions ..
00083       DOUBLE PRECISION   DDOT, DLAMCH, DLANSY
00084       LOGICAL            LSAME
00085       EXTERNAL           DDOT, DLAMCH, DLANSY, LSAME
00086 *     ..
00087 *     .. External Subroutines ..
00088       EXTERNAL           DSCAL, DSYR, DTRMV
00089 *     ..
00090 *     .. Intrinsic Functions ..
00091       INTRINSIC          DBLE
00092 *     ..
00093 *     .. Executable Statements ..
00094 *
00095 *     Quick exit if N = 0.
00096 *
00097       IF( N.LE.0 ) THEN
00098          RESID = ZERO
00099          RETURN
00100       END IF
00101 *
00102 *     Exit with RESID = 1/EPS if ANORM = 0.
00103 *
00104       EPS = DLAMCH( 'Epsilon' )
00105       ANORM = DLANSY( '1', UPLO, N, A, LDA, RWORK )
00106       IF( ANORM.LE.ZERO ) THEN
00107          RESID = ONE / EPS
00108          RETURN
00109       END IF
00110 *
00111 *     Compute the product U'*U, overwriting U.
00112 *
00113       IF( LSAME( UPLO, 'U' ) ) THEN
00114 *
00115          IF( RANK.LT.N ) THEN
00116             DO 110 J = RANK + 1, N
00117                DO 100 I = RANK + 1, J
00118                   AFAC( I, J ) = ZERO
00119   100          CONTINUE
00120   110       CONTINUE
00121          END IF
00122 *
00123          DO 120 K = N, 1, -1
00124 *
00125 *           Compute the (K,K) element of the result.
00126 *
00127             T = DDOT( K, AFAC( 1, K ), 1, AFAC( 1, K ), 1 )
00128             AFAC( K, K ) = T
00129 *
00130 *           Compute the rest of column K.
00131 *
00132             CALL DTRMV( 'Upper', 'Transpose', 'Non-unit', K-1, AFAC,
00133      $                  LDAFAC, AFAC( 1, K ), 1 )
00134 *
00135   120    CONTINUE
00136 *
00137 *     Compute the product L*L', overwriting L.
00138 *
00139       ELSE
00140 *
00141          IF( RANK.LT.N ) THEN
00142             DO 140 J = RANK + 1, N
00143                DO 130 I = J, N
00144                   AFAC( I, J ) = ZERO
00145   130          CONTINUE
00146   140       CONTINUE
00147          END IF
00148 *
00149          DO 150 K = N, 1, -1
00150 *           Add a multiple of column K of the factor L to each of
00151 *           columns K+1 through N.
00152 *
00153             IF( K+1.LE.N )
00154      $         CALL DSYR( 'Lower', N-K, ONE, AFAC( K+1, K ), 1,
00155      $                    AFAC( K+1, K+1 ), LDAFAC )
00156 *
00157 *           Scale column K by the diagonal element.
00158 *
00159             T = AFAC( K, K )
00160             CALL DSCAL( N-K+1, T, AFAC( K, K ), 1 )
00161   150    CONTINUE
00162 *
00163       END IF
00164 *
00165 *        Form P*L*L'*P' or P*U'*U*P'
00166 *
00167       IF( LSAME( UPLO, 'U' ) ) THEN
00168 *
00169          DO 170 J = 1, N
00170             DO 160 I = 1, N
00171                IF( PIV( I ).LE.PIV( J ) ) THEN
00172                   IF( I.LE.J ) THEN
00173                      PERM( PIV( I ), PIV( J ) ) = AFAC( I, J )
00174                   ELSE
00175                      PERM( PIV( I ), PIV( J ) ) = AFAC( J, I )
00176                   END IF
00177                END IF
00178   160       CONTINUE
00179   170    CONTINUE
00180 *
00181 *
00182       ELSE
00183 *
00184          DO 190 J = 1, N
00185             DO 180 I = 1, N
00186                IF( PIV( I ).GE.PIV( J ) ) THEN
00187                   IF( I.GE.J ) THEN
00188                      PERM( PIV( I ), PIV( J ) ) = AFAC( I, J )
00189                   ELSE
00190                      PERM( PIV( I ), PIV( J ) ) = AFAC( J, I )
00191                   END IF
00192                END IF
00193   180       CONTINUE
00194   190    CONTINUE
00195 *
00196       END IF
00197 *
00198 *     Compute the difference  P*L*L'*P' - A (or P*U'*U*P' - A).
00199 *
00200       IF( LSAME( UPLO, 'U' ) ) THEN
00201          DO 210 J = 1, N
00202             DO 200 I = 1, J
00203                PERM( I, J ) = PERM( I, J ) - A( I, J )
00204   200       CONTINUE
00205   210    CONTINUE
00206       ELSE
00207          DO 230 J = 1, N
00208             DO 220 I = J, N
00209                PERM( I, J ) = PERM( I, J ) - A( I, J )
00210   220       CONTINUE
00211   230    CONTINUE
00212       END IF
00213 *
00214 *     Compute norm( P*L*L'P - A ) / ( N * norm(A) * EPS ), or
00215 *     ( P*U'*U*P' - A )/ ( N * norm(A) * EPS ).
00216 *
00217       RESID = DLANSY( '1', UPLO, N, PERM, LDAFAC, RWORK )
00218 *
00219       RESID = ( ( RESID / DBLE( N ) ) / ANORM ) / EPS
00220 *
00221       RETURN
00222 *
00223 *     End of DPST01
00224 *
00225       END