LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
cpst01.f
Go to the documentation of this file.
00001 *> \brief \b CPST01
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 CPST01( UPLO, N, A, LDA, AFAC, LDAFAC, PERM, LDPERM,
00012 *                          PIV, RWORK, RESID, RANK )
00013 * 
00014 *       .. Scalar Arguments ..
00015 *       REAL               RESID
00016 *       INTEGER            LDA, LDAFAC, LDPERM, N, RANK
00017 *       CHARACTER          UPLO
00018 *       ..
00019 *       .. Array Arguments ..
00020 *       COMPLEX            A( LDA, * ), AFAC( LDAFAC, * ),
00021 *      $                   PERM( LDPERM, * )
00022 *       REAL               RWORK( * )
00023 *       INTEGER            PIV( * )
00024 *       ..
00025 *  
00026 *
00027 *> \par Purpose:
00028 *  =============
00029 *>
00030 *> \verbatim
00031 *>
00032 *> CPST01 reconstructs an Hermitian positive semidefinite matrix A
00033 *> from its L or U factors and the permutation matrix P and computes
00034 *> the residual
00035 *>    norm( P*L*L'*P' - A ) / ( N * norm(A) * EPS ) or
00036 *>    norm( P*U'*U*P' - A ) / ( N * norm(A) * EPS ),
00037 *> where EPS is the machine epsilon, L' is the conjugate transpose of L,
00038 *> and U' is the conjugate transpose of U.
00039 *> \endverbatim
00040 *
00041 *  Arguments:
00042 *  ==========
00043 *
00044 *> \param[in] UPLO
00045 *> \verbatim
00046 *>          UPLO is CHARACTER*1
00047 *>          Specifies whether the upper or lower triangular part of the
00048 *>          Hermitian matrix A is stored:
00049 *>          = 'U':  Upper triangular
00050 *>          = 'L':  Lower triangular
00051 *> \endverbatim
00052 *>
00053 *> \param[in] N
00054 *> \verbatim
00055 *>          N is INTEGER
00056 *>          The number of rows and columns of the matrix A.  N >= 0.
00057 *> \endverbatim
00058 *>
00059 *> \param[in] A
00060 *> \verbatim
00061 *>          A is COMPLEX array, dimension (LDA,N)
00062 *>          The original Hermitian matrix A.
00063 *> \endverbatim
00064 *>
00065 *> \param[in] LDA
00066 *> \verbatim
00067 *>          LDA is INTEGER
00068 *>          The leading dimension of the array A.  LDA >= max(1,N)
00069 *> \endverbatim
00070 *>
00071 *> \param[in] AFAC
00072 *> \verbatim
00073 *>          AFAC is COMPLEX array, dimension (LDAFAC,N)
00074 *>          The factor L or U from the L*L' or U'*U
00075 *>          factorization of A.
00076 *> \endverbatim
00077 *>
00078 *> \param[in] LDAFAC
00079 *> \verbatim
00080 *>          LDAFAC is INTEGER
00081 *>          The leading dimension of the array AFAC.  LDAFAC >= max(1,N).
00082 *> \endverbatim
00083 *>
00084 *> \param[out] PERM
00085 *> \verbatim
00086 *>          PERM is COMPLEX array, dimension (LDPERM,N)
00087 *>          Overwritten with the reconstructed matrix, and then with the
00088 *>          difference P*L*L'*P' - A (or P*U'*U*P' - A)
00089 *> \endverbatim
00090 *>
00091 *> \param[in] LDPERM
00092 *> \verbatim
00093 *>          LDPERM is INTEGER
00094 *>          The leading dimension of the array PERM.
00095 *>          LDAPERM >= max(1,N).
00096 *> \endverbatim
00097 *>
00098 *> \param[in] PIV
00099 *> \verbatim
00100 *>          PIV is INTEGER array, dimension (N)
00101 *>          PIV is such that the nonzero entries are
00102 *>          P( PIV( K ), K ) = 1.
00103 *> \endverbatim
00104 *>
00105 *> \param[out] RWORK
00106 *> \verbatim
00107 *>          RWORK is REAL array, dimension (N)
00108 *> \endverbatim
00109 *>
00110 *> \param[out] RESID
00111 *> \verbatim
00112 *>          RESID is REAL
00113 *>          If UPLO = 'L', norm(L*L' - A) / ( N * norm(A) * EPS )
00114 *>          If UPLO = 'U', norm(U'*U - A) / ( N * norm(A) * EPS )
00115 *> \endverbatim
00116 *>
00117 *> \param[in] RANK
00118 *> \verbatim
00119 *>          RANK is INTEGER
00120 *>          number of nonzero singular values of A.
00121 *> \endverbatim
00122 *
00123 *  Authors:
00124 *  ========
00125 *
00126 *> \author Univ. of Tennessee 
00127 *> \author Univ. of California Berkeley 
00128 *> \author Univ. of Colorado Denver 
00129 *> \author NAG Ltd. 
00130 *
00131 *> \date November 2011
00132 *
00133 *> \ingroup complex_lin
00134 *
00135 *  =====================================================================
00136       SUBROUTINE CPST01( UPLO, N, A, LDA, AFAC, LDAFAC, PERM, LDPERM,
00137      $                   PIV, RWORK, RESID, RANK )
00138 *
00139 *  -- LAPACK test routine (version 3.4.0) --
00140 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00141 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00142 *     November 2011
00143 *
00144 *     .. Scalar Arguments ..
00145       REAL               RESID
00146       INTEGER            LDA, LDAFAC, LDPERM, N, RANK
00147       CHARACTER          UPLO
00148 *     ..
00149 *     .. Array Arguments ..
00150       COMPLEX            A( LDA, * ), AFAC( LDAFAC, * ),
00151      $                   PERM( LDPERM, * )
00152       REAL               RWORK( * )
00153       INTEGER            PIV( * )
00154 *     ..
00155 *
00156 *  =====================================================================
00157 *
00158 *     .. Parameters ..
00159       REAL               ZERO, ONE
00160       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
00161       COMPLEX            CZERO
00162       PARAMETER          ( CZERO = ( 0.0E+0, 0.0E+0 ) )
00163 *     ..
00164 *     .. Local Scalars ..
00165       COMPLEX            TC
00166       REAL               ANORM, EPS, TR
00167       INTEGER            I, J, K
00168 *     ..
00169 *     .. External Functions ..
00170       COMPLEX            CDOTC
00171       REAL               CLANHE, SLAMCH
00172       LOGICAL            LSAME
00173       EXTERNAL           CDOTC, CLANHE, SLAMCH, LSAME
00174 *     ..
00175 *     .. External Subroutines ..
00176       EXTERNAL           CHER, CSCAL, CTRMV
00177 *     ..
00178 *     .. Intrinsic Functions ..
00179       INTRINSIC          AIMAG, CONJG, REAL
00180 *     ..
00181 *     .. Executable Statements ..
00182 *
00183 *     Quick exit if N = 0.
00184 *
00185       IF( N.LE.0 ) THEN
00186          RESID = ZERO
00187          RETURN
00188       END IF
00189 *
00190 *     Exit with RESID = 1/EPS if ANORM = 0.
00191 *
00192       EPS = SLAMCH( 'Epsilon' )
00193       ANORM = CLANHE( '1', UPLO, N, A, LDA, RWORK )
00194       IF( ANORM.LE.ZERO ) THEN
00195          RESID = ONE / EPS
00196          RETURN
00197       END IF
00198 *
00199 *     Check the imaginary parts of the diagonal elements and return with
00200 *     an error code if any are nonzero.
00201 *
00202       DO 100 J = 1, N
00203          IF( AIMAG( AFAC( J, J ) ).NE.ZERO ) THEN
00204             RESID = ONE / EPS
00205             RETURN
00206          END IF
00207   100 CONTINUE
00208 *
00209 *     Compute the product U'*U, overwriting U.
00210 *
00211       IF( LSAME( UPLO, 'U' ) ) THEN
00212 *
00213          IF( RANK.LT.N ) THEN
00214             DO 120 J = RANK + 1, N
00215                DO 110 I = RANK + 1, J
00216                   AFAC( I, J ) = CZERO
00217   110          CONTINUE
00218   120       CONTINUE
00219          END IF
00220 *
00221          DO 130 K = N, 1, -1
00222 *
00223 *           Compute the (K,K) element of the result.
00224 *
00225             TR = CDOTC( K, AFAC( 1, K ), 1, AFAC( 1, K ), 1 )
00226             AFAC( K, K ) = TR
00227 *
00228 *           Compute the rest of column K.
00229 *
00230             CALL CTRMV( 'Upper', 'Conjugate', 'Non-unit', K-1, AFAC,
00231      $                  LDAFAC, AFAC( 1, K ), 1 )
00232 *
00233   130    CONTINUE
00234 *
00235 *     Compute the product L*L', overwriting L.
00236 *
00237       ELSE
00238 *
00239          IF( RANK.LT.N ) THEN
00240             DO 150 J = RANK + 1, N
00241                DO 140 I = J, N
00242                   AFAC( I, J ) = CZERO
00243   140          CONTINUE
00244   150       CONTINUE
00245          END IF
00246 *
00247          DO 160 K = N, 1, -1
00248 *           Add a multiple of column K of the factor L to each of
00249 *           columns K+1 through N.
00250 *
00251             IF( K+1.LE.N )
00252      $         CALL CHER( 'Lower', N-K, ONE, AFAC( K+1, K ), 1,
00253      $                    AFAC( K+1, K+1 ), LDAFAC )
00254 *
00255 *           Scale column K by the diagonal element.
00256 *
00257             TC = AFAC( K, K )
00258             CALL CSCAL( N-K+1, TC, AFAC( K, K ), 1 )
00259   160    CONTINUE
00260 *
00261       END IF
00262 *
00263 *        Form P*L*L'*P' or P*U'*U*P'
00264 *
00265       IF( LSAME( UPLO, 'U' ) ) THEN
00266 *
00267          DO 180 J = 1, N
00268             DO 170 I = 1, N
00269                IF( PIV( I ).LE.PIV( J ) ) THEN
00270                   IF( I.LE.J ) THEN
00271                      PERM( PIV( I ), PIV( J ) ) = AFAC( I, J )
00272                   ELSE
00273                      PERM( PIV( I ), PIV( J ) ) = CONJG( AFAC( J, I ) )
00274                   END IF
00275                END IF
00276   170       CONTINUE
00277   180    CONTINUE
00278 *
00279 *
00280       ELSE
00281 *
00282          DO 200 J = 1, N
00283             DO 190 I = 1, N
00284                IF( PIV( I ).GE.PIV( J ) ) THEN
00285                   IF( I.GE.J ) THEN
00286                      PERM( PIV( I ), PIV( J ) ) = AFAC( I, J )
00287                   ELSE
00288                      PERM( PIV( I ), PIV( J ) ) = CONJG( AFAC( J, I ) )
00289                   END IF
00290                END IF
00291   190       CONTINUE
00292   200    CONTINUE
00293 *
00294       END IF
00295 *
00296 *     Compute the difference  P*L*L'*P' - A (or P*U'*U*P' - A).
00297 *
00298       IF( LSAME( UPLO, 'U' ) ) THEN
00299          DO 220 J = 1, N
00300             DO 210 I = 1, J - 1
00301                PERM( I, J ) = PERM( I, J ) - A( I, J )
00302   210       CONTINUE
00303             PERM( J, J ) = PERM( J, J ) - REAL( A( J, J ) )
00304   220    CONTINUE
00305       ELSE
00306          DO 240 J = 1, N
00307             PERM( J, J ) = PERM( J, J ) - REAL( A( J, J ) )
00308             DO 230 I = J + 1, N
00309                PERM( I, J ) = PERM( I, J ) - A( I, J )
00310   230       CONTINUE
00311   240    CONTINUE
00312       END IF
00313 *
00314 *     Compute norm( P*L*L'P - A ) / ( N * norm(A) * EPS ), or
00315 *     ( P*U'*U*P' - A )/ ( N * norm(A) * EPS ).
00316 *
00317       RESID = CLANHE( '1', UPLO, N, PERM, LDAFAC, RWORK )
00318 *
00319       RESID = ( ( RESID / REAL( N ) ) / ANORM ) / EPS
00320 *
00321       RETURN
00322 *
00323 *     End of CPST01
00324 *
00325       END
 All Files Functions