LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
ddrvgt.f
Go to the documentation of this file.
1 *> \brief \b DDRVGT
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 DDRVGT( DOTYPE, NN, NVAL, NRHS, THRESH, TSTERR, A, AF,
12 * B, X, XACT, WORK, RWORK, IWORK, NOUT )
13 *
14 * .. Scalar Arguments ..
15 * LOGICAL TSTERR
16 * INTEGER NN, NOUT, NRHS
17 * DOUBLE PRECISION THRESH
18 * ..
19 * .. Array Arguments ..
20 * LOGICAL DOTYPE( * )
21 * INTEGER IWORK( * ), NVAL( * )
22 * DOUBLE PRECISION A( * ), AF( * ), B( * ), RWORK( * ), WORK( * ),
23 * $ X( * ), XACT( * )
24 * ..
25 *
26 *
27 *> \par Purpose:
28 * =============
29 *>
30 *> \verbatim
31 *>
32 *> DDRVGT tests DGTSV and -SVX.
33 *> \endverbatim
34 *
35 * Arguments:
36 * ==========
37 *
38 *> \param[in] DOTYPE
39 *> \verbatim
40 *> DOTYPE is LOGICAL array, dimension (NTYPES)
41 *> The matrix types to be used for testing. Matrices of type j
42 *> (for 1 <= j <= NTYPES) are used for testing if DOTYPE(j) =
43 *> .TRUE.; if DOTYPE(j) = .FALSE., then type j is not used.
44 *> \endverbatim
45 *>
46 *> \param[in] NN
47 *> \verbatim
48 *> NN is INTEGER
49 *> The number of values of N contained in the vector NVAL.
50 *> \endverbatim
51 *>
52 *> \param[in] NVAL
53 *> \verbatim
54 *> NVAL is INTEGER array, dimension (NN)
55 *> The values of the matrix dimension N.
56 *> \endverbatim
57 *>
58 *> \param[in] NRHS
59 *> \verbatim
60 *> NRHS is INTEGER
61 *> The number of right hand sides, NRHS >= 0.
62 *> \endverbatim
63 *>
64 *> \param[in] THRESH
65 *> \verbatim
66 *> THRESH is DOUBLE PRECISION
67 *> The threshold value for the test ratios. A result is
68 *> included in the output file if RESULT >= THRESH. To have
69 *> every test ratio printed, use THRESH = 0.
70 *> \endverbatim
71 *>
72 *> \param[in] TSTERR
73 *> \verbatim
74 *> TSTERR is LOGICAL
75 *> Flag that indicates whether error exits are to be tested.
76 *> \endverbatim
77 *>
78 *> \param[out] A
79 *> \verbatim
80 *> A is DOUBLE PRECISION array, dimension (NMAX*4)
81 *> \endverbatim
82 *>
83 *> \param[out] AF
84 *> \verbatim
85 *> AF is DOUBLE PRECISION array, dimension (NMAX*4)
86 *> \endverbatim
87 *>
88 *> \param[out] B
89 *> \verbatim
90 *> B is DOUBLE PRECISION array, dimension (NMAX*NRHS)
91 *> \endverbatim
92 *>
93 *> \param[out] X
94 *> \verbatim
95 *> X is DOUBLE PRECISION array, dimension (NMAX*NRHS)
96 *> \endverbatim
97 *>
98 *> \param[out] XACT
99 *> \verbatim
100 *> XACT is DOUBLE PRECISION array, dimension (NMAX*NRHS)
101 *> \endverbatim
102 *>
103 *> \param[out] WORK
104 *> \verbatim
105 *> WORK is DOUBLE PRECISION array, dimension
106 *> (NMAX*max(3,NRHS))
107 *> \endverbatim
108 *>
109 *> \param[out] RWORK
110 *> \verbatim
111 *> RWORK is DOUBLE PRECISION array, dimension
112 *> (max(NMAX,2*NRHS))
113 *> \endverbatim
114 *>
115 *> \param[out] IWORK
116 *> \verbatim
117 *> IWORK is INTEGER array, dimension (2*NMAX)
118 *> \endverbatim
119 *>
120 *> \param[in] NOUT
121 *> \verbatim
122 *> NOUT is INTEGER
123 *> The unit number for output.
124 *> \endverbatim
125 *
126 * Authors:
127 * ========
128 *
129 *> \author Univ. of Tennessee
130 *> \author Univ. of California Berkeley
131 *> \author Univ. of Colorado Denver
132 *> \author NAG Ltd.
133 *
134 *> \date November 2011
135 *
136 *> \ingroup double_lin
137 *
138 * =====================================================================
139  SUBROUTINE ddrvgt( DOTYPE, NN, NVAL, NRHS, THRESH, TSTERR, A, AF,
140  $ b, x, xact, work, rwork, iwork, nout )
141 *
142 * -- LAPACK test routine (version 3.4.0) --
143 * -- LAPACK is a software package provided by Univ. of Tennessee, --
144 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
145 * November 2011
146 *
147 * .. Scalar Arguments ..
148  LOGICAL tsterr
149  INTEGER nn, nout, nrhs
150  DOUBLE PRECISION thresh
151 * ..
152 * .. Array Arguments ..
153  LOGICAL dotype( * )
154  INTEGER iwork( * ), nval( * )
155  DOUBLE PRECISION a( * ), af( * ), b( * ), rwork( * ), work( * ),
156  $ x( * ), xact( * )
157 * ..
158 *
159 * =====================================================================
160 *
161 * .. Parameters ..
162  DOUBLE PRECISION one, zero
163  parameter( one = 1.0d+0, zero = 0.0d+0 )
164  INTEGER ntypes
165  parameter( ntypes = 12 )
166  INTEGER ntests
167  parameter( ntests = 6 )
168 * ..
169 * .. Local Scalars ..
170  LOGICAL trfcon, zerot
171  CHARACTER dist, fact, trans, type
172  CHARACTER*3 path
173  INTEGER i, ifact, imat, in, info, itran, ix, izero, j,
174  $ k, k1, kl, koff, ku, lda, m, mode, n, nerrs,
175  $ nfail, nimat, nrun, nt
176  DOUBLE PRECISION ainvnm, anorm, anormi, anormo, cond, rcond,
177  $ rcondc, rcondi, rcondo
178 * ..
179 * .. Local Arrays ..
180  CHARACTER transs( 3 )
181  INTEGER iseed( 4 ), iseedy( 4 )
182  DOUBLE PRECISION result( ntests ), z( 3 )
183 * ..
184 * .. External Functions ..
185  DOUBLE PRECISION dasum, dget06, dlangt
186  EXTERNAL dasum, dget06, dlangt
187 * ..
188 * .. External Subroutines ..
189  EXTERNAL aladhd, alaerh, alasvm, dcopy, derrvx, dget04,
190  $ dgtsv, dgtsvx, dgtt01, dgtt02, dgtt05, dgttrf,
191  $ dgttrs, dlacpy, dlagtm, dlarnv, dlaset, dlatb4,
192  $ dlatms, dscal
193 * ..
194 * .. Intrinsic Functions ..
195  INTRINSIC max
196 * ..
197 * .. Scalars in Common ..
198  LOGICAL lerr, ok
199  CHARACTER*32 srnamt
200  INTEGER infot, nunit
201 * ..
202 * .. Common blocks ..
203  common / infoc / infot, nunit, ok, lerr
204  common / srnamc / srnamt
205 * ..
206 * .. Data statements ..
207  DATA iseedy / 0, 0, 0, 1 / , transs / 'N', 'T',
208  $ 'C' /
209 * ..
210 * .. Executable Statements ..
211 *
212  path( 1: 1 ) = 'Double precision'
213  path( 2: 3 ) = 'GT'
214  nrun = 0
215  nfail = 0
216  nerrs = 0
217  DO 10 i = 1, 4
218  iseed( i ) = iseedy( i )
219  10 continue
220 *
221 * Test the error exits
222 *
223  IF( tsterr )
224  $ CALL derrvx( path, nout )
225  infot = 0
226 *
227  DO 140 in = 1, nn
228 *
229 * Do for each value of N in NVAL.
230 *
231  n = nval( in )
232  m = max( n-1, 0 )
233  lda = max( 1, n )
234  nimat = ntypes
235  IF( n.LE.0 )
236  $ nimat = 1
237 *
238  DO 130 imat = 1, nimat
239 *
240 * Do the tests only if DOTYPE( IMAT ) is true.
241 *
242  IF( .NOT.dotype( imat ) )
243  $ go to 130
244 *
245 * Set up parameters with DLATB4.
246 *
247  CALL dlatb4( path, imat, n, n, type, kl, ku, anorm, mode,
248  $ cond, dist )
249 *
250  zerot = imat.GE.8 .AND. imat.LE.10
251  IF( imat.LE.6 ) THEN
252 *
253 * Types 1-6: generate matrices of known condition number.
254 *
255  koff = max( 2-ku, 3-max( 1, n ) )
256  srnamt = 'DLATMS'
257  CALL dlatms( n, n, dist, iseed, type, rwork, mode, cond,
258  $ anorm, kl, ku, 'Z', af( koff ), 3, work,
259  $ info )
260 *
261 * Check the error code from DLATMS.
262 *
263  IF( info.NE.0 ) THEN
264  CALL alaerh( path, 'DLATMS', info, 0, ' ', n, n, kl,
265  $ ku, -1, imat, nfail, nerrs, nout )
266  go to 130
267  END IF
268  izero = 0
269 *
270  IF( n.GT.1 ) THEN
271  CALL dcopy( n-1, af( 4 ), 3, a, 1 )
272  CALL dcopy( n-1, af( 3 ), 3, a( n+m+1 ), 1 )
273  END IF
274  CALL dcopy( n, af( 2 ), 3, a( m+1 ), 1 )
275  ELSE
276 *
277 * Types 7-12: generate tridiagonal matrices with
278 * unknown condition numbers.
279 *
280  IF( .NOT.zerot .OR. .NOT.dotype( 7 ) ) THEN
281 *
282 * Generate a matrix with elements from [-1,1].
283 *
284  CALL dlarnv( 2, iseed, n+2*m, a )
285  IF( anorm.NE.one )
286  $ CALL dscal( n+2*m, anorm, a, 1 )
287  ELSE IF( izero.GT.0 ) THEN
288 *
289 * Reuse the last matrix by copying back the zeroed out
290 * elements.
291 *
292  IF( izero.EQ.1 ) THEN
293  a( n ) = z( 2 )
294  IF( n.GT.1 )
295  $ a( 1 ) = z( 3 )
296  ELSE IF( izero.EQ.n ) THEN
297  a( 3*n-2 ) = z( 1 )
298  a( 2*n-1 ) = z( 2 )
299  ELSE
300  a( 2*n-2+izero ) = z( 1 )
301  a( n-1+izero ) = z( 2 )
302  a( izero ) = z( 3 )
303  END IF
304  END IF
305 *
306 * If IMAT > 7, set one column of the matrix to 0.
307 *
308  IF( .NOT.zerot ) THEN
309  izero = 0
310  ELSE IF( imat.EQ.8 ) THEN
311  izero = 1
312  z( 2 ) = a( n )
313  a( n ) = zero
314  IF( n.GT.1 ) THEN
315  z( 3 ) = a( 1 )
316  a( 1 ) = zero
317  END IF
318  ELSE IF( imat.EQ.9 ) THEN
319  izero = n
320  z( 1 ) = a( 3*n-2 )
321  z( 2 ) = a( 2*n-1 )
322  a( 3*n-2 ) = zero
323  a( 2*n-1 ) = zero
324  ELSE
325  izero = ( n+1 ) / 2
326  DO 20 i = izero, n - 1
327  a( 2*n-2+i ) = zero
328  a( n-1+i ) = zero
329  a( i ) = zero
330  20 continue
331  a( 3*n-2 ) = zero
332  a( 2*n-1 ) = zero
333  END IF
334  END IF
335 *
336  DO 120 ifact = 1, 2
337  IF( ifact.EQ.1 ) THEN
338  fact = 'F'
339  ELSE
340  fact = 'N'
341  END IF
342 *
343 * Compute the condition number for comparison with
344 * the value returned by DGTSVX.
345 *
346  IF( zerot ) THEN
347  IF( ifact.EQ.1 )
348  $ go to 120
349  rcondo = zero
350  rcondi = zero
351 *
352  ELSE IF( ifact.EQ.1 ) THEN
353  CALL dcopy( n+2*m, a, 1, af, 1 )
354 *
355 * Compute the 1-norm and infinity-norm of A.
356 *
357  anormo = dlangt( '1', n, a, a( m+1 ), a( n+m+1 ) )
358  anormi = dlangt( 'I', n, a, a( m+1 ), a( n+m+1 ) )
359 *
360 * Factor the matrix A.
361 *
362  CALL dgttrf( n, af, af( m+1 ), af( n+m+1 ),
363  $ af( n+2*m+1 ), iwork, info )
364 *
365 * Use DGTTRS to solve for one column at a time of
366 * inv(A), computing the maximum column sum as we go.
367 *
368  ainvnm = zero
369  DO 40 i = 1, n
370  DO 30 j = 1, n
371  x( j ) = zero
372  30 continue
373  x( i ) = one
374  CALL dgttrs( 'No transpose', n, 1, af, af( m+1 ),
375  $ af( n+m+1 ), af( n+2*m+1 ), iwork, x,
376  $ lda, info )
377  ainvnm = max( ainvnm, dasum( n, x, 1 ) )
378  40 continue
379 *
380 * Compute the 1-norm condition number of A.
381 *
382  IF( anormo.LE.zero .OR. ainvnm.LE.zero ) THEN
383  rcondo = one
384  ELSE
385  rcondo = ( one / anormo ) / ainvnm
386  END IF
387 *
388 * Use DGTTRS to solve for one column at a time of
389 * inv(A'), computing the maximum column sum as we go.
390 *
391  ainvnm = zero
392  DO 60 i = 1, n
393  DO 50 j = 1, n
394  x( j ) = zero
395  50 continue
396  x( i ) = one
397  CALL dgttrs( 'Transpose', n, 1, af, af( m+1 ),
398  $ af( n+m+1 ), af( n+2*m+1 ), iwork, x,
399  $ lda, info )
400  ainvnm = max( ainvnm, dasum( n, x, 1 ) )
401  60 continue
402 *
403 * Compute the infinity-norm condition number of A.
404 *
405  IF( anormi.LE.zero .OR. ainvnm.LE.zero ) THEN
406  rcondi = one
407  ELSE
408  rcondi = ( one / anormi ) / ainvnm
409  END IF
410  END IF
411 *
412  DO 110 itran = 1, 3
413  trans = transs( itran )
414  IF( itran.EQ.1 ) THEN
415  rcondc = rcondo
416  ELSE
417  rcondc = rcondi
418  END IF
419 *
420 * Generate NRHS random solution vectors.
421 *
422  ix = 1
423  DO 70 j = 1, nrhs
424  CALL dlarnv( 2, iseed, n, xact( ix ) )
425  ix = ix + lda
426  70 continue
427 *
428 * Set the right hand side.
429 *
430  CALL dlagtm( trans, n, nrhs, one, a, a( m+1 ),
431  $ a( n+m+1 ), xact, lda, zero, b, lda )
432 *
433  IF( ifact.EQ.2 .AND. itran.EQ.1 ) THEN
434 *
435 * --- Test DGTSV ---
436 *
437 * Solve the system using Gaussian elimination with
438 * partial pivoting.
439 *
440  CALL dcopy( n+2*m, a, 1, af, 1 )
441  CALL dlacpy( 'Full', n, nrhs, b, lda, x, lda )
442 *
443  srnamt = 'DGTSV '
444  CALL dgtsv( n, nrhs, af, af( m+1 ), af( n+m+1 ), x,
445  $ lda, info )
446 *
447 * Check error code from DGTSV .
448 *
449  IF( info.NE.izero )
450  $ CALL alaerh( path, 'DGTSV ', info, izero, ' ',
451  $ n, n, 1, 1, nrhs, imat, nfail,
452  $ nerrs, nout )
453  nt = 1
454  IF( izero.EQ.0 ) THEN
455 *
456 * Check residual of computed solution.
457 *
458  CALL dlacpy( 'Full', n, nrhs, b, lda, work,
459  $ lda )
460  CALL dgtt02( trans, n, nrhs, a, a( m+1 ),
461  $ a( n+m+1 ), x, lda, work, lda,
462  $ result( 2 ) )
463 *
464 * Check solution from generated exact solution.
465 *
466  CALL dget04( n, nrhs, x, lda, xact, lda, rcondc,
467  $ result( 3 ) )
468  nt = 3
469  END IF
470 *
471 * Print information about the tests that did not pass
472 * the threshold.
473 *
474  DO 80 k = 2, nt
475  IF( result( k ).GE.thresh ) THEN
476  IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
477  $ CALL aladhd( nout, path )
478  WRITE( nout, fmt = 9999 )'DGTSV ', n, imat,
479  $ k, result( k )
480  nfail = nfail + 1
481  END IF
482  80 continue
483  nrun = nrun + nt - 1
484  END IF
485 *
486 * --- Test DGTSVX ---
487 *
488  IF( ifact.GT.1 ) THEN
489 *
490 * Initialize AF to zero.
491 *
492  DO 90 i = 1, 3*n - 2
493  af( i ) = zero
494  90 continue
495  END IF
496  CALL dlaset( 'Full', n, nrhs, zero, zero, x, lda )
497 *
498 * Solve the system and compute the condition number and
499 * error bounds using DGTSVX.
500 *
501  srnamt = 'DGTSVX'
502  CALL dgtsvx( fact, trans, n, nrhs, a, a( m+1 ),
503  $ a( n+m+1 ), af, af( m+1 ), af( n+m+1 ),
504  $ af( n+2*m+1 ), iwork, b, lda, x, lda,
505  $ rcond, rwork, rwork( nrhs+1 ), work,
506  $ iwork( n+1 ), info )
507 *
508 * Check the error code from DGTSVX.
509 *
510  IF( info.NE.izero )
511  $ CALL alaerh( path, 'DGTSVX', info, izero,
512  $ fact // trans, n, n, 1, 1, nrhs, imat,
513  $ nfail, nerrs, nout )
514 *
515  IF( ifact.GE.2 ) THEN
516 *
517 * Reconstruct matrix from factors and compute
518 * residual.
519 *
520  CALL dgtt01( n, a, a( m+1 ), a( n+m+1 ), af,
521  $ af( m+1 ), af( n+m+1 ), af( n+2*m+1 ),
522  $ iwork, work, lda, rwork, result( 1 ) )
523  k1 = 1
524  ELSE
525  k1 = 2
526  END IF
527 *
528  IF( info.EQ.0 ) THEN
529  trfcon = .false.
530 *
531 * Check residual of computed solution.
532 *
533  CALL dlacpy( 'Full', n, nrhs, b, lda, work, lda )
534  CALL dgtt02( trans, n, nrhs, a, a( m+1 ),
535  $ a( n+m+1 ), x, lda, work, lda,
536  $ result( 2 ) )
537 *
538 * Check solution from generated exact solution.
539 *
540  CALL dget04( n, nrhs, x, lda, xact, lda, rcondc,
541  $ result( 3 ) )
542 *
543 * Check the error bounds from iterative refinement.
544 *
545  CALL dgtt05( trans, n, nrhs, a, a( m+1 ),
546  $ a( n+m+1 ), b, lda, x, lda, xact, lda,
547  $ rwork, rwork( nrhs+1 ), result( 4 ) )
548  nt = 5
549  END IF
550 *
551 * Print information about the tests that did not pass
552 * the threshold.
553 *
554  DO 100 k = k1, nt
555  IF( result( k ).GE.thresh ) THEN
556  IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
557  $ CALL aladhd( nout, path )
558  WRITE( nout, fmt = 9998 )'DGTSVX', fact, trans,
559  $ n, imat, k, result( k )
560  nfail = nfail + 1
561  END IF
562  100 continue
563 *
564 * Check the reciprocal of the condition number.
565 *
566  result( 6 ) = dget06( rcond, rcondc )
567  IF( result( 6 ).GE.thresh ) THEN
568  IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
569  $ CALL aladhd( nout, path )
570  WRITE( nout, fmt = 9998 )'DGTSVX', fact, trans, n,
571  $ imat, k, result( k )
572  nfail = nfail + 1
573  END IF
574  nrun = nrun + nt - k1 + 2
575 *
576  110 continue
577  120 continue
578  130 continue
579  140 continue
580 *
581 * Print a summary of the results.
582 *
583  CALL alasvm( path, nout, nfail, nrun, nerrs )
584 *
585  9999 format( 1x, a, ', N =', i5, ', type ', i2, ', test ', i2,
586  $ ', ratio = ', g12.5 )
587  9998 format( 1x, a, ', FACT=''', a1, ''', TRANS=''', a1, ''', N =',
588  $ i5, ', type ', i2, ', test ', i2, ', ratio = ', g12.5 )
589  return
590 *
591 * End of DDRVGT
592 *
593  END