LAPACK  3.10.1
LAPACK: Linear Algebra PACKage
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 *> \ingroup complex_lin
136 *
137 * =====================================================================
138  SUBROUTINE cdrvpt( DOTYPE, NN, NVAL, NRHS, THRESH, TSTERR, A, D,
139  $ E, B, X, XACT, WORK, RWORK, NOUT )
140 *
141 * -- LAPACK test routine --
142 * -- LAPACK is a software package provided by Univ. of Tennessee, --
143 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
144 *
145 * .. Scalar Arguments ..
146  LOGICAL TSTERR
147  INTEGER NN, NOUT, NRHS
148  REAL THRESH
149 * ..
150 * .. Array Arguments ..
151  LOGICAL DOTYPE( * )
152  INTEGER NVAL( * )
153  REAL D( * ), RWORK( * )
154  COMPLEX A( * ), B( * ), E( * ), WORK( * ), X( * ),
155  $ xact( * )
156 * ..
157 *
158 * =====================================================================
159 *
160 * .. Parameters ..
161  REAL ONE, ZERO
162  parameter( one = 1.0e+0, zero = 0.0e+0 )
163  INTEGER NTYPES
164  parameter( ntypes = 12 )
165  INTEGER NTESTS
166  parameter( ntests = 6 )
167 * ..
168 * .. Local Scalars ..
169  LOGICAL ZEROT
170  CHARACTER DIST, FACT, TYPE
171  CHARACTER*3 PATH
172  INTEGER I, IA, IFACT, IMAT, IN, INFO, IX, IZERO, J, K,
173  $ k1, kl, ku, lda, mode, n, nerrs, nfail, nimat,
174  $ nrun, nt
175  REAL AINVNM, ANORM, COND, DMAX, RCOND, RCONDC
176 * ..
177 * .. Local Arrays ..
178  INTEGER ISEED( 4 ), ISEEDY( 4 )
179  REAL RESULT( NTESTS ), Z( 3 )
180 * ..
181 * .. External Functions ..
182  INTEGER ISAMAX
183  REAL CLANHT, SCASUM, SGET06
184  EXTERNAL isamax, clanht, scasum, sget06
185 * ..
186 * .. External Subroutines ..
187  EXTERNAL aladhd, alaerh, alasvm, ccopy, cerrvx, cget04,
191 * ..
192 * .. Intrinsic Functions ..
193  INTRINSIC abs, cmplx, max
194 * ..
195 * .. Scalars in Common ..
196  LOGICAL LERR, OK
197  CHARACTER*32 SRNAMT
198  INTEGER INFOT, NUNIT
199 * ..
200 * .. Common blocks ..
201  COMMON / infoc / infot, nunit, ok, lerr
202  COMMON / srnamc / srnamt
203 * ..
204 * .. Data statements ..
205  DATA iseedy / 0, 0, 0, 1 /
206 * ..
207 * .. Executable Statements ..
208 *
209  path( 1: 1 ) = 'Complex precision'
210  path( 2: 3 ) = 'PT'
211  nrun = 0
212  nfail = 0
213  nerrs = 0
214  DO 10 i = 1, 4
215  iseed( i ) = iseedy( i )
216  10 CONTINUE
217 *
218 * Test the error exits
219 *
220  IF( tsterr )
221  $ CALL cerrvx( path, nout )
222  infot = 0
223 *
224  DO 120 in = 1, nn
225 *
226 * Do for each value of N in NVAL.
227 *
228  n = nval( in )
229  lda = max( 1, n )
230  nimat = ntypes
231  IF( n.LE.0 )
232  $ nimat = 1
233 *
234  DO 110 imat = 1, nimat
235 *
236 * Do the tests only if DOTYPE( IMAT ) is true.
237 *
238  IF( n.GT.0 .AND. .NOT.dotype( imat ) )
239  $ GO TO 110
240 *
241 * Set up parameters with CLATB4.
242 *
243  CALL clatb4( path, imat, n, n, TYPE, kl, ku, anorm, mode,
244  $ cond, dist )
245 *
246  zerot = imat.GE.8 .AND. imat.LE.10
247  IF( imat.LE.6 ) THEN
248 *
249 * Type 1-6: generate a symmetric tridiagonal matrix of
250 * known condition number in lower triangular band storage.
251 *
252  srnamt = 'CLATMS'
253  CALL clatms( n, n, dist, iseed, TYPE, rwork, mode, cond,
254  $ anorm, kl, ku, 'B', a, 2, work, info )
255 *
256 * Check the error code from CLATMS.
257 *
258  IF( info.NE.0 ) THEN
259  CALL alaerh( path, 'CLATMS', info, 0, ' ', n, n, kl,
260  $ ku, -1, imat, nfail, nerrs, nout )
261  GO TO 110
262  END IF
263  izero = 0
264 *
265 * Copy the matrix to D and E.
266 *
267  ia = 1
268  DO 20 i = 1, n - 1
269  d( i ) = a( ia )
270  e( i ) = a( ia+1 )
271  ia = ia + 2
272  20 CONTINUE
273  IF( n.GT.0 )
274  $ d( n ) = a( ia )
275  ELSE
276 *
277 * Type 7-12: generate a diagonally dominant matrix with
278 * unknown condition number in the vectors D and E.
279 *
280  IF( .NOT.zerot .OR. .NOT.dotype( 7 ) ) THEN
281 *
282 * Let D and E have values from [-1,1].
283 *
284  CALL slarnv( 2, iseed, n, d )
285  CALL clarnv( 2, iseed, n-1, e )
286 *
287 * Make the tridiagonal matrix diagonally dominant.
288 *
289  IF( n.EQ.1 ) THEN
290  d( 1 ) = abs( d( 1 ) )
291  ELSE
292  d( 1 ) = abs( d( 1 ) ) + abs( e( 1 ) )
293  d( n ) = abs( d( n ) ) + abs( e( n-1 ) )
294  DO 30 i = 2, n - 1
295  d( i ) = abs( d( i ) ) + abs( e( i ) ) +
296  $ abs( e( i-1 ) )
297  30 CONTINUE
298  END IF
299 *
300 * Scale D and E so the maximum element is ANORM.
301 *
302  ix = isamax( n, d, 1 )
303  dmax = d( ix )
304  CALL sscal( n, anorm / dmax, d, 1 )
305  IF( n.GT.1 )
306  $ CALL csscal( n-1, anorm / dmax, e, 1 )
307 *
308  ELSE IF( izero.GT.0 ) THEN
309 *
310 * Reuse the last matrix by copying back the zeroed out
311 * elements.
312 *
313  IF( izero.EQ.1 ) THEN
314  d( 1 ) = z( 2 )
315  IF( n.GT.1 )
316  $ e( 1 ) = z( 3 )
317  ELSE IF( izero.EQ.n ) THEN
318  e( n-1 ) = z( 1 )
319  d( n ) = z( 2 )
320  ELSE
321  e( izero-1 ) = z( 1 )
322  d( izero ) = z( 2 )
323  e( izero ) = z( 3 )
324  END IF
325  END IF
326 *
327 * For types 8-10, set one row and column of the matrix to
328 * zero.
329 *
330  izero = 0
331  IF( imat.EQ.8 ) THEN
332  izero = 1
333  z( 2 ) = d( 1 )
334  d( 1 ) = zero
335  IF( n.GT.1 ) THEN
336  z( 3 ) = e( 1 )
337  e( 1 ) = zero
338  END IF
339  ELSE IF( imat.EQ.9 ) THEN
340  izero = n
341  IF( n.GT.1 ) THEN
342  z( 1 ) = e( n-1 )
343  e( n-1 ) = zero
344  END IF
345  z( 2 ) = d( n )
346  d( n ) = zero
347  ELSE IF( imat.EQ.10 ) THEN
348  izero = ( n+1 ) / 2
349  IF( izero.GT.1 ) THEN
350  z( 1 ) = e( izero-1 )
351  e( izero-1 ) = zero
352  z( 3 ) = e( izero )
353  e( izero ) = zero
354  END IF
355  z( 2 ) = d( izero )
356  d( izero ) = zero
357  END IF
358  END IF
359 *
360 * Generate NRHS random solution vectors.
361 *
362  ix = 1
363  DO 40 j = 1, nrhs
364  CALL clarnv( 2, iseed, n, xact( ix ) )
365  ix = ix + lda
366  40 CONTINUE
367 *
368 * Set the right hand side.
369 *
370  CALL claptm( 'Lower', n, nrhs, one, d, e, xact, lda, zero,
371  $ b, lda )
372 *
373  DO 100 ifact = 1, 2
374  IF( ifact.EQ.1 ) THEN
375  fact = 'F'
376  ELSE
377  fact = 'N'
378  END IF
379 *
380 * Compute the condition number for comparison with
381 * the value returned by CPTSVX.
382 *
383  IF( zerot ) THEN
384  IF( ifact.EQ.1 )
385  $ GO TO 100
386  rcondc = zero
387 *
388  ELSE IF( ifact.EQ.1 ) THEN
389 *
390 * Compute the 1-norm of A.
391 *
392  anorm = clanht( '1', n, d, e )
393 *
394  CALL scopy( n, d, 1, d( n+1 ), 1 )
395  IF( n.GT.1 )
396  $ CALL ccopy( n-1, e, 1, e( n+1 ), 1 )
397 *
398 * Factor the matrix A.
399 *
400  CALL cpttrf( n, d( n+1 ), e( n+1 ), info )
401 *
402 * Use CPTTRS to solve for one column at a time of
403 * inv(A), computing the maximum column sum as we go.
404 *
405  ainvnm = zero
406  DO 60 i = 1, n
407  DO 50 j = 1, n
408  x( j ) = zero
409  50 CONTINUE
410  x( i ) = one
411  CALL cpttrs( 'Lower', n, 1, d( n+1 ), e( n+1 ), x,
412  $ lda, info )
413  ainvnm = max( ainvnm, scasum( n, x, 1 ) )
414  60 CONTINUE
415 *
416 * Compute the 1-norm condition number of A.
417 *
418  IF( anorm.LE.zero .OR. ainvnm.LE.zero ) THEN
419  rcondc = one
420  ELSE
421  rcondc = ( one / anorm ) / ainvnm
422  END IF
423  END IF
424 *
425  IF( ifact.EQ.2 ) THEN
426 *
427 * --- Test CPTSV --
428 *
429  CALL scopy( n, d, 1, d( n+1 ), 1 )
430  IF( n.GT.1 )
431  $ CALL ccopy( n-1, e, 1, e( n+1 ), 1 )
432  CALL clacpy( 'Full', n, nrhs, b, lda, x, lda )
433 *
434 * Factor A as L*D*L' and solve the system A*X = B.
435 *
436  srnamt = 'CPTSV '
437  CALL cptsv( n, nrhs, d( n+1 ), e( n+1 ), x, lda,
438  $ info )
439 *
440 * Check error code from CPTSV .
441 *
442  IF( info.NE.izero )
443  $ CALL alaerh( path, 'CPTSV ', info, izero, ' ', n,
444  $ n, 1, 1, nrhs, imat, nfail, nerrs,
445  $ nout )
446  nt = 0
447  IF( izero.EQ.0 ) THEN
448 *
449 * Check the factorization by computing the ratio
450 * norm(L*D*L' - A) / (n * norm(A) * EPS )
451 *
452  CALL cptt01( n, d, e, d( n+1 ), e( n+1 ), work,
453  $ result( 1 ) )
454 *
455 * Compute the residual in the solution.
456 *
457  CALL clacpy( 'Full', n, nrhs, b, lda, work, lda )
458  CALL cptt02( 'Lower', n, nrhs, d, e, x, lda, work,
459  $ lda, result( 2 ) )
460 *
461 * Check solution from generated exact solution.
462 *
463  CALL cget04( n, nrhs, x, lda, xact, lda, rcondc,
464  $ result( 3 ) )
465  nt = 3
466  END IF
467 *
468 * Print information about the tests that did not pass
469 * the threshold.
470 *
471  DO 70 k = 1, nt
472  IF( result( k ).GE.thresh ) THEN
473  IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
474  $ CALL aladhd( nout, path )
475  WRITE( nout, fmt = 9999 )'CPTSV ', n, imat, k,
476  $ result( k )
477  nfail = nfail + 1
478  END IF
479  70 CONTINUE
480  nrun = nrun + nt
481  END IF
482 *
483 * --- Test CPTSVX ---
484 *
485  IF( ifact.GT.1 ) THEN
486 *
487 * Initialize D( N+1:2*N ) and E( N+1:2*N ) to zero.
488 *
489  DO 80 i = 1, n - 1
490  d( n+i ) = zero
491  e( n+i ) = zero
492  80 CONTINUE
493  IF( n.GT.0 )
494  $ d( n+n ) = zero
495  END IF
496 *
497  CALL claset( 'Full', n, nrhs, cmplx( zero ),
498  $ cmplx( zero ), x, lda )
499 *
500 * Solve the system and compute the condition number and
501 * error bounds using CPTSVX.
502 *
503  srnamt = 'CPTSVX'
504  CALL cptsvx( fact, n, nrhs, d, e, d( n+1 ), e( n+1 ), b,
505  $ lda, x, lda, rcond, rwork, rwork( nrhs+1 ),
506  $ work, rwork( 2*nrhs+1 ), info )
507 *
508 * Check the error code from CPTSVX.
509 *
510  IF( info.NE.izero )
511  $ CALL alaerh( path, 'CPTSVX', 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 cptt01( 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 clacpy( 'Full', n, nrhs, b, lda, work, lda )
529  CALL cptt02( 'Lower', n, nrhs, d, e, x, lda, work,
530  $ lda, result( 2 ) )
531 *
532 * Check solution from generated exact solution.
533 *
534  CALL cget04( n, nrhs, x, lda, xact, lda, rcondc,
535  $ result( 3 ) )
536 *
537 * Check error bounds from iterative refinement.
538 *
539  CALL cptt05( 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 ) = sget06( 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 )'CPTSVX', 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 CDRVPT
577 *
578  END
subroutine slarnv(IDIST, ISEED, N, X)
SLARNV returns a vector of random numbers from a uniform or normal distribution.
Definition: slarnv.f:97
subroutine alasvm(TYPE, NOUT, NFAIL, NRUN, NERRS)
ALASVM
Definition: alasvm.f:73
subroutine aladhd(IOUNIT, PATH)
ALADHD
Definition: aladhd.f:90
subroutine alaerh(PATH, SUBNAM, INFO, INFOE, OPTS, M, N, KL, KU, N5, IMAT, NFAIL, NERRS, NOUT)
ALAERH
Definition: alaerh.f:147
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
Definition: ccopy.f:81
subroutine csscal(N, SA, CX, INCX)
CSSCAL
Definition: csscal.f:78
subroutine claptm(UPLO, N, NRHS, ALPHA, D, E, X, LDX, BETA, B, LDB)
CLAPTM
Definition: claptm.f:129
subroutine clatb4(PATH, IMAT, M, N, TYPE, KL, KU, ANORM, MODE, CNDNUM, DIST)
CLATB4
Definition: clatb4.f:121
subroutine cget04(N, NRHS, X, LDX, XACT, LDXACT, RCOND, RESID)
CGET04
Definition: cget04.f:102
subroutine cptt05(N, NRHS, D, E, B, LDB, X, LDX, XACT, LDXACT, FERR, BERR, RESLTS)
CPTT05
Definition: cptt05.f:150
subroutine cdrvpt(DOTYPE, NN, NVAL, NRHS, THRESH, TSTERR, A, D, E, B, X, XACT, WORK, RWORK, NOUT)
CDRVPT
Definition: cdrvpt.f:140
subroutine cerrvx(PATH, NUNIT)
CERRVX
Definition: cerrvx.f:55
subroutine cptt02(UPLO, N, NRHS, D, E, X, LDX, B, LDB, RESID)
CPTT02
Definition: cptt02.f:115
subroutine cptt01(N, D, E, DF, EF, WORK, RESID)
CPTT01
Definition: cptt01.f:92
subroutine clatms(M, N, DIST, ISEED, SYM, D, MODE, COND, DMAX, KL, KU, PACK, A, LDA, WORK, INFO)
CLATMS
Definition: clatms.f:332
subroutine claset(UPLO, M, N, ALPHA, BETA, A, LDA)
CLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition: claset.f:106
subroutine clarnv(IDIST, ISEED, N, X)
CLARNV returns a vector of random numbers from a uniform or normal distribution.
Definition: clarnv.f:99
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
Definition: clacpy.f:103
subroutine cpttrs(UPLO, N, NRHS, D, E, B, LDB, INFO)
CPTTRS
Definition: cpttrs.f:121
subroutine cpttrf(N, D, E, INFO)
CPTTRF
Definition: cpttrf.f:92
subroutine cptsvx(FACT, N, NRHS, D, E, DF, EF, B, LDB, X, LDX, RCOND, FERR, BERR, WORK, RWORK, INFO)
CPTSVX computes the solution to system of linear equations A * X = B for PT matrices
Definition: cptsvx.f:234
subroutine cptsv(N, NRHS, D, E, B, LDB, INFO)
CPTSV computes the solution to system of linear equations A * X = B for PT matrices
Definition: cptsv.f:115
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:82
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:79