1*> \brief \b DGET03
2*
3*  =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6*            http://www.netlib.org/lapack/explore-html/
7*
8*  Definition:
9*  ===========
10*
11*       SUBROUTINE DGET03( N, A, LDA, AINV, LDAINV, WORK, LDWORK, RWORK,
12*                          RCOND, RESID )
13*
14*       .. Scalar Arguments ..
15*       INTEGER            LDA, LDAINV, LDWORK, N
16*       DOUBLE PRECISION   RCOND, RESID
17*       ..
18*       .. Array Arguments ..
19*       DOUBLE PRECISION   A( LDA, * ), AINV( LDAINV, * ), RWORK( * ),
20*      $                   WORK( LDWORK, * )
21*       ..
22*
23*
24*> \par Purpose:
25*  =============
26*>
27*> \verbatim
28*>
29*> DGET03 computes the residual for a general matrix times its inverse:
30*>    norm( I - AINV*A ) / ( N * norm(A) * norm(AINV) * EPS ),
31*> where EPS is the machine epsilon.
32*> \endverbatim
33*
34*  Arguments:
35*  ==========
36*
37*> \param[in] N
38*> \verbatim
39*>          N is INTEGER
40*>          The number of rows and columns of the matrix A.  N >= 0.
41*> \endverbatim
42*>
43*> \param[in] A
44*> \verbatim
45*>          A is DOUBLE PRECISION array, dimension (LDA,N)
46*>          The original N x N matrix A.
47*> \endverbatim
48*>
49*> \param[in] LDA
50*> \verbatim
51*>          LDA is INTEGER
52*>          The leading dimension of the array A.  LDA >= max(1,N).
53*> \endverbatim
54*>
55*> \param[in] AINV
56*> \verbatim
57*>          AINV is DOUBLE PRECISION array, dimension (LDAINV,N)
58*>          The inverse of the matrix A.
59*> \endverbatim
60*>
61*> \param[in] LDAINV
62*> \verbatim
63*>          LDAINV is INTEGER
64*>          The leading dimension of the array AINV.  LDAINV >= max(1,N).
65*> \endverbatim
66*>
67*> \param[out] WORK
68*> \verbatim
69*>          WORK is DOUBLE PRECISION array, dimension (LDWORK,N)
70*> \endverbatim
71*>
72*> \param[in] LDWORK
73*> \verbatim
74*>          LDWORK is INTEGER
75*>          The leading dimension of the array WORK.  LDWORK >= max(1,N).
76*> \endverbatim
77*>
78*> \param[out] RWORK
79*> \verbatim
80*>          RWORK is DOUBLE PRECISION array, dimension (N)
81*> \endverbatim
82*>
83*> \param[out] RCOND
84*> \verbatim
85*>          RCOND is DOUBLE PRECISION
86*>          The reciprocal of the condition number of A, computed as
87*>          ( 1/norm(A) ) / norm(AINV).
88*> \endverbatim
89*>
90*> \param[out] RESID
91*> \verbatim
92*>          RESID is DOUBLE PRECISION
93*>          norm(I - AINV*A) / ( N * norm(A) * norm(AINV) * EPS )
94*> \endverbatim
95*
96*  Authors:
97*  ========
98*
99*> \author Univ. of Tennessee
100*> \author Univ. of California Berkeley
101*> \author Univ. of Colorado Denver
102*> \author NAG Ltd.
103*
104*> \ingroup double_lin
105*
106*  =====================================================================
107      SUBROUTINE DGET03( N, A, LDA, AINV, LDAINV, WORK, LDWORK, RWORK,
108     $                   RCOND, RESID )
109*
110*  -- LAPACK test routine --
111*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
112*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
113*
114*     .. Scalar Arguments ..
115      INTEGER            LDA, LDAINV, LDWORK, N
116      DOUBLE PRECISION   RCOND, RESID
117*     ..
118*     .. Array Arguments ..
119      DOUBLE PRECISION   A( LDA, * ), AINV( LDAINV, * ), RWORK( * ),
120     $                   WORK( LDWORK, * )
121*     ..
122*
123*  =====================================================================
124*
125*     .. Parameters ..
126      DOUBLE PRECISION   ZERO, ONE
127      PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
128*     ..
129*     .. Local Scalars ..
130      INTEGER            I
131      DOUBLE PRECISION   AINVNM, ANORM, EPS
132*     ..
133*     .. External Functions ..
134      DOUBLE PRECISION   DLAMCH, DLANGE
135      EXTERNAL           DLAMCH, DLANGE
136*     ..
137*     .. External Subroutines ..
138      EXTERNAL           DGEMM
139*     ..
140*     .. Intrinsic Functions ..
141      INTRINSIC          DBLE
142*     ..
143*     .. Executable Statements ..
144*
145*     Quick exit if N = 0.
146*
147      IF( N.LE.0 ) THEN
148         RCOND = ONE
149         RESID = ZERO
150         RETURN
151      END IF
152*
153*     Exit with RESID = 1/EPS if ANORM = 0 or AINVNM = 0.
154*
155      EPS = DLAMCH( 'Epsilon' )
156      ANORM = DLANGE( '1', N, N, A, LDA, RWORK )
157      AINVNM = DLANGE( '1', N, N, AINV, LDAINV, RWORK )
158      IF( ANORM.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN
159         RCOND = ZERO
160         RESID = ONE / EPS
161         RETURN
162      END IF
163      RCOND = ( ONE / ANORM ) / AINVNM
164*
165*     Compute I - A * AINV
166*
167      CALL DGEMM( 'No transpose', 'No transpose', N, N, N, -ONE, AINV,
168     $            LDAINV, A, LDA, ZERO, WORK, LDWORK )
169      DO 10 I = 1, N
170         WORK( I, I ) = ONE + WORK( I, I )
171   10 CONTINUE
172*
173*     Compute norm(I - AINV*A) / (N * norm(A) * norm(AINV) * EPS)
174*
175      RESID = DLANGE( '1', N, N, WORK, LDWORK, RWORK )
176*
177      RESID = ( ( RESID*RCOND ) / EPS ) / DBLE( N )
178*
179      RETURN
180*
181*     End of DGET03
182*
183      END
184