LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
cdrvpt.f
Go to the documentation of this file.
1 *> \brief \b CDRVPT
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 CDRVPT( 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 * REAL THRESH
18 * ..
19 * .. Array Arguments ..
20 * LOGICAL DOTYPE( * )
21 * INTEGER NVAL( * )
22 * REAL D( * ), RWORK( * )
23 * COMPLEX A( * ), B( * ), E( * ), WORK( * ), X( * ),
24 * $ XACT( * )
25 * ..
26 *
27 *
28 *> \par Purpose:
29 * =============
30 *>
31 *> \verbatim
32 *>
33 *> CDRVPT tests CPTSV and -SVX.
34 *> \endverbatim
35 *
36 * Arguments:
37 * ==========
38 *
39 *> \param[in] DOTYPE
40 *> \verbatim
41 *> DOTYPE is LOGICAL array, dimension (NTYPES)
42 *> The matrix types to be used for testing. Matrices of type j
43 *> (for 1 <= j <= NTYPES) are used for testing if DOTYPE(j) =
44 *> .TRUE.; if DOTYPE(j) = .FALSE., then type j is not used.
45 *> \endverbatim
46 *>
47 *> \param[in] NN
48 *> \verbatim
49 *> NN is INTEGER
50 *> The number of values of N contained in the vector NVAL.
51 *> \endverbatim
52 *>
53 *> \param[in] NVAL
54 *> \verbatim
55 *> NVAL is INTEGER array, dimension (NN)
56 *> The values of the matrix dimension N.
57 *> \endverbatim
58 *>
59 *> \param[in] NRHS
60 *> \verbatim
61 *> NRHS is INTEGER
62 *> The number of right hand side vectors to be generated for
63 *> each linear system.
64 *> \endverbatim
65 *>
66 *> \param[in] THRESH
67 *> \verbatim
68 *> THRESH is REAL
69 *> The threshold value for the test ratios. A result is
70 *> included in the output file if RESULT >= THRESH. To have
71 *> every test ratio printed, use THRESH = 0.
72 *> \endverbatim
73 *>
74 *> \param[in] TSTERR
75 *> \verbatim
76 *> TSTERR is LOGICAL
77 *> Flag that indicates whether error exits are to be tested.
78 *> \endverbatim
79 *>
80 *> \param[out] A
81 *> \verbatim
82 *> A is COMPLEX array, dimension (NMAX*2)
83 *> \endverbatim
84 *>
85 *> \param[out] D
86 *> \verbatim
87 *> D is REAL array, dimension (NMAX*2)
88 *> \endverbatim
89 *>
90 *> \param[out] E
91 *> \verbatim
92 *> E is COMPLEX array, dimension (NMAX*2)
93 *> \endverbatim
94 *>
95 *> \param[out] B
96 *> \verbatim
97 *> B is COMPLEX array, dimension (NMAX*NRHS)
98 *> \endverbatim
99 *>
100 *> \param[out] X
101 *> \verbatim
102 *> X is COMPLEX array, dimension (NMAX*NRHS)
103 *> \endverbatim
104 *>
105 *> \param[out] XACT
106 *> \verbatim
107 *> XACT is COMPLEX array, dimension (NMAX*NRHS)
108 *> \endverbatim
109 *>
110 *> \param[out] WORK
111 *> \verbatim
112 *> WORK is COMPLEX array, dimension
113 *> (NMAX*max(3,NRHS))
114 *> \endverbatim
115 *>
116 *> \param[out] RWORK
117 *> \verbatim
118 *> RWORK is REAL array, dimension (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 complex_lin
138 *
139 * =====================================================================
140  SUBROUTINE cdrvpt( 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  REAL thresh
152 * ..
153 * .. Array Arguments ..
154  LOGICAL dotype( * )
155  INTEGER nval( * )
156  REAL d( * ), rwork( * )
157  COMPLEX a( * ), b( * ), e( * ), work( * ), x( * ),
158  $ xact( * )
159 * ..
160 *
161 * =====================================================================
162 *
163 * .. Parameters ..
164  REAL one, zero
165  parameter( one = 1.0e+0, zero = 0.0e+0 )
166  INTEGER ntypes
167  parameter( ntypes = 12 )
168  INTEGER ntests
169  parameter( ntests = 6 )
170 * ..
171 * .. Local Scalars ..
172  LOGICAL zerot
173  CHARACTER dist, fact, type
174  CHARACTER*3 path
175  INTEGER i, ia, ifact, imat, in, info, ix, izero, j, k,
176  $ k1, kl, ku, lda, mode, n, nerrs, nfail, nimat,
177  $ nrun, nt
178  REAL ainvnm, anorm, cond, dmax, rcond, rcondc
179 * ..
180 * .. Local Arrays ..
181  INTEGER iseed( 4 ), iseedy( 4 )
182  REAL result( ntests ), z( 3 )
183 * ..
184 * .. External Functions ..
185  INTEGER isamax
186  REAL clanht, scasum, sget06
187  EXTERNAL isamax, clanht, scasum, sget06
188 * ..
189 * .. External Subroutines ..
190  EXTERNAL aladhd, alaerh, alasvm, ccopy, cerrvx, cget04,
191  $ clacpy, claptm, clarnv, claset, clatb4, clatms,
192  $ cptsv, cptsvx, cptt01, cptt02, cptt05, cpttrf,
193  $ cpttrs, csscal, scopy, slarnv, sscal
194 * ..
195 * .. Intrinsic Functions ..
196  INTRINSIC abs, cmplx, max
197 * ..
198 * .. Scalars in Common ..
199  LOGICAL lerr, ok
200  CHARACTER*32 srnamt
201  INTEGER infot, nunit
202 * ..
203 * .. Common blocks ..
204  common / infoc / infot, nunit, ok, lerr
205  common / srnamc / srnamt
206 * ..
207 * .. Data statements ..
208  DATA iseedy / 0, 0, 0, 1 /
209 * ..
210 * .. Executable Statements ..
211 *
212  path( 1: 1 ) = 'Complex precision'
213  path( 2: 3 ) = 'PT'
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 cerrvx( path, nout )
225  infot = 0
226 *
227  DO 120 in = 1, nn
228 *
229 * Do for each value of N in NVAL.
230 *
231  n = nval( in )
232  lda = max( 1, n )
233  nimat = ntypes
234  IF( n.LE.0 )
235  $ nimat = 1
236 *
237  DO 110 imat = 1, nimat
238 *
239 * Do the tests only if DOTYPE( IMAT ) is true.
240 *
241  IF( n.GT.0 .AND. .NOT.dotype( imat ) )
242  $ go to 110
243 *
244 * Set up parameters with CLATB4.
245 *
246  CALL clatb4( path, imat, n, n, type, kl, ku, anorm, mode,
247  $ cond, dist )
248 *
249  zerot = imat.GE.8 .AND. imat.LE.10
250  IF( imat.LE.6 ) THEN
251 *
252 * Type 1-6: generate a symmetric tridiagonal matrix of
253 * known condition number in lower triangular band storage.
254 *
255  srnamt = 'CLATMS'
256  CALL clatms( n, n, dist, iseed, type, rwork, mode, cond,
257  $ anorm, kl, ku, 'B', a, 2, work, info )
258 *
259 * Check the error code from CLATMS.
260 *
261  IF( info.NE.0 ) THEN
262  CALL alaerh( path, 'CLATMS', info, 0, ' ', n, n, kl,
263  $ ku, -1, imat, nfail, nerrs, nout )
264  go to 110
265  END IF
266  izero = 0
267 *
268 * Copy the matrix to D and E.
269 *
270  ia = 1
271  DO 20 i = 1, n - 1
272  d( i ) = a( ia )
273  e( i ) = a( ia+1 )
274  ia = ia + 2
275  20 continue
276  IF( n.GT.0 )
277  $ d( n ) = a( ia )
278  ELSE
279 *
280 * Type 7-12: generate a diagonally dominant matrix with
281 * unknown condition number in the vectors D and E.
282 *
283  IF( .NOT.zerot .OR. .NOT.dotype( 7 ) ) THEN
284 *
285 * Let D and E have values from [-1,1].
286 *
287  CALL slarnv( 2, iseed, n, d )
288  CALL clarnv( 2, iseed, n-1, e )
289 *
290 * Make the tridiagonal matrix diagonally dominant.
291 *
292  IF( n.EQ.1 ) THEN
293  d( 1 ) = abs( d( 1 ) )
294  ELSE
295  d( 1 ) = abs( d( 1 ) ) + abs( e( 1 ) )
296  d( n ) = abs( d( n ) ) + abs( e( n-1 ) )
297  DO 30 i = 2, n - 1
298  d( i ) = abs( d( i ) ) + abs( e( i ) ) +
299  $ abs( e( i-1 ) )
300  30 continue
301  END IF
302 *
303 * Scale D and E so the maximum element is ANORM.
304 *
305  ix = isamax( n, d, 1 )
306  dmax = d( ix )
307  CALL sscal( n, anorm / dmax, d, 1 )
308  IF( n.GT.1 )
309  $ CALL csscal( n-1, anorm / dmax, e, 1 )
310 *
311  ELSE IF( izero.GT.0 ) THEN
312 *
313 * Reuse the last matrix by copying back the zeroed out
314 * elements.
315 *
316  IF( izero.EQ.1 ) THEN
317  d( 1 ) = z( 2 )
318  IF( n.GT.1 )
319  $ e( 1 ) = z( 3 )
320  ELSE IF( izero.EQ.n ) THEN
321  e( n-1 ) = z( 1 )
322  d( n ) = z( 2 )
323  ELSE
324  e( izero-1 ) = z( 1 )
325  d( izero ) = z( 2 )
326  e( izero ) = z( 3 )
327  END IF
328  END IF
329 *
330 * For types 8-10, set one row and column of the matrix to
331 * zero.
332 *
333  izero = 0
334  IF( imat.EQ.8 ) THEN
335  izero = 1
336  z( 2 ) = d( 1 )
337  d( 1 ) = zero
338  IF( n.GT.1 ) THEN
339  z( 3 ) = e( 1 )
340  e( 1 ) = zero
341  END IF
342  ELSE IF( imat.EQ.9 ) THEN
343  izero = n
344  IF( n.GT.1 ) THEN
345  z( 1 ) = e( n-1 )
346  e( n-1 ) = zero
347  END IF
348  z( 2 ) = d( n )
349  d( n ) = zero
350  ELSE IF( imat.EQ.10 ) THEN
351  izero = ( n+1 ) / 2
352  IF( izero.GT.1 ) THEN
353  z( 1 ) = e( izero-1 )
354  e( izero-1 ) = zero
355  z( 3 ) = e( izero )
356  e( izero ) = zero
357  END IF
358  z( 2 ) = d( izero )
359  d( izero ) = zero
360  END IF
361  END IF
362 *
363 * Generate NRHS random solution vectors.
364 *
365  ix = 1
366  DO 40 j = 1, nrhs
367  CALL clarnv( 2, iseed, n, xact( ix ) )
368  ix = ix + lda
369  40 continue
370 *
371 * Set the right hand side.
372 *
373  CALL claptm( 'Lower', n, nrhs, one, d, e, xact, lda, zero,
374  $ b, lda )
375 *
376  DO 100 ifact = 1, 2
377  IF( ifact.EQ.1 ) THEN
378  fact = 'F'
379  ELSE
380  fact = 'N'
381  END IF
382 *
383 * Compute the condition number for comparison with
384 * the value returned by CPTSVX.
385 *
386  IF( zerot ) THEN
387  IF( ifact.EQ.1 )
388  $ go to 100
389  rcondc = zero
390 *
391  ELSE IF( ifact.EQ.1 ) THEN
392 *
393 * Compute the 1-norm of A.
394 *
395  anorm = clanht( '1', n, d, e )
396 *
397  CALL scopy( n, d, 1, d( n+1 ), 1 )
398  IF( n.GT.1 )
399  $ CALL ccopy( n-1, e, 1, e( n+1 ), 1 )
400 *
401 * Factor the matrix A.
402 *
403  CALL cpttrf( n, d( n+1 ), e( n+1 ), info )
404 *
405 * Use CPTTRS to solve for one column at a time of
406 * inv(A), computing the maximum column sum as we go.
407 *
408  ainvnm = zero
409  DO 60 i = 1, n
410  DO 50 j = 1, n
411  x( j ) = zero
412  50 continue
413  x( i ) = one
414  CALL cpttrs( 'Lower', n, 1, d( n+1 ), e( n+1 ), x,
415  $ lda, info )
416  ainvnm = max( ainvnm, scasum( n, x, 1 ) )
417  60 continue
418 *
419 * Compute the 1-norm condition number of A.
420 *
421  IF( anorm.LE.zero .OR. ainvnm.LE.zero ) THEN
422  rcondc = one
423  ELSE
424  rcondc = ( one / anorm ) / ainvnm
425  END IF
426  END IF
427 *
428  IF( ifact.EQ.2 ) THEN
429 *
430 * --- Test CPTSV --
431 *
432  CALL scopy( n, d, 1, d( n+1 ), 1 )
433  IF( n.GT.1 )
434  $ CALL ccopy( n-1, e, 1, e( n+1 ), 1 )
435  CALL clacpy( 'Full', n, nrhs, b, lda, x, lda )
436 *
437 * Factor A as L*D*L' and solve the system A*X = B.
438 *
439  srnamt = 'CPTSV '
440  CALL cptsv( n, nrhs, d( n+1 ), e( n+1 ), x, lda,
441  $ info )
442 *
443 * Check error code from CPTSV .
444 *
445  IF( info.NE.izero )
446  $ CALL alaerh( path, 'CPTSV ', info, izero, ' ', n,
447  $ n, 1, 1, nrhs, imat, nfail, nerrs,
448  $ nout )
449  nt = 0
450  IF( izero.EQ.0 ) THEN
451 *
452 * Check the factorization by computing the ratio
453 * norm(L*D*L' - A) / (n * norm(A) * EPS )
454 *
455  CALL cptt01( n, d, e, d( n+1 ), e( n+1 ), work,
456  $ result( 1 ) )
457 *
458 * Compute the residual in the solution.
459 *
460  CALL clacpy( 'Full', n, nrhs, b, lda, work, lda )
461  CALL cptt02( 'Lower', n, nrhs, d, e, x, lda, work,
462  $ lda, result( 2 ) )
463 *
464 * Check solution from generated exact solution.
465 *
466  CALL cget04( 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 70 k = 1, 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 )'CPTSV ', n, imat, k,
479  $ result( k )
480  nfail = nfail + 1
481  END IF
482  70 continue
483  nrun = nrun + nt
484  END IF
485 *
486 * --- Test CPTSVX ---
487 *
488  IF( ifact.GT.1 ) THEN
489 *
490 * Initialize D( N+1:2*N ) and E( N+1:2*N ) to zero.
491 *
492  DO 80 i = 1, n - 1
493  d( n+i ) = zero
494  e( n+i ) = zero
495  80 continue
496  IF( n.GT.0 )
497  $ d( n+n ) = zero
498  END IF
499 *
500  CALL claset( 'Full', n, nrhs, cmplx( zero ),
501  $ cmplx( zero ), x, lda )
502 *
503 * Solve the system and compute the condition number and
504 * error bounds using CPTSVX.
505 *
506  srnamt = 'CPTSVX'
507  CALL cptsvx( fact, n, nrhs, d, e, d( n+1 ), e( n+1 ), b,
508  $ lda, x, lda, rcond, rwork, rwork( nrhs+1 ),
509  $ work, rwork( 2*nrhs+1 ), info )
510 *
511 * Check the error code from CPTSVX.
512 *
513  IF( info.NE.izero )
514  $ CALL alaerh( path, 'CPTSVX', info, izero, fact, n, n,
515  $ 1, 1, nrhs, imat, nfail, nerrs, nout )
516  IF( izero.EQ.0 ) THEN
517  IF( ifact.EQ.2 ) THEN
518 *
519 * Check the factorization by computing the ratio
520 * norm(L*D*L' - A) / (n * norm(A) * EPS )
521 *
522  k1 = 1
523  CALL cptt01( n, d, e, d( n+1 ), e( n+1 ), work,
524  $ result( 1 ) )
525  ELSE
526  k1 = 2
527  END IF
528 *
529 * Compute the residual in the solution.
530 *
531  CALL clacpy( 'Full', n, nrhs, b, lda, work, lda )
532  CALL cptt02( 'Lower', n, nrhs, d, e, x, lda, work,
533  $ lda, result( 2 ) )
534 *
535 * Check solution from generated exact solution.
536 *
537  CALL cget04( n, nrhs, x, lda, xact, lda, rcondc,
538  $ result( 3 ) )
539 *
540 * Check error bounds from iterative refinement.
541 *
542  CALL cptt05( n, nrhs, d, e, b, lda, x, lda, xact, lda,
543  $ rwork, rwork( nrhs+1 ), result( 4 ) )
544  ELSE
545  k1 = 6
546  END IF
547 *
548 * Check the reciprocal of the condition number.
549 *
550  result( 6 ) = sget06( rcond, rcondc )
551 *
552 * Print information about the tests that did not pass
553 * the threshold.
554 *
555  DO 90 k = k1, 6
556  IF( result( k ).GE.thresh ) THEN
557  IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
558  $ CALL aladhd( nout, path )
559  WRITE( nout, fmt = 9998 )'CPTSVX', fact, n, imat,
560  $ k, result( k )
561  nfail = nfail + 1
562  END IF
563  90 continue
564  nrun = nrun + 7 - k1
565  100 continue
566  110 continue
567  120 continue
568 *
569 * Print a summary of the results.
570 *
571  CALL alasvm( path, nout, nfail, nrun, nerrs )
572 *
573  9999 format( 1x, a, ', N =', i5, ', type ', i2, ', test ', i2,
574  $ ', ratio = ', g12.5 )
575  9998 format( 1x, a, ', FACT=''', a1, ''', N =', i5, ', type ', i2,
576  $ ', test ', i2, ', ratio = ', g12.5 )
577  return
578 *
579 * End of CDRVPT
580 *
581  END