LAPACK: TESTING/LIN/cpst01.f Source File

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