1*> \brief \b ZSYT03
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 ZSYT03( 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*> ZSYT03 computes the residual for a complex symmetric matrix times
32*> its 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*>          complex symmetric 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 complex symmetric 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 symmetric
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*>          RCOND = 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 December 2016
122*
123*> \ingroup complex16_lin
124*
125*  =====================================================================
126      SUBROUTINE ZSYT03( UPLO, N, A, LDA, AINV, LDAINV, WORK, LDWORK,
127     $                   RWORK, RCOND, RESID )
128*
129*  -- LAPACK test routine (version 3.7.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*     December 2016
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*
148*     .. Parameters ..
149      DOUBLE PRECISION   ZERO, ONE
150      PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
151      COMPLEX*16         CZERO, CONE
152      PARAMETER          ( CZERO = ( 0.0D+0, 0.0D+0 ),
153     $                   CONE = ( 1.0D+0, 0.0D+0 ) )
154*     ..
155*     .. Local Scalars ..
156      INTEGER            I, J
157      DOUBLE PRECISION   AINVNM, ANORM, EPS
158*     ..
159*     .. External Functions ..
160      LOGICAL            LSAME
161      DOUBLE PRECISION   DLAMCH, ZLANGE, ZLANSY
162      EXTERNAL           LSAME, DLAMCH, ZLANGE, ZLANSY
163*     ..
164*     .. External Subroutines ..
165      EXTERNAL           ZSYMM
166*     ..
167*     .. Intrinsic Functions ..
168      INTRINSIC          DBLE
169*     ..
170*     .. Executable Statements ..
171*
172*     Quick exit if N = 0
173*
174      IF( N.LE.0 ) THEN
175         RCOND = ONE
176         RESID = ZERO
177         RETURN
178      END IF
179*
180*     Exit with RESID = 1/EPS if ANORM = 0 or AINVNM = 0.
181*
182      EPS = DLAMCH( 'Epsilon' )
183      ANORM = ZLANSY( '1', UPLO, N, A, LDA, RWORK )
184      AINVNM = ZLANSY( '1', UPLO, N, AINV, LDAINV, RWORK )
185      IF( ANORM.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN
186         RCOND = ZERO
187         RESID = ONE / EPS
188         RETURN
189      END IF
190      RCOND = ( ONE / ANORM ) / AINVNM
191*
192*     Expand AINV into a full matrix and call ZSYMM to multiply
193*     AINV on the left by A (store the result in WORK).
194*
195      IF( LSAME( UPLO, 'U' ) ) THEN
196         DO 20 J = 1, N
197            DO 10 I = 1, J - 1
198               AINV( J, I ) = AINV( I, J )
199   10       CONTINUE
200   20    CONTINUE
201      ELSE
202         DO 40 J = 1, N
203            DO 30 I = J + 1, N
204               AINV( J, I ) = AINV( I, J )
205   30       CONTINUE
206   40    CONTINUE
207      END IF
208      CALL ZSYMM( 'Left', UPLO, N, N, -CONE, A, LDA, AINV, LDAINV,
209     $            CZERO, WORK, LDWORK )
210*
211*     Add the identity matrix to WORK .
212*
213      DO 50 I = 1, N
214         WORK( I, I ) = WORK( I, I ) + CONE
215   50 CONTINUE
216*
217*     Compute norm(I - A*AINV) / (N * norm(A) * norm(AINV) * EPS)
218*
219      RESID = ZLANGE( '1', N, N, WORK, LDWORK, RWORK )
220*
221      RESID = ( ( RESID*RCOND ) / EPS ) / DBLE( N )
222*
223      RETURN
224*
225*     End of ZSYT03
226*
227      END
228