LAPACK: TESTING/LIN/dppt03.f Source File

Go to the documentation of this file.
00001       SUBROUTINE DPPT03( UPLO, N, A, AINV, WORK, LDWORK, RWORK, RCOND,
00002      $                   RESID )
00003 *
00004 *  -- LAPACK test routine (version 3.1) --
00005 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
00006 *     November 2006
00007 *
00008 *     .. Scalar Arguments ..
00009       CHARACTER          UPLO
00010       INTEGER            LDWORK, N
00011       DOUBLE PRECISION   RCOND, RESID
00012 *     ..
00013 *     .. Array Arguments ..
00014       DOUBLE PRECISION   A( * ), AINV( * ), RWORK( * ),
00015      $                   WORK( LDWORK, * )
00016 *     ..
00017 *
00018 *  Purpose
00019 *  =======
00020 *
00021 *  DPPT03 computes the residual for a symmetric packed matrix times its
00022 *  inverse:
00023 *     norm( I - A*AINV ) / ( N * norm(A) * norm(AINV) * EPS ),
00024 *  where EPS is the machine epsilon.
00025 *
00026 *  Arguments
00027 *  ==========
00028 *
00029 *  UPLO    (input) CHARACTER*1
00030 *          Specifies whether the upper or lower triangular part of the
00031 *          symmetric matrix A is stored:
00032 *          = 'U':  Upper triangular
00033 *          = 'L':  Lower triangular
00034 *
00035 *  N       (input) INTEGER
00036 *          The number of rows and columns of the matrix A.  N >= 0.
00037 *
00038 *  A       (input) DOUBLE PRECISION array, dimension (N*(N+1)/2)
00039 *          The original symmetric matrix A, stored as a packed
00040 *          triangular matrix.
00041 *
00042 *  AINV    (input) DOUBLE PRECISION array, dimension (N*(N+1)/2)
00043 *          The (symmetric) inverse of the matrix A, stored as a packed
00044 *          triangular matrix.
00045 *
00046 *  WORK    (workspace) DOUBLE PRECISION array, dimension (LDWORK,N)
00047 *
00048 *  LDWORK  (input) INTEGER
00049 *          The leading dimension of the array WORK.  LDWORK >= max(1,N).
00050 *
00051 *  RWORK   (workspace) DOUBLE PRECISION array, dimension (N)
00052 *
00053 *  RCOND   (output) DOUBLE PRECISION
00054 *          The reciprocal of the condition number of A, computed as
00055 *          ( 1/norm(A) ) / norm(AINV).
00056 *
00057 *  RESID   (output) DOUBLE PRECISION
00058 *          norm(I - A*AINV) / ( N * norm(A) * norm(AINV) * EPS )
00059 *
00060 *  =====================================================================
00061 *
00062 *     .. Parameters ..
00063       DOUBLE PRECISION   ZERO, ONE
00064       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
00065 *     ..
00066 *     .. Local Scalars ..
00067       INTEGER            I, J, JJ
00068       DOUBLE PRECISION   AINVNM, ANORM, EPS
00069 *     ..
00070 *     .. External Functions ..
00071       LOGICAL            LSAME
00072       DOUBLE PRECISION   DLAMCH, DLANGE, DLANSP
00073       EXTERNAL           LSAME, DLAMCH, DLANGE, DLANSP
00074 *     ..
00075 *     .. Intrinsic Functions ..
00076       INTRINSIC          DBLE
00077 *     ..
00078 *     .. External Subroutines ..
00079       EXTERNAL           DCOPY, DSPMV
00080 *     ..
00081 *     .. Executable Statements ..
00082 *
00083 *     Quick exit if N = 0.
00084 *
00085       IF( N.LE.0 ) THEN
00086          RCOND = ONE
00087          RESID = ZERO
00088          RETURN
00089       END IF
00090 *
00091 *     Exit with RESID = 1/EPS if ANORM = 0 or AINVNM = 0.
00092 *
00093       EPS = DLAMCH( 'Epsilon' )
00094       ANORM = DLANSP( '1', UPLO, N, A, RWORK )
00095       AINVNM = DLANSP( '1', UPLO, N, AINV, RWORK )
00096       IF( ANORM.LE.ZERO .OR. AINVNM.EQ.ZERO ) THEN
00097          RCOND = ZERO
00098          RESID = ONE / EPS
00099          RETURN
00100       END IF
00101       RCOND = ( ONE / ANORM ) / AINVNM
00102 *
00103 *     UPLO = 'U':
00104 *     Copy the leading N-1 x N-1 submatrix of AINV to WORK(1:N,2:N) and
00105 *     expand it to a full matrix, then multiply by A one column at a
00106 *     time, moving the result one column to the left.
00107 *
00108       IF( LSAME( UPLO, 'U' ) ) THEN
00109 *
00110 *        Copy AINV
00111 *
00112          JJ = 1
00113          DO 10 J = 1, N - 1
00114             CALL DCOPY( J, AINV( JJ ), 1, WORK( 1, J+1 ), 1 )
00115             CALL DCOPY( J-1, AINV( JJ ), 1, WORK( J, 2 ), LDWORK )
00116             JJ = JJ + J
00117    10    CONTINUE
00118          JJ = ( ( N-1 )*N ) / 2 + 1
00119          CALL DCOPY( N-1, AINV( JJ ), 1, WORK( N, 2 ), LDWORK )
00120 *
00121 *        Multiply by A
00122 *
00123          DO 20 J = 1, N - 1
00124             CALL DSPMV( 'Upper', N, -ONE, A, WORK( 1, J+1 ), 1, ZERO,
00125      $                  WORK( 1, J ), 1 )
00126    20    CONTINUE
00127          CALL DSPMV( 'Upper', N, -ONE, A, AINV( JJ ), 1, ZERO,
00128      $               WORK( 1, N ), 1 )
00129 *
00130 *     UPLO = 'L':
00131 *     Copy the trailing N-1 x N-1 submatrix of AINV to WORK(1:N,1:N-1)
00132 *     and multiply by A, moving each column to the right.
00133 *
00134       ELSE
00135 *
00136 *        Copy AINV
00137 *
00138          CALL DCOPY( N-1, AINV( 2 ), 1, WORK( 1, 1 ), LDWORK )
00139          JJ = N + 1
00140          DO 30 J = 2, N
00141             CALL DCOPY( N-J+1, AINV( JJ ), 1, WORK( J, J-1 ), 1 )
00142             CALL DCOPY( N-J, AINV( JJ+1 ), 1, WORK( J, J ), LDWORK )
00143             JJ = JJ + N - J + 1
00144    30    CONTINUE
00145 *
00146 *        Multiply by A
00147 *
00148          DO 40 J = N, 2, -1
00149             CALL DSPMV( 'Lower', N, -ONE, A, WORK( 1, J-1 ), 1, ZERO,
00150      $                  WORK( 1, J ), 1 )
00151    40    CONTINUE
00152          CALL DSPMV( 'Lower', N, -ONE, A, AINV( 1 ), 1, ZERO,
00153      $               WORK( 1, 1 ), 1 )
00154 *
00155       END IF
00156 *
00157 *     Add the identity matrix to WORK .
00158 *
00159       DO 50 I = 1, N
00160          WORK( I, I ) = WORK( I, I ) + ONE
00161    50 CONTINUE
00162 *
00163 *     Compute norm(I - A*AINV) / (N * norm(A) * norm(AINV) * EPS)
00164 *
00165       RESID = DLANGE( '1', N, N, WORK, LDWORK, RWORK )
00166 *
00167       RESID = ( ( RESID*RCOND ) / EPS ) / DBLE( N )
00168 *
00169       RETURN
00170 *
00171 *     End of DPPT03
00172 *
00173       END