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