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