LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
ssyt01.f
Go to the documentation of this file.
00001 *> \brief \b SSYT01
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 SSYT01( UPLO, N, A, LDA, AFAC, LDAFAC, IPIV, C, LDC,
00012 *                          RWORK, RESID )
00013 * 
00014 *       .. Scalar Arguments ..
00015 *       CHARACTER          UPLO
00016 *       INTEGER            LDA, LDAFAC, LDC, N
00017 *       REAL               RESID
00018 *       ..
00019 *       .. Array Arguments ..
00020 *       INTEGER            IPIV( * )
00021 *       REAL               A( LDA, * ), AFAC( LDAFAC, * ), C( LDC, * ),
00022 *      $                   RWORK( * )
00023 *       ..
00024 *  
00025 *
00026 *> \par Purpose:
00027 *  =============
00028 *>
00029 *> \verbatim
00030 *>
00031 *> SSYT01 reconstructs a symmetric indefinite matrix A from its
00032 *> block L*D*L' or U*D*U' factorization and computes the residual
00033 *>    norm( C - A ) / ( N * norm(A) * EPS ),
00034 *> where C is the reconstructed matrix and EPS is the machine epsilon.
00035 *> \endverbatim
00036 *
00037 *  Arguments:
00038 *  ==========
00039 *
00040 *> \param[in] UPLO
00041 *> \verbatim
00042 *>          UPLO is CHARACTER*1
00043 *>          Specifies whether the upper or lower triangular part of the
00044 *>          symmetric matrix A is stored:
00045 *>          = 'U':  Upper triangular
00046 *>          = 'L':  Lower triangular
00047 *> \endverbatim
00048 *>
00049 *> \param[in] N
00050 *> \verbatim
00051 *>          N is INTEGER
00052 *>          The number of rows and columns of the matrix A.  N >= 0.
00053 *> \endverbatim
00054 *>
00055 *> \param[in] A
00056 *> \verbatim
00057 *>          A is REAL array, dimension (LDA,N)
00058 *>          The original symmetric matrix A.
00059 *> \endverbatim
00060 *>
00061 *> \param[in] LDA
00062 *> \verbatim
00063 *>          LDA is INTEGER
00064 *>          The leading dimension of the array A.  LDA >= max(1,N)
00065 *> \endverbatim
00066 *>
00067 *> \param[in] AFAC
00068 *> \verbatim
00069 *>          AFAC is REAL array, dimension (LDAFAC,N)
00070 *>          The factored form of the matrix A.  AFAC contains the block
00071 *>          diagonal matrix D and the multipliers used to obtain the
00072 *>          factor L or U from the block L*D*L' or U*D*U' factorization
00073 *>          as computed by SSYTRF.
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[in] IPIV
00083 *> \verbatim
00084 *>          IPIV is INTEGER array, dimension (N)
00085 *>          The pivot indices from SSYTRF.
00086 *> \endverbatim
00087 *>
00088 *> \param[out] C
00089 *> \verbatim
00090 *>          C is REAL array, dimension (LDC,N)
00091 *> \endverbatim
00092 *>
00093 *> \param[in] LDC
00094 *> \verbatim
00095 *>          LDC is INTEGER
00096 *>          The leading dimension of the array C.  LDC >= max(1,N).
00097 *> \endverbatim
00098 *>
00099 *> \param[out] RWORK
00100 *> \verbatim
00101 *>          RWORK is REAL array, dimension (N)
00102 *> \endverbatim
00103 *>
00104 *> \param[out] RESID
00105 *> \verbatim
00106 *>          RESID is REAL
00107 *>          If UPLO = 'L', norm(L*D*L' - A) / ( N * norm(A) * EPS )
00108 *>          If UPLO = 'U', norm(U*D*U' - A) / ( N * norm(A) * EPS )
00109 *> \endverbatim
00110 *
00111 *  Authors:
00112 *  ========
00113 *
00114 *> \author Univ. of Tennessee 
00115 *> \author Univ. of California Berkeley 
00116 *> \author Univ. of Colorado Denver 
00117 *> \author NAG Ltd. 
00118 *
00119 *> \date April 2012
00120 *
00121 *> \ingroup single_lin
00122 *
00123 *  =====================================================================
00124       SUBROUTINE SSYT01( UPLO, N, A, LDA, AFAC, LDAFAC, IPIV, C, LDC,
00125      $                   RWORK, RESID )
00126 *
00127 *  -- LAPACK test routine (version 3.4.1) --
00128 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00129 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00130 *     April 2012
00131 *
00132 *     .. Scalar Arguments ..
00133       CHARACTER          UPLO
00134       INTEGER            LDA, LDAFAC, LDC, N
00135       REAL               RESID
00136 *     ..
00137 *     .. Array Arguments ..
00138       INTEGER            IPIV( * )
00139       REAL               A( LDA, * ), AFAC( LDAFAC, * ), C( LDC, * ),
00140      $                   RWORK( * )
00141 *     ..
00142 *
00143 *  =====================================================================
00144 *
00145 *     .. Parameters ..
00146       REAL               ZERO, ONE
00147       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
00148 *     ..
00149 *     .. Local Scalars ..
00150       INTEGER            I, INFO, J
00151       REAL               ANORM, EPS
00152 *     ..
00153 *     .. External Functions ..
00154       LOGICAL            LSAME
00155       REAL               SLAMCH, SLANSY
00156       EXTERNAL           LSAME, SLAMCH, SLANSY
00157 *     ..
00158 *     .. External Subroutines ..
00159       EXTERNAL           SLASET, SLAVSY
00160 *     ..
00161 *     .. Intrinsic Functions ..
00162       INTRINSIC          REAL
00163 *     ..
00164 *     .. Executable Statements ..
00165 *
00166 *     Quick exit if N = 0.
00167 *
00168       IF( N.LE.0 ) THEN
00169          RESID = ZERO
00170          RETURN
00171       END IF
00172 *
00173 *     Determine EPS and the norm of A.
00174 *
00175       EPS = SLAMCH( 'Epsilon' )
00176       ANORM = SLANSY( '1', UPLO, N, A, LDA, RWORK )
00177 *
00178 *     Initialize C to the identity matrix.
00179 *
00180       CALL SLASET( 'Full', N, N, ZERO, ONE, C, LDC )
00181 *
00182 *     Call SLAVSY to form the product D * U' (or D * L' ).
00183 *
00184       CALL SLAVSY( UPLO, 'Transpose', 'Non-unit', N, N, AFAC, LDAFAC,
00185      $             IPIV, C, LDC, INFO )
00186 *
00187 *     Call SLAVSY again to multiply by U (or L ).
00188 *
00189       CALL SLAVSY( UPLO, 'No transpose', 'Unit', N, N, AFAC, LDAFAC,
00190      $             IPIV, C, LDC, INFO )
00191 *
00192 *     Compute the difference  C - A .
00193 *
00194       IF( LSAME( UPLO, 'U' ) ) THEN
00195          DO 20 J = 1, N
00196             DO 10 I = 1, J
00197                C( I, J ) = C( I, J ) - A( I, J )
00198    10       CONTINUE
00199    20    CONTINUE
00200       ELSE
00201          DO 40 J = 1, N
00202             DO 30 I = J, N
00203                C( I, J ) = C( I, J ) - A( I, J )
00204    30       CONTINUE
00205    40    CONTINUE
00206       END IF
00207 *
00208 *     Compute norm( C - A ) / ( N * norm(A) * EPS )
00209 *
00210       RESID = SLANSY( '1', UPLO, N, C, LDC, RWORK )
00211 *
00212       IF( ANORM.LE.ZERO ) THEN
00213          IF( RESID.NE.ZERO )
00214      $      RESID = ONE / EPS
00215       ELSE
00216          RESID = ( ( RESID / REAL( N ) ) / ANORM ) / EPS
00217       END IF
00218 *
00219       RETURN
00220 *
00221 *     End of SSYT01
00222 *
00223       END
 All Files Functions