d4/de1/dget03_8f_source.html

00001       SUBROUTINE DGET03( N, A, LDA, AINV, LDAINV, WORK, LDWORK, RWORK,
00002      $                   RCOND, 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       INTEGER            LDA, LDAINV, LDWORK, N
00010       DOUBLE PRECISION   RCOND, RESID
00011 *     ..
00012 *     .. Array Arguments ..
00013       DOUBLE PRECISION   A( LDA, * ), AINV( LDAINV, * ), RWORK( * ),
00014      $                   WORK( LDWORK, * )
00015 *     ..
00016 *
00017 *  Purpose
00018 *  =======
00019 *
00020 *  DGET03 computes the residual for a general matrix times its inverse:
00021 *     norm( I - AINV*A ) / ( N * norm(A) * norm(AINV) * EPS ),
00022 *  where EPS is the machine epsilon.
00023 *
00024 *  Arguments
00025 *  ==========
00026 *
00027 *  N       (input) INTEGER
00028 *          The number of rows and columns of the matrix A.  N >= 0.
00029 *
00030 *  A       (input) DOUBLE PRECISION array, dimension (LDA,N)
00031 *          The original N x N matrix A.
00032 *
00033 *  LDA     (input) INTEGER
00034 *          The leading dimension of the array A.  LDA >= max(1,N).
00035 *
00036 *  AINV    (input) DOUBLE PRECISION array, dimension (LDAINV,N)
00037 *          The inverse of the matrix A.
00038 *
00039 *  LDAINV  (input) INTEGER
00040 *          The leading dimension of the array AINV.  LDAINV >= max(1,N).
00041 *
00042 *  WORK    (workspace) DOUBLE PRECISION array, dimension (LDWORK,N)
00043 *
00044 *  LDWORK  (input) INTEGER
00045 *          The leading dimension of the array WORK.  LDWORK >= max(1,N).
00046 *
00047 *  RWORK   (workspace) DOUBLE PRECISION array, dimension (N)
00048 *
00049 *  RCOND   (output) DOUBLE PRECISION
00050 *          The reciprocal of the condition number of A, computed as
00051 *          ( 1/norm(A) ) / norm(AINV).
00052 *
00053 *  RESID   (output) DOUBLE PRECISION
00054 *          norm(I - AINV*A) / ( N * norm(A) * norm(AINV) * EPS )
00055 *
00056 *  =====================================================================
00057 *
00058 *     .. Parameters ..
00059       DOUBLE PRECISION   ZERO, ONE
00060       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
00061 *     ..
00062 *     .. Local Scalars ..
00063       INTEGER            I
00064       DOUBLE PRECISION   AINVNM, ANORM, EPS
00065 *     ..
00066 *     .. External Functions ..
00067       DOUBLE PRECISION   DLAMCH, DLANGE
00068       EXTERNAL           DLAMCH, DLANGE
00069 *     ..
00070 *     .. External Subroutines ..
00071       EXTERNAL           DGEMM
00072 *     ..
00073 *     .. Intrinsic Functions ..
00074       INTRINSIC          DBLE
00075 *     ..
00076 *     .. Executable Statements ..
00077 *
00078 *     Quick exit if N = 0.
00079 *
00080       IF( N.LE.0 ) THEN
00081          RCOND = ONE
00082          RESID = ZERO
00083          RETURN
00084       END IF
00085 *
00086 *     Exit with RESID = 1/EPS if ANORM = 0 or AINVNM = 0.
00087 *
00088       EPS = DLAMCH( 'Epsilon' )
00089       ANORM = DLANGE( '1', N, N, A, LDA, RWORK )
00090       AINVNM = DLANGE( '1', N, N, AINV, LDAINV, RWORK )
00091       IF( ANORM.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN
00092          RCOND = ZERO
00093          RESID = ONE / EPS
00094          RETURN
00095       END IF
00096       RCOND = ( ONE / ANORM ) / AINVNM
00097 *
00098 *     Compute I - A * AINV
00099 *
00100       CALL DGEMM( 'No transpose', 'No transpose', N, N, N, -ONE, AINV,
00101      $            LDAINV, A, LDA, ZERO, WORK, LDWORK )
00102       DO 10 I = 1, N
00103          WORK( I, I ) = ONE + WORK( I, I )
00104    10 CONTINUE
00105 *
00106 *     Compute norm(I - AINV*A) / (N * norm(A) * norm(AINV) * EPS)
00107 *
00108       RESID = DLANGE( '1', N, N, WORK, LDWORK, RWORK )
00109 *
00110       RESID = ( ( RESID*RCOND ) / EPS ) / DBLE( N )
00111 *
00112       RETURN
00113 *
00114 *     End of DGET03
00115 *
00116       END