TESTING/LIN/spst01.f

   1 *> \brief \b SPST01
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *  Definition:
   9 *  ===========
  10 *
  11 *       SUBROUTINE SPST01( UPLO, N, A, LDA, AFAC, LDAFAC, PERM, LDPERM,
  12 *                          PIV, RWORK, RESID, RANK )
  13 *
  14 *       .. Scalar Arguments ..
  15 *       REAL               RESID
  16 *       INTEGER            LDA, LDAFAC, LDPERM, N, RANK
  17 *       CHARACTER          UPLO
  18 *       ..
  19 *       .. Array Arguments ..
  20 *       REAL               A( LDA, * ), AFAC( LDAFAC, * ),
  21 *      $                   PERM( LDPERM, * ), RWORK( * )
  22 *       INTEGER            PIV( * )
  23 *       ..
  24 *
  25 *
  26 *> \par Purpose:
  27 *  =============
  28 *>
  29 *> \verbatim
  30 *>
  31 *> SPST01 reconstructs a symmetric positive semidefinite matrix A
  32 *> from its L or U factors and the permutation matrix P and computes
  33 *> the residual
  34 *>    norm( P*L*L'*P' - A ) / ( N * norm(A) * EPS ) or
  35 *>    norm( P*U'*U*P' - A ) / ( N * norm(A) * EPS ),
  36 *> where EPS is the machine epsilon.
  37 *> \endverbatim
  38 *
  39 *  Arguments:
  40 *  ==========
  41 *
  42 *> \param[in] UPLO
  43 *> \verbatim
  44 *>          UPLO is CHARACTER*1
  45 *>          Specifies whether the upper or lower triangular part of the
  46 *>          symmetric matrix A is stored:
  47 *>          = 'U':  Upper triangular
  48 *>          = 'L':  Lower triangular
  49 *> \endverbatim
  50 *>
  51 *> \param[in] N
  52 *> \verbatim
  53 *>          N is INTEGER
  54 *>          The number of rows and columns of the matrix A.  N >= 0.
  55 *> \endverbatim
  56 *>
  57 *> \param[in] A
  58 *> \verbatim
  59 *>          A is REAL array, dimension (LDA,N)
  60 *>          The original symmetric matrix A.
  61 *> \endverbatim
  62 *>
  63 *> \param[in] LDA
  64 *> \verbatim
  65 *>          LDA is INTEGER
  66 *>          The leading dimension of the array A.  LDA >= max(1,N)
  67 *> \endverbatim
  68 *>
  69 *> \param[in] AFAC
  70 *> \verbatim
  71 *>          AFAC is REAL array, dimension (LDAFAC,N)
  72 *>          The factor L or U from the L*L' or U'*U
  73 *>          factorization of A.
  74 *> \endverbatim
  75 *>
  76 *> \param[in] LDAFAC
  77 *> \verbatim
  78 *>          LDAFAC is INTEGER
  79 *>          The leading dimension of the array AFAC.  LDAFAC >= max(1,N).
  80 *> \endverbatim
  81 *>
  82 *> \param[out] PERM
  83 *> \verbatim
  84 *>          PERM is REAL array, dimension (LDPERM,N)
  85 *>          Overwritten with the reconstructed matrix, and then with the
  86 *>          difference P*L*L'*P' - A (or P*U'*U*P' - A)
  87 *> \endverbatim
  88 *>
  89 *> \param[in] LDPERM
  90 *> \verbatim
  91 *>          LDPERM is INTEGER
  92 *>          The leading dimension of the array PERM.
  93 *>          LDAPERM >= max(1,N).
  94 *> \endverbatim
  95 *>
  96 *> \param[in] PIV
  97 *> \verbatim
  98 *>          PIV is INTEGER array, dimension (N)
  99 *>          PIV is such that the nonzero entries are
 100 *>          P( PIV( K ), K ) = 1.
 101 *> \endverbatim
 102 *>
 103 *> \param[out] RWORK
 104 *> \verbatim
 105 *>          RWORK is REAL array, dimension (N)
 106 *> \endverbatim
 107 *>
 108 *> \param[out] RESID
 109 *> \verbatim
 110 *>          RESID is REAL
 111 *>          If UPLO = 'L', norm(L*L' - A) / ( N * norm(A) * EPS )
 112 *>          If UPLO = 'U', norm(U'*U - A) / ( N * norm(A) * EPS )
 113 *> \endverbatim
 114 *>
 115 *> \param[in] RANK
 116 *> \verbatim
 117 *>          RANK is INTEGER
 118 *>          number of nonzero singular values of A.
 119 *> \endverbatim
 120 *
 121 *  Authors:
 122 *  ========
 123 *
 124 *> \author Univ. of Tennessee
 125 *> \author Univ. of California Berkeley
 126 *> \author Univ. of Colorado Denver
 127 *> \author NAG Ltd.
 128 *
 129 *> \date November 2011
 130 *
 131 *> \ingroup single_lin
 132 *
 133 *  =====================================================================
 134       SUBROUTINE SPST01( UPLO, N, A, LDA, AFAC, LDAFAC, PERM, LDPERM,
 135      $                   PIV, RWORK, RESID, RANK )
 136 *
 137 *  -- LAPACK test routine (version 3.4.0) --
 138 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 139 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 140 *     November 2011
 141 *
 142 *     .. Scalar Arguments ..
 143       REAL               RESID
 144       INTEGER            LDA, LDAFAC, LDPERM, N, RANK
 145       CHARACTER          UPLO
 146 *     ..
 147 *     .. Array Arguments ..
 148       REAL               A( LDA, * ), AFAC( LDAFAC, * ),
 149      $                   PERM( LDPERM, * ), RWORK( * )
 150       INTEGER            PIV( * )
 151 *     ..
 152 *
 153 *  =====================================================================
 154 *
 155 *     .. Parameters ..
 156       REAL               ZERO, ONE
 157       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
 158 *     ..
 159 *     .. Local Scalars ..
 160       REAL               ANORM, EPS, T
 161       INTEGER            I, J, K
 162 *     ..
 163 *     .. External Functions ..
 164       REAL               SDOT, SLAMCH, SLANSY
 165       LOGICAL            LSAME
 166       EXTERNAL           SDOT, SLAMCH, SLANSY, LSAME
 167 *     ..
 168 *     .. External Subroutines ..
 169       EXTERNAL           SSCAL, SSYR, STRMV
 170 *     ..
 171 *     .. Intrinsic Functions ..
 172       INTRINSIC          REAL
 173 *     ..
 174 *     .. Executable Statements ..
 175 *
 176 *     Quick exit if N = 0.
 177 *
 178       IF( N.LE.0 ) THEN
 179          RESID = ZERO
 180          RETURN
 181       END IF
 182 *
 183 *     Exit with RESID = 1/EPS if ANORM = 0.
 184 *
 185       EPS = SLAMCH( 'Epsilon' )
 186       ANORM = SLANSY( '1', UPLO, N, A, LDA, RWORK )
 187       IF( ANORM.LE.ZERO ) THEN
 188          RESID = ONE / EPS
 189          RETURN
 190       END IF
 191 *
 192 *     Compute the product U'*U, overwriting U.
 193 *
 194       IF( LSAME( UPLO, 'U' ) ) THEN
 195 *
 196          IF( RANK.LT.N ) THEN
 197             DO 110 J = RANK + 1, N
 198                DO 100 I = RANK + 1, J
 199                   AFAC( I, J ) = ZERO
 200   100          CONTINUE
 201   110       CONTINUE
 202          END IF
 203 *
 204          DO 120 K = N, 1, -1
 205 *
 206 *           Compute the (K,K) element of the result.
 207 *
 208             T = SDOT( K, AFAC( 1, K ), 1, AFAC( 1, K ), 1 )
 209             AFAC( K, K ) = T
 210 *
 211 *           Compute the rest of column K.
 212 *
 213             CALL STRMV( 'Upper', 'Transpose', 'Non-unit', K-1, AFAC,
 214      $                  LDAFAC, AFAC( 1, K ), 1 )
 215 *
 216   120    CONTINUE
 217 *
 218 *     Compute the product L*L', overwriting L.
 219 *
 220       ELSE
 221 *
 222          IF( RANK.LT.N ) THEN
 223             DO 140 J = RANK + 1, N
 224                DO 130 I = J, N
 225                   AFAC( I, J ) = ZERO
 226   130          CONTINUE
 227   140       CONTINUE
 228          END IF
 229 *
 230          DO 150 K = N, 1, -1
 231 *           Add a multiple of column K of the factor L to each of
 232 *           columns K+1 through N.
 233 *
 234             IF( K+1.LE.N )
 235      $         CALL SSYR( 'Lower', N-K, ONE, AFAC( K+1, K ), 1,
 236      $                    AFAC( K+1, K+1 ), LDAFAC )
 237 *
 238 *           Scale column K by the diagonal element.
 239 *
 240             T = AFAC( K, K )
 241             CALL SSCAL( N-K+1, T, AFAC( K, K ), 1 )
 242   150    CONTINUE
 243 *
 244       END IF
 245 *
 246 *        Form P*L*L'*P' or P*U'*U*P'
 247 *
 248       IF( LSAME( UPLO, 'U' ) ) THEN
 249 *
 250          DO 170 J = 1, N
 251             DO 160 I = 1, N
 252                IF( PIV( I ).LE.PIV( J ) ) THEN
 253                   IF( I.LE.J ) THEN
 254                      PERM( PIV( I ), PIV( J ) ) = AFAC( I, J )
 255                   ELSE
 256                      PERM( PIV( I ), PIV( J ) ) = AFAC( J, I )
 257                   END IF
 258                END IF
 259   160       CONTINUE
 260   170    CONTINUE
 261 *
 262 *
 263       ELSE
 264 *
 265          DO 190 J = 1, N
 266             DO 180 I = 1, N
 267                IF( PIV( I ).GE.PIV( J ) ) THEN
 268                   IF( I.GE.J ) THEN
 269                      PERM( PIV( I ), PIV( J ) ) = AFAC( I, J )
 270                   ELSE
 271                      PERM( PIV( I ), PIV( J ) ) = AFAC( J, I )
 272                   END IF
 273                END IF
 274   180       CONTINUE
 275   190    CONTINUE
 276 *
 277       END IF
 278 *
 279 *     Compute the difference  P*L*L'*P' - A (or P*U'*U*P' - A).
 280 *
 281       IF( LSAME( UPLO, 'U' ) ) THEN
 282          DO 210 J = 1, N
 283             DO 200 I = 1, J
 284                PERM( I, J ) = PERM( I, J ) - A( I, J )
 285   200       CONTINUE
 286   210    CONTINUE
 287       ELSE
 288          DO 230 J = 1, N
 289             DO 220 I = J, N
 290                PERM( I, J ) = PERM( I, J ) - A( I, J )
 291   220       CONTINUE
 292   230    CONTINUE
 293       END IF
 294 *
 295 *     Compute norm( P*L*L'P - A ) / ( N * norm(A) * EPS ), or
 296 *     ( P*U'*U*P' - A )/ ( N * norm(A) * EPS ).
 297 *
 298       RESID = SLANSY( '1', UPLO, N, PERM, LDAFAC, RWORK )
 299 *
 300       RESID = ( ( RESID / REAL( N ) ) / ANORM ) / EPS
 301 *
 302       RETURN
 303 *
 304 *     End of SPST01
 305 *
 306       END