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