1*> \brief \b ZPOT03
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 ZPOT03( UPLO, N, A, LDA, AINV, LDAINV, WORK, LDWORK,
12*                          RWORK, RCOND, RESID )
13*
14*       .. Scalar Arguments ..
15*       CHARACTER          UPLO
16*       INTEGER            LDA, LDAINV, LDWORK, N
17*       DOUBLE PRECISION   RCOND, RESID
18*       ..
19*       .. Array Arguments ..
20*       DOUBLE PRECISION   RWORK( * )
21*       COMPLEX*16         A( LDA, * ), AINV( LDAINV, * ),
22*      $                   WORK( LDWORK, * )
23*       ..
24*
25*
26*> \par Purpose:
27*  =============
28*>
29*> \verbatim
30*>
31*> ZPOT03 computes the residual for a Hermitian matrix times its
32*> inverse:
33*>    norm( I - A*AINV ) / ( N * norm(A) * norm(AINV) * EPS ),
34*> where EPS is the machine epsilon.
35*> \endverbatim
36*
37*  Arguments:
38*  ==========
39*
40*> \param[in] UPLO
41*> \verbatim
42*>          UPLO is CHARACTER*1
43*>          Specifies whether the upper or lower triangular part of the
44*>          Hermitian matrix A is stored:
45*>          = 'U':  Upper triangular
46*>          = 'L':  Lower triangular
47*> \endverbatim
48*>
49*> \param[in] N
50*> \verbatim
51*>          N is INTEGER
52*>          The number of rows and columns of the matrix A.  N >= 0.
53*> \endverbatim
54*>
55*> \param[in] A
56*> \verbatim
57*>          A is COMPLEX*16 array, dimension (LDA,N)
58*>          The original Hermitian matrix A.
59*> \endverbatim
60*>
61*> \param[in] LDA
62*> \verbatim
63*>          LDA is INTEGER
64*>          The leading dimension of the array A.  LDA >= max(1,N)
65*> \endverbatim
66*>
67*> \param[in,out] AINV
68*> \verbatim
69*>          AINV is COMPLEX*16 array, dimension (LDAINV,N)
70*>          On entry, the inverse of the matrix A, stored as a Hermitian
71*>          matrix in the same format as A.
72*>          In this version, AINV is expanded into a full matrix and
73*>          multiplied by A, so the opposing triangle of AINV will be
74*>          changed; i.e., if the upper triangular part of AINV is
75*>          stored, the lower triangular part will be used as work space.
76*> \endverbatim
77*>
78*> \param[in] LDAINV
79*> \verbatim
80*>          LDAINV is INTEGER
81*>          The leading dimension of the array AINV.  LDAINV >= max(1,N).
82*> \endverbatim
83*>
84*> \param[out] WORK
85*> \verbatim
86*>          WORK is COMPLEX*16 array, dimension (LDWORK,N)
87*> \endverbatim
88*>
89*> \param[in] LDWORK
90*> \verbatim
91*>          LDWORK is INTEGER
92*>          The leading dimension of the array WORK.  LDWORK >= max(1,N).
93*> \endverbatim
94*>
95*> \param[out] RWORK
96*> \verbatim
97*>          RWORK is DOUBLE PRECISION array, dimension (N)
98*> \endverbatim
99*>
100*> \param[out] RCOND
101*> \verbatim
102*>          RCOND is DOUBLE PRECISION
103*>          The reciprocal of the condition number of A, computed as
104*>          ( 1/norm(A) ) / norm(AINV).
105*> \endverbatim
106*>
107*> \param[out] RESID
108*> \verbatim
109*>          RESID is DOUBLE PRECISION
110*>          norm(I - A*AINV) / ( N * norm(A) * norm(AINV) * EPS )
111*> \endverbatim
112*
113*  Authors:
114*  ========
115*
116*> \author Univ. of Tennessee
117*> \author Univ. of California Berkeley
118*> \author Univ. of Colorado Denver
119*> \author NAG Ltd.
120*
121*> \date November 2011
122*
123*> \ingroup complex16_lin
124*
125*  =====================================================================
126      SUBROUTINE ZPOT03( UPLO, N, A, LDA, AINV, LDAINV, WORK, LDWORK,
127     $                   RWORK, RCOND, RESID )
128*
129*  -- LAPACK test routine (version 3.4.0) --
130*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
131*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
132*     November 2011
133*
134*     .. Scalar Arguments ..
135      CHARACTER          UPLO
136      INTEGER            LDA, LDAINV, LDWORK, N
137      DOUBLE PRECISION   RCOND, RESID
138*     ..
139*     .. Array Arguments ..
140      DOUBLE PRECISION   RWORK( * )
141      COMPLEX*16         A( LDA, * ), AINV( LDAINV, * ),
142     $                   WORK( LDWORK, * )
143*     ..
144*
145*  =====================================================================
146*
147*     .. Parameters ..
148      DOUBLE PRECISION   ZERO, ONE
149      PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
150      COMPLEX*16         CZERO, CONE
151      PARAMETER          ( CZERO = ( 0.0D+0, 0.0D+0 ),
152     $                   CONE = ( 1.0D+0, 0.0D+0 ) )
153*     ..
154*     .. Local Scalars ..
155      INTEGER            I, J
156      DOUBLE PRECISION   AINVNM, ANORM, EPS
157*     ..
158*     .. External Functions ..
159      LOGICAL            LSAME
160      DOUBLE PRECISION   DLAMCH, ZLANGE, ZLANHE
161      EXTERNAL           LSAME, DLAMCH, ZLANGE, ZLANHE
162*     ..
163*     .. External Subroutines ..
164      EXTERNAL           ZHEMM
165*     ..
166*     .. Intrinsic Functions ..
167      INTRINSIC          DBLE, DCONJG
168*     ..
169*     .. Executable Statements ..
170*
171*     Quick exit if N = 0.
172*
173      IF( N.LE.0 ) THEN
174         RCOND = ONE
175         RESID = ZERO
176         RETURN
177      END IF
178*
179*     Exit with RESID = 1/EPS if ANORM = 0 or AINVNM = 0.
180*
181      EPS = DLAMCH( 'Epsilon' )
182      ANORM = ZLANHE( '1', UPLO, N, A, LDA, RWORK )
183      AINVNM = ZLANHE( '1', UPLO, N, AINV, LDAINV, RWORK )
184      IF( ANORM.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN
185         RCOND = ZERO
186         RESID = ONE / EPS
187         RETURN
188      END IF
189      RCOND = ( ONE / ANORM ) / AINVNM
190*
191*     Expand AINV into a full matrix and call ZHEMM to multiply
192*     AINV on the left by A.
193*
194      IF( LSAME( UPLO, 'U' ) ) THEN
195         DO 20 J = 1, N
196            DO 10 I = 1, J - 1
197               AINV( J, I ) = DCONJG( AINV( I, J ) )
198   10       CONTINUE
199   20    CONTINUE
200      ELSE
201         DO 40 J = 1, N
202            DO 30 I = J + 1, N
203               AINV( J, I ) = DCONJG( AINV( I, J ) )
204   30       CONTINUE
205   40    CONTINUE
206      END IF
207      CALL ZHEMM( 'Left', UPLO, N, N, -CONE, A, LDA, AINV, LDAINV,
208     $            CZERO, WORK, LDWORK )
209*
210*     Add the identity matrix to WORK .
211*
212      DO 50 I = 1, N
213         WORK( I, I ) = WORK( I, I ) + CONE
214   50 CONTINUE
215*
216*     Compute norm(I - A*AINV) / (N * norm(A) * norm(AINV) * EPS)
217*
218      RESID = ZLANGE( '1', N, N, WORK, LDWORK, RWORK )
219*
220      RESID = ( ( RESID*RCOND ) / EPS ) / DBLE( N )
221*
222      RETURN
223*
224*     End of ZPOT03
225*
226      END
227