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