1*> \brief \b ZPPT05
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 ZPPT05( UPLO, N, NRHS, AP, B, LDB, X, LDX, XACT,
12*                          LDXACT, FERR, BERR, RESLTS )
13*
14*       .. Scalar Arguments ..
15*       CHARACTER          UPLO
16*       INTEGER            LDB, LDX, LDXACT, N, NRHS
17*       ..
18*       .. Array Arguments ..
19*       DOUBLE PRECISION   BERR( * ), FERR( * ), RESLTS( * )
20*       COMPLEX*16         AP( * ), B( LDB, * ), X( LDX, * ),
21*      $                   XACT( LDXACT, * )
22*       ..
23*
24*
25*> \par Purpose:
26*  =============
27*>
28*> \verbatim
29*>
30*> ZPPT05 tests the error bounds from iterative refinement for the
31*> computed solution to a system of equations A*X = B, where A is a
32*> Hermitian matrix in packed storage format.
33*>
34*> RESLTS(1) = test of the error bound
35*>           = norm(X - XACT) / ( norm(X) * FERR )
36*>
37*> A large value is returned if this ratio is not less than one.
38*>
39*> RESLTS(2) = residual from the iterative refinement routine
40*>           = the maximum of BERR / ( (n+1)*EPS + (*) ), where
41*>             (*) = (n+1)*UNFL / (min_i (abs(A)*abs(X) +abs(b))_i )
42*> \endverbatim
43*
44*  Arguments:
45*  ==========
46*
47*> \param[in] UPLO
48*> \verbatim
49*>          UPLO is CHARACTER*1
50*>          Specifies whether the upper or lower triangular part of the
51*>          Hermitian matrix A is stored.
52*>          = 'U':  Upper triangular
53*>          = 'L':  Lower triangular
54*> \endverbatim
55*>
56*> \param[in] N
57*> \verbatim
58*>          N is INTEGER
59*>          The number of rows of the matrices X, B, and XACT, and the
60*>          order of the matrix A.  N >= 0.
61*> \endverbatim
62*>
63*> \param[in] NRHS
64*> \verbatim
65*>          NRHS is INTEGER
66*>          The number of columns of the matrices X, B, and XACT.
67*>          NRHS >= 0.
68*> \endverbatim
69*>
70*> \param[in] AP
71*> \verbatim
72*>          AP is COMPLEX*16 array, dimension (N*(N+1)/2)
73*>          The upper or lower triangle of the Hermitian matrix A, packed
74*>          columnwise in a linear array.  The j-th column of A is stored
75*>          in the array AP as follows:
76*>          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
77*>          if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.
78*> \endverbatim
79*>
80*> \param[in] B
81*> \verbatim
82*>          B is COMPLEX*16 array, dimension (LDB,NRHS)
83*>          The right hand side vectors for the system of linear
84*>          equations.
85*> \endverbatim
86*>
87*> \param[in] LDB
88*> \verbatim
89*>          LDB is INTEGER
90*>          The leading dimension of the array B.  LDB >= max(1,N).
91*> \endverbatim
92*>
93*> \param[in] X
94*> \verbatim
95*>          X is COMPLEX*16 array, dimension (LDX,NRHS)
96*>          The computed solution vectors.  Each vector is stored as a
97*>          column of the matrix X.
98*> \endverbatim
99*>
100*> \param[in] LDX
101*> \verbatim
102*>          LDX is INTEGER
103*>          The leading dimension of the array X.  LDX >= max(1,N).
104*> \endverbatim
105*>
106*> \param[in] XACT
107*> \verbatim
108*>          XACT is COMPLEX*16 array, dimension (LDX,NRHS)
109*>          The exact solution vectors.  Each vector is stored as a
110*>          column of the matrix XACT.
111*> \endverbatim
112*>
113*> \param[in] LDXACT
114*> \verbatim
115*>          LDXACT is INTEGER
116*>          The leading dimension of the array XACT.  LDXACT >= max(1,N).
117*> \endverbatim
118*>
119*> \param[in] FERR
120*> \verbatim
121*>          FERR is DOUBLE PRECISION array, dimension (NRHS)
122*>          The estimated forward error bounds for each solution vector
123*>          X.  If XTRUE is the true solution, FERR bounds the magnitude
124*>          of the largest entry in (X - XTRUE) divided by the magnitude
125*>          of the largest entry in X.
126*> \endverbatim
127*>
128*> \param[in] BERR
129*> \verbatim
130*>          BERR is DOUBLE PRECISION array, dimension (NRHS)
131*>          The componentwise relative backward error of each solution
132*>          vector (i.e., the smallest relative change in any entry of A
133*>          or B that makes X an exact solution).
134*> \endverbatim
135*>
136*> \param[out] RESLTS
137*> \verbatim
138*>          RESLTS is DOUBLE PRECISION array, dimension (2)
139*>          The maximum over the NRHS solution vectors of the ratios:
140*>          RESLTS(1) = norm(X - XACT) / ( norm(X) * FERR )
141*>          RESLTS(2) = BERR / ( (n+1)*EPS + (*) )
142*> \endverbatim
143*
144*  Authors:
145*  ========
146*
147*> \author Univ. of Tennessee
148*> \author Univ. of California Berkeley
149*> \author Univ. of Colorado Denver
150*> \author NAG Ltd.
151*
152*> \date December 2016
153*
154*> \ingroup complex16_lin
155*
156*  =====================================================================
157      SUBROUTINE ZPPT05( UPLO, N, NRHS, AP, B, LDB, X, LDX, XACT,
158     $                   LDXACT, FERR, BERR, RESLTS )
159*
160*  -- LAPACK test routine (version 3.7.0) --
161*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
162*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
163*     December 2016
164*
165*     .. Scalar Arguments ..
166      CHARACTER          UPLO
167      INTEGER            LDB, LDX, LDXACT, N, NRHS
168*     ..
169*     .. Array Arguments ..
170      DOUBLE PRECISION   BERR( * ), FERR( * ), RESLTS( * )
171      COMPLEX*16         AP( * ), B( LDB, * ), X( LDX, * ),
172     $                   XACT( LDXACT, * )
173*     ..
174*
175*  =====================================================================
176*
177*     .. Parameters ..
178      DOUBLE PRECISION   ZERO, ONE
179      PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
180*     ..
181*     .. Local Scalars ..
182      LOGICAL            UPPER
183      INTEGER            I, IMAX, J, JC, K
184      DOUBLE PRECISION   AXBI, DIFF, EPS, ERRBND, OVFL, TMP, UNFL, XNORM
185      COMPLEX*16         ZDUM
186*     ..
187*     .. External Functions ..
188      LOGICAL            LSAME
189      INTEGER            IZAMAX
190      DOUBLE PRECISION   DLAMCH
191      EXTERNAL           LSAME, IZAMAX, DLAMCH
192*     ..
193*     .. Intrinsic Functions ..
194      INTRINSIC          ABS, DBLE, DIMAG, MAX, MIN
195*     ..
196*     .. Statement Functions ..
197      DOUBLE PRECISION   CABS1
198*     ..
199*     .. Statement Function definitions ..
200      CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
201*     ..
202*     .. Executable Statements ..
203*
204*     Quick exit if N = 0 or NRHS = 0.
205*
206      IF( N.LE.0 .OR. NRHS.LE.0 ) THEN
207         RESLTS( 1 ) = ZERO
208         RESLTS( 2 ) = ZERO
209         RETURN
210      END IF
211*
212      EPS = DLAMCH( 'Epsilon' )
213      UNFL = DLAMCH( 'Safe minimum' )
214      OVFL = ONE / UNFL
215      UPPER = LSAME( UPLO, 'U' )
216*
217*     Test 1:  Compute the maximum of
218*        norm(X - XACT) / ( norm(X) * FERR )
219*     over all the vectors X and XACT using the infinity-norm.
220*
221      ERRBND = ZERO
222      DO 30 J = 1, NRHS
223         IMAX = IZAMAX( N, X( 1, J ), 1 )
224         XNORM = MAX( CABS1( X( IMAX, J ) ), UNFL )
225         DIFF = ZERO
226         DO 10 I = 1, N
227            DIFF = MAX( DIFF, CABS1( X( I, J )-XACT( I, J ) ) )
228   10    CONTINUE
229*
230         IF( XNORM.GT.ONE ) THEN
231            GO TO 20
232         ELSE IF( DIFF.LE.OVFL*XNORM ) THEN
233            GO TO 20
234         ELSE
235            ERRBND = ONE / EPS
236            GO TO 30
237         END IF
238*
239   20    CONTINUE
240         IF( DIFF / XNORM.LE.FERR( J ) ) THEN
241            ERRBND = MAX( ERRBND, ( DIFF / XNORM ) / FERR( J ) )
242         ELSE
243            ERRBND = ONE / EPS
244         END IF
245   30 CONTINUE
246      RESLTS( 1 ) = ERRBND
247*
248*     Test 2:  Compute the maximum of BERR / ( (n+1)*EPS + (*) ), where
249*     (*) = (n+1)*UNFL / (min_i (abs(A)*abs(X) +abs(b))_i )
250*
251      DO 90 K = 1, NRHS
252         DO 80 I = 1, N
253            TMP = CABS1( B( I, K ) )
254            IF( UPPER ) THEN
255               JC = ( ( I-1 )*I ) / 2
256               DO 40 J = 1, I - 1
257                  TMP = TMP + CABS1( AP( JC+J ) )*CABS1( X( J, K ) )
258   40          CONTINUE
259               TMP = TMP + ABS( DBLE( AP( JC+I ) ) )*CABS1( X( I, K ) )
260               JC = JC + I + I
261               DO 50 J = I + 1, N
262                  TMP = TMP + CABS1( AP( JC ) )*CABS1( X( J, K ) )
263                  JC = JC + J
264   50          CONTINUE
265            ELSE
266               JC = I
267               DO 60 J = 1, I - 1
268                  TMP = TMP + CABS1( AP( JC ) )*CABS1( X( J, K ) )
269                  JC = JC + N - J
270   60          CONTINUE
271               TMP = TMP + ABS( DBLE( AP( JC ) ) )*CABS1( X( I, K ) )
272               DO 70 J = I + 1, N
273                  TMP = TMP + CABS1( AP( JC+J-I ) )*CABS1( X( J, K ) )
274   70          CONTINUE
275            END IF
276            IF( I.EQ.1 ) THEN
277               AXBI = TMP
278            ELSE
279               AXBI = MIN( AXBI, TMP )
280            END IF
281   80    CONTINUE
282         TMP = BERR( K ) / ( ( N+1 )*EPS+( N+1 )*UNFL /
283     $         MAX( AXBI, ( N+1 )*UNFL ) )
284         IF( K.EQ.1 ) THEN
285            RESLTS( 2 ) = TMP
286         ELSE
287            RESLTS( 2 ) = MAX( RESLTS( 2 ), TMP )
288         END IF
289   90 CONTINUE
290*
291      RETURN
292*
293*     End of ZPPT05
294*
295      END
296