LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
sdrvpt.f
Go to the documentation of this file.
1 *> \brief \b SDRVPT
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 SDRVPT( 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 A( * ), B( * ), D( * ), E( * ), RWORK( * ),
23 * $ WORK( * ), X( * ), XACT( * )
24 * ..
25 *
26 *
27 *> \par Purpose:
28 * =============
29 *>
30 *> \verbatim
31 *>
32 *> SDRVPT tests SPTSV 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 REAL
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 REAL array, dimension (NMAX*2)
82 *> \endverbatim
83 *>
84 *> \param[out] D
85 *> \verbatim
86 *> D is REAL array, dimension (NMAX*2)
87 *> \endverbatim
88 *>
89 *> \param[out] E
90 *> \verbatim
91 *> E is REAL array, dimension (NMAX*2)
92 *> \endverbatim
93 *>
94 *> \param[out] B
95 *> \verbatim
96 *> B is REAL array, dimension (NMAX*NRHS)
97 *> \endverbatim
98 *>
99 *> \param[out] X
100 *> \verbatim
101 *> X is REAL array, dimension (NMAX*NRHS)
102 *> \endverbatim
103 *>
104 *> \param[out] XACT
105 *> \verbatim
106 *> XACT is REAL array, dimension (NMAX*NRHS)
107 *> \endverbatim
108 *>
109 *> \param[out] WORK
110 *> \verbatim
111 *> WORK is REAL array, dimension
112 *> (NMAX*max(3,NRHS))
113 *> \endverbatim
114 *>
115 *> \param[out] RWORK
116 *> \verbatim
117 *> RWORK is REAL 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 *> \ingroup single_lin
136 *
137 * =====================================================================
138  SUBROUTINE sdrvpt( 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 A( * ), B( * ), D( * ), E( * ), RWORK( * ),
154  $ work( * ), x( * ), xact( * )
155 * ..
156 *
157 * =====================================================================
158 *
159 * .. Parameters ..
160  REAL ONE, ZERO
161  parameter( one = 1.0e+0, zero = 0.0e+0 )
162  INTEGER NTYPES
163  parameter( ntypes = 12 )
164  INTEGER NTESTS
165  parameter( ntests = 6 )
166 * ..
167 * .. Local Scalars ..
168  LOGICAL ZEROT
169  CHARACTER DIST, FACT, TYPE
170  CHARACTER*3 PATH
171  INTEGER I, IA, IFACT, IMAT, IN, INFO, IX, IZERO, J, K,
172  $ k1, kl, ku, lda, mode, n, nerrs, nfail, nimat,
173  $ nrun, nt
174  REAL AINVNM, ANORM, COND, DMAX, RCOND, RCONDC
175 * ..
176 * .. Local Arrays ..
177  INTEGER ISEED( 4 ), ISEEDY( 4 )
178  REAL RESULT( NTESTS ), Z( 3 )
179 * ..
180 * .. External Functions ..
181  INTEGER ISAMAX
182  REAL SASUM, SGET06, SLANST
183  EXTERNAL isamax, sasum, sget06, slanst
184 * ..
185 * .. External Subroutines ..
186  EXTERNAL aladhd, alaerh, alasvm, scopy, serrvx, sget04,
189  $ spttrs, sscal
190 * ..
191 * .. Intrinsic Functions ..
192  INTRINSIC abs, max
193 * ..
194 * .. Scalars in Common ..
195  LOGICAL LERR, OK
196  CHARACTER*32 SRNAMT
197  INTEGER INFOT, NUNIT
198 * ..
199 * .. Common blocks ..
200  COMMON / infoc / infot, nunit, ok, lerr
201  COMMON / srnamc / srnamt
202 * ..
203 * .. Data statements ..
204  DATA iseedy / 0, 0, 0, 1 /
205 * ..
206 * .. Executable Statements ..
207 *
208  path( 1: 1 ) = 'Single precision'
209  path( 2: 3 ) = 'PT'
210  nrun = 0
211  nfail = 0
212  nerrs = 0
213  DO 10 i = 1, 4
214  iseed( i ) = iseedy( i )
215  10 CONTINUE
216 *
217 * Test the error exits
218 *
219  IF( tsterr )
220  $ CALL serrvx( path, nout )
221  infot = 0
222 *
223  DO 120 in = 1, nn
224 *
225 * Do for each value of N in NVAL.
226 *
227  n = nval( in )
228  lda = max( 1, n )
229  nimat = ntypes
230  IF( n.LE.0 )
231  $ nimat = 1
232 *
233  DO 110 imat = 1, nimat
234 *
235 * Do the tests only if DOTYPE( IMAT ) is true.
236 *
237  IF( n.GT.0 .AND. .NOT.dotype( imat ) )
238  $ GO TO 110
239 *
240 * Set up parameters with SLATB4.
241 *
242  CALL slatb4( path, imat, n, n, TYPE, kl, ku, anorm, mode,
243  $ cond, dist )
244 *
245  zerot = imat.GE.8 .AND. imat.LE.10
246  IF( imat.LE.6 ) THEN
247 *
248 * Type 1-6: generate a symmetric tridiagonal matrix of
249 * known condition number in lower triangular band storage.
250 *
251  srnamt = 'SLATMS'
252  CALL slatms( n, n, dist, iseed, TYPE, rwork, mode, cond,
253  $ anorm, kl, ku, 'B', a, 2, work, info )
254 *
255 * Check the error code from SLATMS.
256 *
257  IF( info.NE.0 ) THEN
258  CALL alaerh( path, 'SLATMS', info, 0, ' ', n, n, kl,
259  $ ku, -1, imat, nfail, nerrs, nout )
260  GO TO 110
261  END IF
262  izero = 0
263 *
264 * Copy the matrix to D and E.
265 *
266  ia = 1
267  DO 20 i = 1, n - 1
268  d( i ) = a( ia )
269  e( i ) = a( ia+1 )
270  ia = ia + 2
271  20 CONTINUE
272  IF( n.GT.0 )
273  $ d( n ) = a( ia )
274  ELSE
275 *
276 * Type 7-12: generate a diagonally dominant matrix with
277 * unknown condition number in the vectors D and E.
278 *
279  IF( .NOT.zerot .OR. .NOT.dotype( 7 ) ) THEN
280 *
281 * Let D and E have values from [-1,1].
282 *
283  CALL slarnv( 2, iseed, n, d )
284  CALL slarnv( 2, iseed, n-1, e )
285 *
286 * Make the tridiagonal matrix diagonally dominant.
287 *
288  IF( n.EQ.1 ) THEN
289  d( 1 ) = abs( d( 1 ) )
290  ELSE
291  d( 1 ) = abs( d( 1 ) ) + abs( e( 1 ) )
292  d( n ) = abs( d( n ) ) + abs( e( n-1 ) )
293  DO 30 i = 2, n - 1
294  d( i ) = abs( d( i ) ) + abs( e( i ) ) +
295  $ abs( e( i-1 ) )
296  30 CONTINUE
297  END IF
298 *
299 * Scale D and E so the maximum element is ANORM.
300 *
301  ix = isamax( n, d, 1 )
302  dmax = d( ix )
303  CALL sscal( n, anorm / dmax, d, 1 )
304  IF( n.GT.1 )
305  $ CALL sscal( n-1, anorm / dmax, e, 1 )
306 *
307  ELSE IF( izero.GT.0 ) THEN
308 *
309 * Reuse the last matrix by copying back the zeroed out
310 * elements.
311 *
312  IF( izero.EQ.1 ) THEN
313  d( 1 ) = z( 2 )
314  IF( n.GT.1 )
315  $ e( 1 ) = z( 3 )
316  ELSE IF( izero.EQ.n ) THEN
317  e( n-1 ) = z( 1 )
318  d( n ) = z( 2 )
319  ELSE
320  e( izero-1 ) = z( 1 )
321  d( izero ) = z( 2 )
322  e( izero ) = z( 3 )
323  END IF
324  END IF
325 *
326 * For types 8-10, set one row and column of the matrix to
327 * zero.
328 *
329  izero = 0
330  IF( imat.EQ.8 ) THEN
331  izero = 1
332  z( 2 ) = d( 1 )
333  d( 1 ) = zero
334  IF( n.GT.1 ) THEN
335  z( 3 ) = e( 1 )
336  e( 1 ) = zero
337  END IF
338  ELSE IF( imat.EQ.9 ) THEN
339  izero = n
340  IF( n.GT.1 ) THEN
341  z( 1 ) = e( n-1 )
342  e( n-1 ) = zero
343  END IF
344  z( 2 ) = d( n )
345  d( n ) = zero
346  ELSE IF( imat.EQ.10 ) THEN
347  izero = ( n+1 ) / 2
348  IF( izero.GT.1 ) THEN
349  z( 1 ) = e( izero-1 )
350  z( 3 ) = e( izero )
351  e( izero-1 ) = zero
352  e( izero ) = zero
353  END IF
354  z( 2 ) = d( izero )
355  d( izero ) = zero
356  END IF
357  END IF
358 *
359 * Generate NRHS random solution vectors.
360 *
361  ix = 1
362  DO 40 j = 1, nrhs
363  CALL slarnv( 2, iseed, n, xact( ix ) )
364  ix = ix + lda
365  40 CONTINUE
366 *
367 * Set the right hand side.
368 *
369  CALL slaptm( n, nrhs, one, d, e, xact, lda, zero, b, lda )
370 *
371  DO 100 ifact = 1, 2
372  IF( ifact.EQ.1 ) THEN
373  fact = 'F'
374  ELSE
375  fact = 'N'
376  END IF
377 *
378 * Compute the condition number for comparison with
379 * the value returned by SPTSVX.
380 *
381  IF( zerot ) THEN
382  IF( ifact.EQ.1 )
383  $ GO TO 100
384  rcondc = zero
385 *
386  ELSE IF( ifact.EQ.1 ) THEN
387 *
388 * Compute the 1-norm of A.
389 *
390  anorm = slanst( '1', n, d, e )
391 *
392  CALL scopy( n, d, 1, d( n+1 ), 1 )
393  IF( n.GT.1 )
394  $ CALL scopy( n-1, e, 1, e( n+1 ), 1 )
395 *
396 * Factor the matrix A.
397 *
398  CALL spttrf( n, d( n+1 ), e( n+1 ), info )
399 *
400 * Use SPTTRS to solve for one column at a time of
401 * inv(A), computing the maximum column sum as we go.
402 *
403  ainvnm = zero
404  DO 60 i = 1, n
405  DO 50 j = 1, n
406  x( j ) = zero
407  50 CONTINUE
408  x( i ) = one
409  CALL spttrs( n, 1, d( n+1 ), e( n+1 ), x, lda,
410  $ info )
411  ainvnm = max( ainvnm, sasum( n, x, 1 ) )
412  60 CONTINUE
413 *
414 * Compute the 1-norm condition number of A.
415 *
416  IF( anorm.LE.zero .OR. ainvnm.LE.zero ) THEN
417  rcondc = one
418  ELSE
419  rcondc = ( one / anorm ) / ainvnm
420  END IF
421  END IF
422 *
423  IF( ifact.EQ.2 ) THEN
424 *
425 * --- Test SPTSV --
426 *
427  CALL scopy( n, d, 1, d( n+1 ), 1 )
428  IF( n.GT.1 )
429  $ CALL scopy( n-1, e, 1, e( n+1 ), 1 )
430  CALL slacpy( 'Full', n, nrhs, b, lda, x, lda )
431 *
432 * Factor A as L*D*L' and solve the system A*X = B.
433 *
434  srnamt = 'SPTSV '
435  CALL sptsv( n, nrhs, d( n+1 ), e( n+1 ), x, lda,
436  $ info )
437 *
438 * Check error code from SPTSV .
439 *
440  IF( info.NE.izero )
441  $ CALL alaerh( path, 'SPTSV ', info, izero, ' ', n,
442  $ n, 1, 1, nrhs, imat, nfail, nerrs,
443  $ nout )
444  nt = 0
445  IF( izero.EQ.0 ) THEN
446 *
447 * Check the factorization by computing the ratio
448 * norm(L*D*L' - A) / (n * norm(A) * EPS )
449 *
450  CALL sptt01( n, d, e, d( n+1 ), e( n+1 ), work,
451  $ result( 1 ) )
452 *
453 * Compute the residual in the solution.
454 *
455  CALL slacpy( 'Full', n, nrhs, b, lda, work, lda )
456  CALL sptt02( n, nrhs, d, e, x, lda, work, lda,
457  $ result( 2 ) )
458 *
459 * Check solution from generated exact solution.
460 *
461  CALL sget04( n, nrhs, x, lda, xact, lda, rcondc,
462  $ result( 3 ) )
463  nt = 3
464  END IF
465 *
466 * Print information about the tests that did not pass
467 * the threshold.
468 *
469  DO 70 k = 1, nt
470  IF( result( k ).GE.thresh ) THEN
471  IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
472  $ CALL aladhd( nout, path )
473  WRITE( nout, fmt = 9999 )'SPTSV ', n, imat, k,
474  $ result( k )
475  nfail = nfail + 1
476  END IF
477  70 CONTINUE
478  nrun = nrun + nt
479  END IF
480 *
481 * --- Test SPTSVX ---
482 *
483  IF( ifact.GT.1 ) THEN
484 *
485 * Initialize D( N+1:2*N ) and E( N+1:2*N ) to zero.
486 *
487  DO 80 i = 1, n - 1
488  d( n+i ) = zero
489  e( n+i ) = zero
490  80 CONTINUE
491  IF( n.GT.0 )
492  $ d( n+n ) = zero
493  END IF
494 *
495  CALL slaset( 'Full', n, nrhs, zero, zero, x, lda )
496 *
497 * Solve the system and compute the condition number and
498 * error bounds using SPTSVX.
499 *
500  srnamt = 'SPTSVX'
501  CALL sptsvx( fact, n, nrhs, d, e, d( n+1 ), e( n+1 ), b,
502  $ lda, x, lda, rcond, rwork, rwork( nrhs+1 ),
503  $ work, info )
504 *
505 * Check the error code from SPTSVX.
506 *
507  IF( info.NE.izero )
508  $ CALL alaerh( path, 'SPTSVX', info, izero, fact, n, n,
509  $ 1, 1, nrhs, imat, nfail, nerrs, nout )
510  IF( izero.EQ.0 ) THEN
511  IF( ifact.EQ.2 ) THEN
512 *
513 * Check the factorization by computing the ratio
514 * norm(L*D*L' - A) / (n * norm(A) * EPS )
515 *
516  k1 = 1
517  CALL sptt01( n, d, e, d( n+1 ), e( n+1 ), work,
518  $ result( 1 ) )
519  ELSE
520  k1 = 2
521  END IF
522 *
523 * Compute the residual in the solution.
524 *
525  CALL slacpy( 'Full', n, nrhs, b, lda, work, lda )
526  CALL sptt02( n, nrhs, d, e, x, lda, work, lda,
527  $ result( 2 ) )
528 *
529 * Check solution from generated exact solution.
530 *
531  CALL sget04( n, nrhs, x, lda, xact, lda, rcondc,
532  $ result( 3 ) )
533 *
534 * Check error bounds from iterative refinement.
535 *
536  CALL sptt05( n, nrhs, d, e, b, lda, x, lda, xact, lda,
537  $ rwork, rwork( nrhs+1 ), result( 4 ) )
538  ELSE
539  k1 = 6
540  END IF
541 *
542 * Check the reciprocal of the condition number.
543 *
544  result( 6 ) = sget06( rcond, rcondc )
545 *
546 * Print information about the tests that did not pass
547 * the threshold.
548 *
549  DO 90 k = k1, 6
550  IF( result( k ).GE.thresh ) THEN
551  IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
552  $ CALL aladhd( nout, path )
553  WRITE( nout, fmt = 9998 )'SPTSVX', fact, n, imat,
554  $ k, result( k )
555  nfail = nfail + 1
556  END IF
557  90 CONTINUE
558  nrun = nrun + 7 - k1
559  100 CONTINUE
560  110 CONTINUE
561  120 CONTINUE
562 *
563 * Print a summary of the results.
564 *
565  CALL alasvm( path, nout, nfail, nrun, nerrs )
566 *
567  9999 FORMAT( 1x, a, ', N =', i5, ', type ', i2, ', test ', i2,
568  $ ', ratio = ', g12.5 )
569  9998 FORMAT( 1x, a, ', FACT=''', a1, ''', N =', i5, ', type ', i2,
570  $ ', test ', i2, ', ratio = ', g12.5 )
571  RETURN
572 *
573 * End of SDRVPT
574 *
575  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 slaset(UPLO, M, N, ALPHA, BETA, A, LDA)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition: slaset.f:110
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
Definition: slacpy.f:103
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 spttrf(N, D, E, INFO)
SPTTRF
Definition: spttrf.f:91
subroutine slatms(M, N, DIST, ISEED, SYM, D, MODE, COND, DMAX, KL, KU, PACK, A, LDA, WORK, INFO)
SLATMS
Definition: slatms.f:321
subroutine spttrs(N, NRHS, D, E, B, LDB, INFO)
SPTTRS
Definition: spttrs.f:109
subroutine sptsv(N, NRHS, D, E, B, LDB, INFO)
SPTSV computes the solution to system of linear equations A * X = B for PT matrices
Definition: sptsv.f:114
subroutine sptsvx(FACT, N, NRHS, D, E, DF, EF, B, LDB, X, LDX, RCOND, FERR, BERR, WORK, INFO)
SPTSVX computes the solution to system of linear equations A * X = B for PT matrices
Definition: sptsvx.f:228
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:82
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:79
subroutine slatb4(PATH, IMAT, M, N, TYPE, KL, KU, ANORM, MODE, CNDNUM, DIST)
SLATB4
Definition: slatb4.f:120
subroutine serrvx(PATH, NUNIT)
SERRVX
Definition: serrvx.f:55
subroutine sptt05(N, NRHS, D, E, B, LDB, X, LDX, XACT, LDXACT, FERR, BERR, RESLTS)
SPTT05
Definition: sptt05.f:150
subroutine sptt02(N, NRHS, D, E, X, LDX, B, LDB, RESID)
SPTT02
Definition: sptt02.f:104
subroutine sget04(N, NRHS, X, LDX, XACT, LDXACT, RCOND, RESID)
SGET04
Definition: sget04.f:102
subroutine sdrvpt(DOTYPE, NN, NVAL, NRHS, THRESH, TSTERR, A, D, E, B, X, XACT, WORK, RWORK, NOUT)
SDRVPT
Definition: sdrvpt.f:140
subroutine sptt01(N, D, E, DF, EF, WORK, RESID)
SPTT01
Definition: sptt01.f:91
subroutine slaptm(N, NRHS, ALPHA, D, E, X, LDX, BETA, B, LDB)
SLAPTM
Definition: slaptm.f:116