LAPACK 3.3.0

sspt01.f

Go to the documentation of this file.
00001       SUBROUTINE SSPT01( UPLO, N, A, AFAC, IPIV, C, LDC, RWORK, RESID )
00002 *
00003 *  -- LAPACK test routine (version 3.1) --
00004 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
00005 *     November 2006
00006 *
00007 *     .. Scalar Arguments ..
00008       CHARACTER          UPLO
00009       INTEGER            LDC, N
00010       REAL               RESID
00011 *     ..
00012 *     .. Array Arguments ..
00013       INTEGER            IPIV( * )
00014       REAL               A( * ), AFAC( * ), C( LDC, * ), RWORK( * )
00015 *     ..
00016 *
00017 *  Purpose
00018 *  =======
00019 *
00020 *  SSPT01 reconstructs a symmetric indefinite packed matrix A from its
00021 *  block L*D*L' or U*D*U' factorization and computes the residual
00022 *       norm( C - A ) / ( N * norm(A) * EPS ),
00023 *  where C is the reconstructed matrix and EPS is the machine epsilon.
00024 *
00025 *  Arguments
00026 *  ==========
00027 *
00028 *  UPLO    (input) CHARACTER*1
00029 *          Specifies whether the upper or lower triangular part of the
00030 *          symmetric matrix A is stored:
00031 *          = 'U':  Upper triangular
00032 *          = 'L':  Lower triangular
00033 *
00034 *  N       (input) INTEGER
00035 *          The number of rows and columns of the matrix A.  N >= 0.
00036 *
00037 *  A       (input) REAL array, dimension (N*(N+1)/2)
00038 *          The original symmetric matrix A, stored as a packed
00039 *          triangular matrix.
00040 *
00041 *  AFAC    (input) REAL array, dimension (N*(N+1)/2)
00042 *          The factored form of the matrix A, stored as a packed
00043 *          triangular matrix.  AFAC contains the block diagonal matrix D
00044 *          and the multipliers used to obtain the factor L or U from the
00045 *          block L*D*L' or U*D*U' factorization as computed by SSPTRF.
00046 *
00047 *  IPIV    (input) INTEGER array, dimension (N)
00048 *          The pivot indices from SSPTRF.
00049 *
00050 *  C       (workspace) REAL array, dimension (LDC,N)
00051 *
00052 *  LDC     (integer) INTEGER
00053 *          The leading dimension of the array C.  LDC >= max(1,N).
00054 *
00055 *  RWORK   (workspace) REAL array, dimension (N)
00056 *
00057 *  RESID   (output) REAL
00058 *          If UPLO = 'L', norm(L*D*L' - A) / ( N * norm(A) * EPS )
00059 *          If UPLO = 'U', norm(U*D*U' - A) / ( N * norm(A) * EPS )
00060 *
00061 *  =====================================================================
00062 *
00063 *     .. Parameters ..
00064       REAL               ZERO, ONE
00065       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
00066 *     ..
00067 *     .. Local Scalars ..
00068       INTEGER            I, INFO, J, JC
00069       REAL               ANORM, EPS
00070 *     ..
00071 *     .. External Functions ..
00072       LOGICAL            LSAME
00073       REAL               SLAMCH, SLANSP, SLANSY
00074       EXTERNAL           LSAME, SLAMCH, SLANSP, SLANSY
00075 *     ..
00076 *     .. External Subroutines ..
00077       EXTERNAL           SLAVSP, SLASET
00078 *     ..
00079 *     .. Intrinsic Functions ..
00080       INTRINSIC          REAL
00081 *     ..
00082 *     .. Executable Statements ..
00083 *
00084 *     Quick exit if N = 0.
00085 *
00086       IF( N.LE.0 ) THEN
00087          RESID = ZERO
00088          RETURN
00089       END IF
00090 *
00091 *     Determine EPS and the norm of A.
00092 *
00093       EPS = SLAMCH( 'Epsilon' )
00094       ANORM = SLANSP( '1', UPLO, N, A, RWORK )
00095 *
00096 *     Initialize C to the identity matrix.
00097 *
00098       CALL SLASET( 'Full', N, N, ZERO, ONE, C, LDC )
00099 *
00100 *     Call SLAVSP to form the product D * U' (or D * L' ).
00101 *
00102       CALL SLAVSP( UPLO, 'Transpose', 'Non-unit', N, N, AFAC, IPIV, C,
00103      $             LDC, INFO )
00104 *
00105 *     Call SLAVSP again to multiply by U ( or L ).
00106 *
00107       CALL SLAVSP( UPLO, 'No transpose', 'Unit', N, N, AFAC, IPIV, C,
00108      $             LDC, INFO )
00109 *
00110 *     Compute the difference  C - A .
00111 *
00112       IF( LSAME( UPLO, 'U' ) ) THEN
00113          JC = 0
00114          DO 20 J = 1, N
00115             DO 10 I = 1, J
00116                C( I, J ) = C( I, J ) - A( JC+I )
00117    10       CONTINUE
00118             JC = JC + J
00119    20    CONTINUE
00120       ELSE
00121          JC = 1
00122          DO 40 J = 1, N
00123             DO 30 I = J, N
00124                C( I, J ) = C( I, J ) - A( JC+I-J )
00125    30       CONTINUE
00126             JC = JC + N - J + 1
00127    40    CONTINUE
00128       END IF
00129 *
00130 *     Compute norm( C - A ) / ( N * norm(A) * EPS )
00131 *
00132       RESID = SLANSY( '1', UPLO, N, C, LDC, RWORK )
00133 *
00134       IF( ANORM.LE.ZERO ) THEN
00135          IF( RESID.NE.ZERO )
00136      $      RESID = ONE / EPS
00137       ELSE
00138          RESID = ( ( RESID / REAL( N ) ) / ANORM ) / EPS
00139       END IF
00140 *
00141       RETURN
00142 *
00143 *     End of SSPT01
00144 *
00145       END
 All Files Functions