1*> \brief \b SGET39
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 SGET39( RMAX, LMAX, NINFO, KNT )
12*
13*       .. Scalar Arguments ..
14*       INTEGER            KNT, LMAX, NINFO
15*       REAL               RMAX
16*       ..
17*
18*
19*> \par Purpose:
20*  =============
21*>
22*> \verbatim
23*>
24*> SGET39 tests SLAQTR, a routine for solving the real or
25*> special complex quasi upper triangular system
26*>
27*>      op(T)*p = scale*c,
28*> or
29*>      op(T + iB)*(p+iq) = scale*(c+id),
30*>
31*> in real arithmetic. T is upper quasi-triangular.
32*> If it is complex, then the first diagonal block of T must be
33*> 1 by 1, B has the special structure
34*>
35*>                B = [ b(1) b(2) ... b(n) ]
36*>                    [       w            ]
37*>                    [           w        ]
38*>                    [              .     ]
39*>                    [                 w  ]
40*>
41*> op(A) = A or A', where A' denotes the conjugate transpose of
42*> the matrix A.
43*>
44*> On input, X = [ c ].  On output, X = [ p ].
45*>               [ d ]                  [ q ]
46*>
47*> Scale is an output less than or equal to 1, chosen to avoid
48*> overflow in X.
49*> This subroutine is specially designed for the condition number
50*> estimation in the eigenproblem routine STRSNA.
51*>
52*> The test code verifies that the following residual is order 1:
53*>
54*>      ||(T+i*B)*(x1+i*x2) - scale*(d1+i*d2)||
55*>    -----------------------------------------
56*>        max(ulp*(||T||+||B||)*(||x1||+||x2||),
57*>            (||T||+||B||)*smlnum/ulp,
58*>            smlnum)
59*>
60*> (The (||T||+||B||)*smlnum/ulp term accounts for possible
61*>  (gradual or nongradual) underflow in x1 and x2.)
62*> \endverbatim
63*
64*  Arguments:
65*  ==========
66*
67*> \param[out] RMAX
68*> \verbatim
69*>          RMAX is REAL
70*>          Value of the largest test ratio.
71*> \endverbatim
72*>
73*> \param[out] LMAX
74*> \verbatim
75*>          LMAX is INTEGER
76*>          Example number where largest test ratio achieved.
77*> \endverbatim
78*>
79*> \param[out] NINFO
80*> \verbatim
81*>          NINFO is INTEGER
82*>          Number of examples where INFO is nonzero.
83*> \endverbatim
84*>
85*> \param[out] KNT
86*> \verbatim
87*>          KNT is INTEGER
88*>          Total number of examples tested.
89*> \endverbatim
90*
91*  Authors:
92*  ========
93*
94*> \author Univ. of Tennessee
95*> \author Univ. of California Berkeley
96*> \author Univ. of Colorado Denver
97*> \author NAG Ltd.
98*
99*> \ingroup single_eig
100*
101*  =====================================================================
102      SUBROUTINE SGET39( RMAX, LMAX, NINFO, KNT )
103*
104*  -- LAPACK test routine --
105*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
106*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
107*
108*     .. Scalar Arguments ..
109      INTEGER            KNT, LMAX, NINFO
110      REAL               RMAX
111*     ..
112*
113*  =====================================================================
114*
115*     .. Parameters ..
116      INTEGER            LDT, LDT2
117      PARAMETER          ( LDT = 10, LDT2 = 2*LDT )
118      REAL               ZERO, ONE
119      PARAMETER          ( ZERO = 0.0, ONE = 1.0 )
120*     ..
121*     .. Local Scalars ..
122      INTEGER            I, INFO, IVM1, IVM2, IVM3, IVM4, IVM5, J, K, N,
123     $                   NDIM
124      REAL               BIGNUM, DOMIN, DUMM, EPS, NORM, NORMTB, RESID,
125     $                   SCALE, SMLNUM, W, XNORM
126*     ..
127*     .. External Functions ..
128      INTEGER            ISAMAX
129      REAL               SASUM, SDOT, SLAMCH, SLANGE
130      EXTERNAL           ISAMAX, SASUM, SDOT, SLAMCH, SLANGE
131*     ..
132*     .. External Subroutines ..
133      EXTERNAL           SCOPY, SGEMV, SLABAD, SLAQTR
134*     ..
135*     .. Intrinsic Functions ..
136      INTRINSIC          ABS, COS, MAX, REAL, SIN, SQRT
137*     ..
138*     .. Local Arrays ..
139      INTEGER            IDIM( 6 ), IVAL( 5, 5, 6 )
140      REAL               B( LDT ), D( LDT2 ), DUM( 1 ), T( LDT, LDT ),
141     $                   VM1( 5 ), VM2( 5 ), VM3( 5 ), VM4( 5 ),
142     $                   VM5( 3 ), WORK( LDT ), X( LDT2 ), Y( LDT2 )
143*     ..
144*     .. Data statements ..
145      DATA               IDIM / 4, 5*5 /
146      DATA               IVAL / 3, 4*0, 1, 1, -1, 0, 0, 3, 2, 1, 0, 0,
147     $                   4, 3, 2, 2, 0, 5*0, 1, 4*0, 2, 2, 3*0, 3, 3, 4,
148     $                   0, 0, 4, 2, 2, 3, 0, 4*1, 5, 1, 4*0, 2, 4, -2,
149     $                   0, 0, 3, 3, 4, 0, 0, 4, 2, 2, 3, 0, 5*1, 1,
150     $                   4*0, 2, 1, -1, 0, 0, 9, 8, 1, 0, 0, 4, 9, 1, 2,
151     $                   -1, 5*2, 9, 4*0, 6, 4, 0, 0, 0, 3, 2, 1, 1, 0,
152     $                   5, 1, -1, 1, 0, 5*2, 4, 4*0, 2, 2, 0, 0, 0, 1,
153     $                   4, 4, 0, 0, 2, 4, 2, 2, -1, 5*2 /
154*     ..
155*     .. Executable Statements ..
156*
157*     Get machine parameters
158*
159      EPS = SLAMCH( 'P' )
160      SMLNUM = SLAMCH( 'S' )
161      BIGNUM = ONE / SMLNUM
162      CALL SLABAD( SMLNUM, BIGNUM )
163*
164*     Set up test case parameters
165*
166      VM1( 1 ) = ONE
167      VM1( 2 ) = SQRT( SMLNUM )
168      VM1( 3 ) = SQRT( VM1( 2 ) )
169      VM1( 4 ) = SQRT( BIGNUM )
170      VM1( 5 ) = SQRT( VM1( 4 ) )
171*
172      VM2( 1 ) = ONE
173      VM2( 2 ) = SQRT( SMLNUM )
174      VM2( 3 ) = SQRT( VM2( 2 ) )
175      VM2( 4 ) = SQRT( BIGNUM )
176      VM2( 5 ) = SQRT( VM2( 4 ) )
177*
178      VM3( 1 ) = ONE
179      VM3( 2 ) = SQRT( SMLNUM )
180      VM3( 3 ) = SQRT( VM3( 2 ) )
181      VM3( 4 ) = SQRT( BIGNUM )
182      VM3( 5 ) = SQRT( VM3( 4 ) )
183*
184      VM4( 1 ) = ONE
185      VM4( 2 ) = SQRT( SMLNUM )
186      VM4( 3 ) = SQRT( VM4( 2 ) )
187      VM4( 4 ) = SQRT( BIGNUM )
188      VM4( 5 ) = SQRT( VM4( 4 ) )
189*
190      VM5( 1 ) = ONE
191      VM5( 2 ) = EPS
192      VM5( 3 ) = SQRT( SMLNUM )
193*
194*     Initialization
195*
196      KNT = 0
197      RMAX = ZERO
198      NINFO = 0
199      SMLNUM = SMLNUM / EPS
200*
201*     Begin test loop
202*
203      DO 140 IVM5 = 1, 3
204         DO 130 IVM4 = 1, 5
205            DO 120 IVM3 = 1, 5
206               DO 110 IVM2 = 1, 5
207                  DO 100 IVM1 = 1, 5
208                     DO 90 NDIM = 1, 6
209*
210                        N = IDIM( NDIM )
211                        DO 20 I = 1, N
212                           DO 10 J = 1, N
213                              T( I, J ) = REAL( IVAL( I, J, NDIM ) )*
214     $                                    VM1( IVM1 )
215                              IF( I.GE.J )
216     $                           T( I, J ) = T( I, J )*VM5( IVM5 )
217   10                      CONTINUE
218   20                   CONTINUE
219*
220                        W = ONE*VM2( IVM2 )
221*
222                        DO 30 I = 1, N
223                           B( I ) = COS( REAL( I ) )*VM3( IVM3 )
224   30                   CONTINUE
225*
226                        DO 40 I = 1, 2*N
227                           D( I ) = SIN( REAL( I ) )*VM4( IVM4 )
228   40                   CONTINUE
229*
230                        NORM = SLANGE( '1', N, N, T, LDT, WORK )
231                        K = ISAMAX( N, B, 1 )
232                        NORMTB = NORM + ABS( B( K ) ) + ABS( W )
233*
234                        CALL SCOPY( N, D, 1, X, 1 )
235                        KNT = KNT + 1
236                        CALL SLAQTR( .FALSE., .TRUE., N, T, LDT, DUM,
237     $                               DUMM, SCALE, X, WORK, INFO )
238                        IF( INFO.NE.0 )
239     $                     NINFO = NINFO + 1
240*
241*                       || T*x - scale*d || /
242*                         max(ulp*||T||*||x||,smlnum/ulp*||T||,smlnum)
243*
244                        CALL SCOPY( N, D, 1, Y, 1 )
245                        CALL SGEMV( 'No transpose', N, N, ONE, T, LDT,
246     $                              X, 1, -SCALE, Y, 1 )
247                        XNORM = SASUM( N, X, 1 )
248                        RESID = SASUM( N, Y, 1 )
249                        DOMIN = MAX( SMLNUM, ( SMLNUM / EPS )*NORM,
250     $                          ( NORM*EPS )*XNORM )
251                        RESID = RESID / DOMIN
252                        IF( RESID.GT.RMAX ) THEN
253                           RMAX = RESID
254                           LMAX = KNT
255                        END IF
256*
257                        CALL SCOPY( N, D, 1, X, 1 )
258                        KNT = KNT + 1
259                        CALL SLAQTR( .TRUE., .TRUE., N, T, LDT, DUM,
260     $                               DUMM, SCALE, X, WORK, INFO )
261                        IF( INFO.NE.0 )
262     $                     NINFO = NINFO + 1
263*
264*                       || T*x - scale*d || /
265*                         max(ulp*||T||*||x||,smlnum/ulp*||T||,smlnum)
266*
267                        CALL SCOPY( N, D, 1, Y, 1 )
268                        CALL SGEMV( 'Transpose', N, N, ONE, T, LDT, X,
269     $                              1, -SCALE, Y, 1 )
270                        XNORM = SASUM( N, X, 1 )
271                        RESID = SASUM( N, Y, 1 )
272                        DOMIN = MAX( SMLNUM, ( SMLNUM / EPS )*NORM,
273     $                          ( NORM*EPS )*XNORM )
274                        RESID = RESID / DOMIN
275                        IF( RESID.GT.RMAX ) THEN
276                           RMAX = RESID
277                           LMAX = KNT
278                        END IF
279*
280                        CALL SCOPY( 2*N, D, 1, X, 1 )
281                        KNT = KNT + 1
282                        CALL SLAQTR( .FALSE., .FALSE., N, T, LDT, B, W,
283     $                               SCALE, X, WORK, INFO )
284                        IF( INFO.NE.0 )
285     $                     NINFO = NINFO + 1
286*
287*                       ||(T+i*B)*(x1+i*x2) - scale*(d1+i*d2)|| /
288*                          max(ulp*(||T||+||B||)*(||x1||+||x2||),
289*                                  smlnum/ulp * (||T||+||B||), smlnum )
290*
291*
292                        CALL SCOPY( 2*N, D, 1, Y, 1 )
293                        Y( 1 ) = SDOT( N, B, 1, X( 1+N ), 1 ) +
294     $                           SCALE*Y( 1 )
295                        DO 50 I = 2, N
296                           Y( I ) = W*X( I+N ) + SCALE*Y( I )
297   50                   CONTINUE
298                        CALL SGEMV( 'No transpose', N, N, ONE, T, LDT,
299     $                              X, 1, -ONE, Y, 1 )
300*
301                        Y( 1+N ) = SDOT( N, B, 1, X, 1 ) -
302     $                             SCALE*Y( 1+N )
303                        DO 60 I = 2, N
304                           Y( I+N ) = W*X( I ) - SCALE*Y( I+N )
305   60                   CONTINUE
306                        CALL SGEMV( 'No transpose', N, N, ONE, T, LDT,
307     $                              X( 1+N ), 1, ONE, Y( 1+N ), 1 )
308*
309                        RESID = SASUM( 2*N, Y, 1 )
310                        DOMIN = MAX( SMLNUM, ( SMLNUM / EPS )*NORMTB,
311     $                          EPS*( NORMTB*SASUM( 2*N, X, 1 ) ) )
312                        RESID = RESID / DOMIN
313                        IF( RESID.GT.RMAX ) THEN
314                           RMAX = RESID
315                           LMAX = KNT
316                        END IF
317*
318                        CALL SCOPY( 2*N, D, 1, X, 1 )
319                        KNT = KNT + 1
320                        CALL SLAQTR( .TRUE., .FALSE., N, T, LDT, B, W,
321     $                               SCALE, X, WORK, INFO )
322                        IF( INFO.NE.0 )
323     $                     NINFO = NINFO + 1
324*
325*                       ||(T+i*B)*(x1+i*x2) - scale*(d1+i*d2)|| /
326*                          max(ulp*(||T||+||B||)*(||x1||+||x2||),
327*                                  smlnum/ulp * (||T||+||B||), smlnum )
328*
329                        CALL SCOPY( 2*N, D, 1, Y, 1 )
330                        Y( 1 ) = B( 1 )*X( 1+N ) - SCALE*Y( 1 )
331                        DO 70 I = 2, N
332                           Y( I ) = B( I )*X( 1+N ) + W*X( I+N ) -
333     $                              SCALE*Y( I )
334   70                   CONTINUE
335                        CALL SGEMV( 'Transpose', N, N, ONE, T, LDT, X,
336     $                              1, ONE, Y, 1 )
337*
338                        Y( 1+N ) = B( 1 )*X( 1 ) + SCALE*Y( 1+N )
339                        DO 80 I = 2, N
340                           Y( I+N ) = B( I )*X( 1 ) + W*X( I ) +
341     $                                SCALE*Y( I+N )
342   80                   CONTINUE
343                        CALL SGEMV( 'Transpose', N, N, ONE, T, LDT,
344     $                              X( 1+N ), 1, -ONE, Y( 1+N ), 1 )
345*
346                        RESID = SASUM( 2*N, Y, 1 )
347                        DOMIN = MAX( SMLNUM, ( SMLNUM / EPS )*NORMTB,
348     $                          EPS*( NORMTB*SASUM( 2*N, X, 1 ) ) )
349                        RESID = RESID / DOMIN
350                        IF( RESID.GT.RMAX ) THEN
351                           RMAX = RESID
352                           LMAX = KNT
353                        END IF
354*
355   90                CONTINUE
356  100             CONTINUE
357  110          CONTINUE
358  120       CONTINUE
359  130    CONTINUE
360  140 CONTINUE
361*
362      RETURN
363*
364*     End of SGET39
365*
366      END
367