LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
sdrvpp.f
Go to the documentation of this file.
1 *> \brief \b SDRVPP
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 SDRVPP( DOTYPE, NN, NVAL, NRHS, THRESH, TSTERR, NMAX,
12 * A, AFAC, ASAV, B, BSAV, X, XACT, S, WORK,
13 * RWORK, IWORK, NOUT )
14 *
15 * .. Scalar Arguments ..
16 * LOGICAL TSTERR
17 * INTEGER NMAX, NN, NOUT, NRHS
18 * REAL THRESH
19 * ..
20 * .. Array Arguments ..
21 * LOGICAL DOTYPE( * )
22 * INTEGER IWORK( * ), NVAL( * )
23 * REAL A( * ), AFAC( * ), ASAV( * ), B( * ),
24 * $ BSAV( * ), RWORK( * ), S( * ), WORK( * ),
25 * $ X( * ), XACT( * )
26 * ..
27 *
28 *
29 *> \par Purpose:
30 * =============
31 *>
32 *> \verbatim
33 *>
34 *> SDRVPP tests the driver routines SPPSV and -SVX.
35 *> \endverbatim
36 *
37 * Arguments:
38 * ==========
39 *
40 *> \param[in] DOTYPE
41 *> \verbatim
42 *> DOTYPE is LOGICAL array, dimension (NTYPES)
43 *> The matrix types to be used for testing. Matrices of type j
44 *> (for 1 <= j <= NTYPES) are used for testing if DOTYPE(j) =
45 *> .TRUE.; if DOTYPE(j) = .FALSE., then type j is not used.
46 *> \endverbatim
47 *>
48 *> \param[in] NN
49 *> \verbatim
50 *> NN is INTEGER
51 *> The number of values of N contained in the vector NVAL.
52 *> \endverbatim
53 *>
54 *> \param[in] NVAL
55 *> \verbatim
56 *> NVAL is INTEGER array, dimension (NN)
57 *> The values of the matrix dimension N.
58 *> \endverbatim
59 *>
60 *> \param[in] NRHS
61 *> \verbatim
62 *> NRHS is INTEGER
63 *> The number of right hand side vectors to be generated for
64 *> each linear system.
65 *> \endverbatim
66 *>
67 *> \param[in] THRESH
68 *> \verbatim
69 *> THRESH is REAL
70 *> The threshold value for the test ratios. A result is
71 *> included in the output file if RESULT >= THRESH. To have
72 *> every test ratio printed, use THRESH = 0.
73 *> \endverbatim
74 *>
75 *> \param[in] TSTERR
76 *> \verbatim
77 *> TSTERR is LOGICAL
78 *> Flag that indicates whether error exits are to be tested.
79 *> \endverbatim
80 *>
81 *> \param[in] NMAX
82 *> \verbatim
83 *> NMAX is INTEGER
84 *> The maximum value permitted for N, used in dimensioning the
85 *> work arrays.
86 *> \endverbatim
87 *>
88 *> \param[out] A
89 *> \verbatim
90 *> A is REAL array, dimension
91 *> (NMAX*(NMAX+1)/2)
92 *> \endverbatim
93 *>
94 *> \param[out] AFAC
95 *> \verbatim
96 *> AFAC is REAL array, dimension
97 *> (NMAX*(NMAX+1)/2)
98 *> \endverbatim
99 *>
100 *> \param[out] ASAV
101 *> \verbatim
102 *> ASAV is REAL array, dimension
103 *> (NMAX*(NMAX+1)/2)
104 *> \endverbatim
105 *>
106 *> \param[out] B
107 *> \verbatim
108 *> B is REAL array, dimension (NMAX*NRHS)
109 *> \endverbatim
110 *>
111 *> \param[out] BSAV
112 *> \verbatim
113 *> BSAV is REAL array, dimension (NMAX*NRHS)
114 *> \endverbatim
115 *>
116 *> \param[out] X
117 *> \verbatim
118 *> X is REAL array, dimension (NMAX*NRHS)
119 *> \endverbatim
120 *>
121 *> \param[out] XACT
122 *> \verbatim
123 *> XACT is REAL array, dimension (NMAX*NRHS)
124 *> \endverbatim
125 *>
126 *> \param[out] S
127 *> \verbatim
128 *> S is REAL array, dimension (NMAX)
129 *> \endverbatim
130 *>
131 *> \param[out] WORK
132 *> \verbatim
133 *> WORK is REAL array, dimension
134 *> (NMAX*max(3,NRHS))
135 *> \endverbatim
136 *>
137 *> \param[out] RWORK
138 *> \verbatim
139 *> RWORK is REAL array, dimension (NMAX+2*NRHS)
140 *> \endverbatim
141 *>
142 *> \param[out] IWORK
143 *> \verbatim
144 *> IWORK is INTEGER array, dimension (NMAX)
145 *> \endverbatim
146 *>
147 *> \param[in] NOUT
148 *> \verbatim
149 *> NOUT is INTEGER
150 *> The unit number for output.
151 *> \endverbatim
152 *
153 * Authors:
154 * ========
155 *
156 *> \author Univ. of Tennessee
157 *> \author Univ. of California Berkeley
158 *> \author Univ. of Colorado Denver
159 *> \author NAG Ltd.
160 *
161 *> \date November 2011
162 *
163 *> \ingroup single_lin
164 *
165 * =====================================================================
166  SUBROUTINE sdrvpp( DOTYPE, NN, NVAL, NRHS, THRESH, TSTERR, NMAX,
167  $ a, afac, asav, b, bsav, x, xact, s, work,
168  $ rwork, iwork, nout )
169 *
170 * -- LAPACK test routine (version 3.4.0) --
171 * -- LAPACK is a software package provided by Univ. of Tennessee, --
172 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
173 * November 2011
174 *
175 * .. Scalar Arguments ..
176  LOGICAL tsterr
177  INTEGER nmax, nn, nout, nrhs
178  REAL thresh
179 * ..
180 * .. Array Arguments ..
181  LOGICAL dotype( * )
182  INTEGER iwork( * ), nval( * )
183  REAL a( * ), afac( * ), asav( * ), b( * ),
184  $ bsav( * ), rwork( * ), s( * ), work( * ),
185  $ x( * ), xact( * )
186 * ..
187 *
188 * =====================================================================
189 *
190 * .. Parameters ..
191  REAL one, zero
192  parameter( one = 1.0e+0, zero = 0.0e+0 )
193  INTEGER ntypes
194  parameter( ntypes = 9 )
195  INTEGER ntests
196  parameter( ntests = 6 )
197 * ..
198 * .. Local Scalars ..
199  LOGICAL equil, nofact, prefac, zerot
200  CHARACTER dist, equed, fact, packit, type, uplo, xtype
201  CHARACTER*3 path
202  INTEGER i, iequed, ifact, imat, in, info, ioff, iuplo,
203  $ izero, k, k1, kl, ku, lda, mode, n, nerrs,
204  $ nfact, nfail, nimat, npp, nrun, nt
205  REAL ainvnm, amax, anorm, cndnum, rcond, rcondc,
206  $ roldc, scond
207 * ..
208 * .. Local Arrays ..
209  CHARACTER equeds( 2 ), facts( 3 ), packs( 2 ), uplos( 2 )
210  INTEGER iseed( 4 ), iseedy( 4 )
211  REAL result( ntests )
212 * ..
213 * .. External Functions ..
214  LOGICAL lsame
215  REAL sget06, slansp
216  EXTERNAL lsame, sget06, slansp
217 * ..
218 * .. External Subroutines ..
219  EXTERNAL aladhd, alaerh, alasvm, scopy, serrvx, sget04,
220  $ slacpy, slaqsp, slarhs, slaset, slatb4, slatms,
221  $ sppequ, sppsv, sppsvx, sppt01, sppt02, sppt05,
222  $ spptrf, spptri
223 * ..
224 * .. Scalars in Common ..
225  LOGICAL lerr, ok
226  CHARACTER*32 srnamt
227  INTEGER infot, nunit
228 * ..
229 * .. Common blocks ..
230  common / infoc / infot, nunit, ok, lerr
231  common / srnamc / srnamt
232 * ..
233 * .. Intrinsic Functions ..
234  INTRINSIC max
235 * ..
236 * .. Data statements ..
237  DATA iseedy / 1988, 1989, 1990, 1991 /
238  DATA uplos / 'U', 'L' / , facts / 'F', 'N', 'E' / ,
239  $ packs / 'C', 'R' / , equeds / 'N', 'Y' /
240 * ..
241 * .. Executable Statements ..
242 *
243 * Initialize constants and the random number seed.
244 *
245  path( 1: 1 ) = 'Single precision'
246  path( 2: 3 ) = 'PP'
247  nrun = 0
248  nfail = 0
249  nerrs = 0
250  DO 10 i = 1, 4
251  iseed( i ) = iseedy( i )
252  10 continue
253 *
254 * Test the error exits
255 *
256  IF( tsterr )
257  $ CALL serrvx( path, nout )
258  infot = 0
259 *
260 * Do for each value of N in NVAL
261 *
262  DO 140 in = 1, nn
263  n = nval( in )
264  lda = max( n, 1 )
265  npp = n*( n+1 ) / 2
266  xtype = 'N'
267  nimat = ntypes
268  IF( n.LE.0 )
269  $ nimat = 1
270 *
271  DO 130 imat = 1, nimat
272 *
273 * Do the tests only if DOTYPE( IMAT ) is true.
274 *
275  IF( .NOT.dotype( imat ) )
276  $ go to 130
277 *
278 * Skip types 3, 4, or 5 if the matrix size is too small.
279 *
280  zerot = imat.GE.3 .AND. imat.LE.5
281  IF( zerot .AND. n.LT.imat-2 )
282  $ go to 130
283 *
284 * Do first for UPLO = 'U', then for UPLO = 'L'
285 *
286  DO 120 iuplo = 1, 2
287  uplo = uplos( iuplo )
288  packit = packs( iuplo )
289 *
290 * Set up parameters with SLATB4 and generate a test matrix
291 * with SLATMS.
292 *
293  CALL slatb4( path, imat, n, n, type, kl, ku, anorm, mode,
294  $ cndnum, dist )
295  rcondc = one / cndnum
296 *
297  srnamt = 'SLATMS'
298  CALL slatms( n, n, dist, iseed, type, rwork, mode,
299  $ cndnum, anorm, kl, ku, packit, a, lda, work,
300  $ info )
301 *
302 * Check error code from SLATMS.
303 *
304  IF( info.NE.0 ) THEN
305  CALL alaerh( path, 'SLATMS', info, 0, uplo, n, n, -1,
306  $ -1, -1, imat, nfail, nerrs, nout )
307  go to 120
308  END IF
309 *
310 * For types 3-5, zero one row and column of the matrix to
311 * test that INFO is returned correctly.
312 *
313  IF( zerot ) THEN
314  IF( imat.EQ.3 ) THEN
315  izero = 1
316  ELSE IF( imat.EQ.4 ) THEN
317  izero = n
318  ELSE
319  izero = n / 2 + 1
320  END IF
321 *
322 * Set row and column IZERO of A to 0.
323 *
324  IF( iuplo.EQ.1 ) THEN
325  ioff = ( izero-1 )*izero / 2
326  DO 20 i = 1, izero - 1
327  a( ioff+i ) = zero
328  20 continue
329  ioff = ioff + izero
330  DO 30 i = izero, n
331  a( ioff ) = zero
332  ioff = ioff + i
333  30 continue
334  ELSE
335  ioff = izero
336  DO 40 i = 1, izero - 1
337  a( ioff ) = zero
338  ioff = ioff + n - i
339  40 continue
340  ioff = ioff - izero
341  DO 50 i = izero, n
342  a( ioff+i ) = zero
343  50 continue
344  END IF
345  ELSE
346  izero = 0
347  END IF
348 *
349 * Save a copy of the matrix A in ASAV.
350 *
351  CALL scopy( npp, a, 1, asav, 1 )
352 *
353  DO 110 iequed = 1, 2
354  equed = equeds( iequed )
355  IF( iequed.EQ.1 ) THEN
356  nfact = 3
357  ELSE
358  nfact = 1
359  END IF
360 *
361  DO 100 ifact = 1, nfact
362  fact = facts( ifact )
363  prefac = lsame( fact, 'F' )
364  nofact = lsame( fact, 'N' )
365  equil = lsame( fact, 'E' )
366 *
367  IF( zerot ) THEN
368  IF( prefac )
369  $ go to 100
370  rcondc = zero
371 *
372  ELSE IF( .NOT.lsame( fact, 'N' ) ) THEN
373 *
374 * Compute the condition number for comparison with
375 * the value returned by SPPSVX (FACT = 'N' reuses
376 * the condition number from the previous iteration
377 * with FACT = 'F').
378 *
379  CALL scopy( npp, asav, 1, afac, 1 )
380  IF( equil .OR. iequed.GT.1 ) THEN
381 *
382 * Compute row and column scale factors to
383 * equilibrate the matrix A.
384 *
385  CALL sppequ( uplo, n, afac, s, scond, amax,
386  $ info )
387  IF( info.EQ.0 .AND. n.GT.0 ) THEN
388  IF( iequed.GT.1 )
389  $ scond = zero
390 *
391 * Equilibrate the matrix.
392 *
393  CALL slaqsp( uplo, n, afac, s, scond,
394  $ amax, equed )
395  END IF
396  END IF
397 *
398 * Save the condition number of the
399 * non-equilibrated system for use in SGET04.
400 *
401  IF( equil )
402  $ roldc = rcondc
403 *
404 * Compute the 1-norm of A.
405 *
406  anorm = slansp( '1', uplo, n, afac, rwork )
407 *
408 * Factor the matrix A.
409 *
410  CALL spptrf( uplo, n, afac, info )
411 *
412 * Form the inverse of A.
413 *
414  CALL scopy( npp, afac, 1, a, 1 )
415  CALL spptri( uplo, n, a, info )
416 *
417 * Compute the 1-norm condition number of A.
418 *
419  ainvnm = slansp( '1', uplo, n, a, rwork )
420  IF( anorm.LE.zero .OR. ainvnm.LE.zero ) THEN
421  rcondc = one
422  ELSE
423  rcondc = ( one / anorm ) / ainvnm
424  END IF
425  END IF
426 *
427 * Restore the matrix A.
428 *
429  CALL scopy( npp, asav, 1, a, 1 )
430 *
431 * Form an exact solution and set the right hand side.
432 *
433  srnamt = 'SLARHS'
434  CALL slarhs( path, xtype, uplo, ' ', n, n, kl, ku,
435  $ nrhs, a, lda, xact, lda, b, lda,
436  $ iseed, info )
437  xtype = 'C'
438  CALL slacpy( 'Full', n, nrhs, b, lda, bsav, lda )
439 *
440  IF( nofact ) THEN
441 *
442 * --- Test SPPSV ---
443 *
444 * Compute the L*L' or U'*U factorization of the
445 * matrix and solve the system.
446 *
447  CALL scopy( npp, a, 1, afac, 1 )
448  CALL slacpy( 'Full', n, nrhs, b, lda, x, lda )
449 *
450  srnamt = 'SPPSV '
451  CALL sppsv( uplo, n, nrhs, afac, x, lda, info )
452 *
453 * Check error code from SPPSV .
454 *
455  IF( info.NE.izero ) THEN
456  CALL alaerh( path, 'SPPSV ', info, izero,
457  $ uplo, n, n, -1, -1, nrhs, imat,
458  $ nfail, nerrs, nout )
459  go to 70
460  ELSE IF( info.NE.0 ) THEN
461  go to 70
462  END IF
463 *
464 * Reconstruct matrix from factors and compute
465 * residual.
466 *
467  CALL sppt01( uplo, n, a, afac, rwork,
468  $ result( 1 ) )
469 *
470 * Compute residual of the computed solution.
471 *
472  CALL slacpy( 'Full', n, nrhs, b, lda, work,
473  $ lda )
474  CALL sppt02( uplo, n, nrhs, a, x, lda, work,
475  $ lda, rwork, result( 2 ) )
476 *
477 * Check solution from generated exact solution.
478 *
479  CALL sget04( n, nrhs, x, lda, xact, lda, rcondc,
480  $ result( 3 ) )
481  nt = 3
482 *
483 * Print information about the tests that did not
484 * pass the threshold.
485 *
486  DO 60 k = 1, nt
487  IF( result( k ).GE.thresh ) THEN
488  IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
489  $ CALL aladhd( nout, path )
490  WRITE( nout, fmt = 9999 )'SPPSV ', uplo,
491  $ n, imat, k, result( k )
492  nfail = nfail + 1
493  END IF
494  60 continue
495  nrun = nrun + nt
496  70 continue
497  END IF
498 *
499 * --- Test SPPSVX ---
500 *
501  IF( .NOT.prefac .AND. npp.GT.0 )
502  $ CALL slaset( 'Full', npp, 1, zero, zero, afac,
503  $ npp )
504  CALL slaset( 'Full', n, nrhs, zero, zero, x, lda )
505  IF( iequed.GT.1 .AND. n.GT.0 ) THEN
506 *
507 * Equilibrate the matrix if FACT='F' and
508 * EQUED='Y'.
509 *
510  CALL slaqsp( uplo, n, a, s, scond, amax, equed )
511  END IF
512 *
513 * Solve the system and compute the condition number
514 * and error bounds using SPPSVX.
515 *
516  srnamt = 'SPPSVX'
517  CALL sppsvx( fact, uplo, n, nrhs, a, afac, equed,
518  $ s, b, lda, x, lda, rcond, rwork,
519  $ rwork( nrhs+1 ), work, iwork, info )
520 *
521 * Check the error code from SPPSVX.
522 *
523  IF( info.NE.izero ) THEN
524  CALL alaerh( path, 'SPPSVX', info, izero,
525  $ fact // uplo, n, n, -1, -1, nrhs,
526  $ imat, nfail, nerrs, nout )
527  go to 90
528  END IF
529 *
530  IF( info.EQ.0 ) THEN
531  IF( .NOT.prefac ) THEN
532 *
533 * Reconstruct matrix from factors and compute
534 * residual.
535 *
536  CALL sppt01( uplo, n, a, afac,
537  $ rwork( 2*nrhs+1 ), result( 1 ) )
538  k1 = 1
539  ELSE
540  k1 = 2
541  END IF
542 *
543 * Compute residual of the computed solution.
544 *
545  CALL slacpy( 'Full', n, nrhs, bsav, lda, work,
546  $ lda )
547  CALL sppt02( uplo, n, nrhs, asav, x, lda, work,
548  $ lda, rwork( 2*nrhs+1 ),
549  $ result( 2 ) )
550 *
551 * Check solution from generated exact solution.
552 *
553  IF( nofact .OR. ( prefac .AND. lsame( equed,
554  $ 'N' ) ) ) THEN
555  CALL sget04( n, nrhs, x, lda, xact, lda,
556  $ rcondc, result( 3 ) )
557  ELSE
558  CALL sget04( n, nrhs, x, lda, xact, lda,
559  $ roldc, result( 3 ) )
560  END IF
561 *
562 * Check the error bounds from iterative
563 * refinement.
564 *
565  CALL sppt05( uplo, n, nrhs, asav, b, lda, x,
566  $ lda, xact, lda, rwork,
567  $ rwork( nrhs+1 ), result( 4 ) )
568  ELSE
569  k1 = 6
570  END IF
571 *
572 * Compare RCOND from SPPSVX with the computed value
573 * in RCONDC.
574 *
575  result( 6 ) = sget06( rcond, rcondc )
576 *
577 * Print information about the tests that did not pass
578 * the threshold.
579 *
580  DO 80 k = k1, 6
581  IF( result( k ).GE.thresh ) THEN
582  IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
583  $ CALL aladhd( nout, path )
584  IF( prefac ) THEN
585  WRITE( nout, fmt = 9997 )'SPPSVX', fact,
586  $ uplo, n, equed, imat, k, result( k )
587  ELSE
588  WRITE( nout, fmt = 9998 )'SPPSVX', fact,
589  $ uplo, n, imat, k, result( k )
590  END IF
591  nfail = nfail + 1
592  END IF
593  80 continue
594  nrun = nrun + 7 - k1
595  90 continue
596  100 continue
597  110 continue
598  120 continue
599  130 continue
600  140 continue
601 *
602 * Print a summary of the results.
603 *
604  CALL alasvm( path, nout, nfail, nrun, nerrs )
605 *
606  9999 format( 1x, a, ', UPLO=''', a1, ''', N =', i5, ', type ', i1,
607  $ ', test(', i1, ')=', g12.5 )
608  9998 format( 1x, a, ', FACT=''', a1, ''', UPLO=''', a1, ''', N=', i5,
609  $ ', type ', i1, ', test(', i1, ')=', g12.5 )
610  9997 format( 1x, a, ', FACT=''', a1, ''', UPLO=''', a1, ''', N=', i5,
611  $ ', EQUED=''', a1, ''', type ', i1, ', test(', i1, ')=',
612  $ g12.5 )
613  return
614 *
615 * End of SDRVPP
616 *
617  END