LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
cdrvhe.f
Go to the documentation of this file.
1 *> \brief \b CDRVHE
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 CDRVHE( DOTYPE, NN, NVAL, NRHS, THRESH, TSTERR, NMAX,
12 * A, AFAC, AINV, B, X, XACT, WORK, RWORK, IWORK,
13 * 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 RWORK( * )
24 * COMPLEX A( * ), AFAC( * ), AINV( * ), B( * ),
25 * $ WORK( * ), X( * ), XACT( * )
26 * ..
27 *
28 *
29 *> \par Purpose:
30 * =============
31 *>
32 *> \verbatim
33 *>
34 *> CDRVHE tests the driver routines CHESV 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 COMPLEX array, dimension (NMAX*NMAX)
91 *> \endverbatim
92 *>
93 *> \param[out] AFAC
94 *> \verbatim
95 *> AFAC is COMPLEX array, dimension (NMAX*NMAX)
96 *> \endverbatim
97 *>
98 *> \param[out] AINV
99 *> \verbatim
100 *> AINV is COMPLEX array, dimension (NMAX*NMAX)
101 *> \endverbatim
102 *>
103 *> \param[out] B
104 *> \verbatim
105 *> B is COMPLEX array, dimension (NMAX*NRHS)
106 *> \endverbatim
107 *>
108 *> \param[out] X
109 *> \verbatim
110 *> X is COMPLEX array, dimension (NMAX*NRHS)
111 *> \endverbatim
112 *>
113 *> \param[out] XACT
114 *> \verbatim
115 *> XACT is COMPLEX array, dimension (NMAX*NRHS)
116 *> \endverbatim
117 *>
118 *> \param[out] WORK
119 *> \verbatim
120 *> WORK is COMPLEX array, dimension
121 *> (NMAX*max(2,NRHS))
122 *> \endverbatim
123 *>
124 *> \param[out] RWORK
125 *> \verbatim
126 *> RWORK is REAL array, dimension (NMAX+2*NRHS)
127 *> \endverbatim
128 *>
129 *> \param[out] IWORK
130 *> \verbatim
131 *> IWORK is INTEGER array, dimension (NMAX)
132 *> \endverbatim
133 *>
134 *> \param[in] NOUT
135 *> \verbatim
136 *> NOUT is INTEGER
137 *> The unit number for output.
138 *> \endverbatim
139 *
140 * Authors:
141 * ========
142 *
143 *> \author Univ. of Tennessee
144 *> \author Univ. of California Berkeley
145 *> \author Univ. of Colorado Denver
146 *> \author NAG Ltd.
147 *
148 *> \date November 2011
149 *
150 *> \ingroup complex_lin
151 *
152 * =====================================================================
153  SUBROUTINE cdrvhe( DOTYPE, NN, NVAL, NRHS, THRESH, TSTERR, NMAX,
154  $ a, afac, ainv, b, x, xact, work, rwork, iwork,
155  $ nout )
156 *
157 * -- LAPACK test routine (version 3.4.0) --
158 * -- LAPACK is a software package provided by Univ. of Tennessee, --
159 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
160 * November 2011
161 *
162 * .. Scalar Arguments ..
163  LOGICAL tsterr
164  INTEGER nmax, nn, nout, nrhs
165  REAL thresh
166 * ..
167 * .. Array Arguments ..
168  LOGICAL dotype( * )
169  INTEGER iwork( * ), nval( * )
170  REAL rwork( * )
171  COMPLEX a( * ), afac( * ), ainv( * ), b( * ),
172  $ work( * ), x( * ), xact( * )
173 * ..
174 *
175 * =====================================================================
176 *
177 * .. Parameters ..
178  REAL one, zero
179  parameter( one = 1.0e+0, zero = 0.0e+0 )
180  INTEGER ntypes, ntests
181  parameter( ntypes = 10, ntests = 6 )
182  INTEGER nfact
183  parameter( nfact = 2 )
184 * ..
185 * .. Local Scalars ..
186  LOGICAL zerot
187  CHARACTER dist, fact, type, uplo, xtype
188  CHARACTER*3 path
189  INTEGER i, i1, i2, ifact, imat, in, info, ioff, iuplo,
190  $ izero, j, k, k1, kl, ku, lda, lwork, mode, n,
191  $ nb, nbmin, nerrs, nfail, nimat, nrun, nt
192  REAL ainvnm, anorm, cndnum, rcond, rcondc
193 * ..
194 * .. Local Arrays ..
195  CHARACTER facts( nfact ), uplos( 2 )
196  INTEGER iseed( 4 ), iseedy( 4 )
197  REAL result( ntests )
198 * ..
199 * .. External Functions ..
200  REAL clanhe, sget06
201  EXTERNAL clanhe, sget06
202 * ..
203 * .. External Subroutines ..
204  EXTERNAL aladhd, alaerh, alasvm, cerrvx, cget04, chesv,
205  $ chesvx, chet01, chetrf, chetri2, clacpy,
206  $ claipd, clarhs, claset, clatb4, clatms, cpot02,
207  $ cpot05, xlaenv
208 * ..
209 * .. Scalars in Common ..
210  LOGICAL lerr, ok
211  CHARACTER*32 srnamt
212  INTEGER infot, nunit
213 * ..
214 * .. Common blocks ..
215  common / infoc / infot, nunit, ok, lerr
216  common / srnamc / srnamt
217 * ..
218 * .. Intrinsic Functions ..
219  INTRINSIC cmplx, max, min
220 * ..
221 * .. Data statements ..
222  DATA iseedy / 1988, 1989, 1990, 1991 /
223  DATA uplos / 'U', 'L' / , facts / 'F', 'N' /
224 * ..
225 * .. Executable Statements ..
226 *
227 * Initialize constants and the random number seed.
228 *
229  path( 1: 1 ) = 'C'
230  path( 2: 3 ) = 'HE'
231  nrun = 0
232  nfail = 0
233  nerrs = 0
234  DO 10 i = 1, 4
235  iseed( i ) = iseedy( i )
236  10 continue
237  lwork = max( 2*nmax, nmax*nrhs )
238 *
239 * Test the error exits
240 *
241  IF( tsterr )
242  $ CALL cerrvx( path, nout )
243  infot = 0
244 *
245 * Set the block size and minimum block size for testing.
246 *
247  nb = 1
248  nbmin = 2
249  CALL xlaenv( 1, nb )
250  CALL xlaenv( 2, nbmin )
251 *
252 * Do for each value of N in NVAL
253 *
254  DO 180 in = 1, nn
255  n = nval( in )
256  lda = max( n, 1 )
257  xtype = 'N'
258  nimat = ntypes
259  IF( n.LE.0 )
260  $ nimat = 1
261 *
262  DO 170 imat = 1, nimat
263 *
264 * Do the tests only if DOTYPE( IMAT ) is true.
265 *
266  IF( .NOT.dotype( imat ) )
267  $ go to 170
268 *
269 * Skip types 3, 4, 5, or 6 if the matrix size is too small.
270 *
271  zerot = imat.GE.3 .AND. imat.LE.6
272  IF( zerot .AND. n.LT.imat-2 )
273  $ go to 170
274 *
275 * Do first for UPLO = 'U', then for UPLO = 'L'
276 *
277  DO 160 iuplo = 1, 2
278  uplo = uplos( iuplo )
279 *
280 * Set up parameters with CLATB4 and generate a test matrix
281 * with CLATMS.
282 *
283  CALL clatb4( path, imat, n, n, type, kl, ku, anorm, mode,
284  $ cndnum, dist )
285 *
286  srnamt = 'CLATMS'
287  CALL clatms( n, n, dist, iseed, type, rwork, mode,
288  $ cndnum, anorm, kl, ku, uplo, a, lda, work,
289  $ info )
290 *
291 * Check error code from CLATMS.
292 *
293  IF( info.NE.0 ) THEN
294  CALL alaerh( path, 'CLATMS', info, 0, uplo, n, n, -1,
295  $ -1, -1, imat, nfail, nerrs, nout )
296  go to 160
297  END IF
298 *
299 * For types 3-6, zero one or more rows and columns of the
300 * matrix to test that INFO is returned correctly.
301 *
302  IF( zerot ) THEN
303  IF( imat.EQ.3 ) THEN
304  izero = 1
305  ELSE IF( imat.EQ.4 ) THEN
306  izero = n
307  ELSE
308  izero = n / 2 + 1
309  END IF
310 *
311  IF( imat.LT.6 ) THEN
312 *
313 * Set row and column IZERO to zero.
314 *
315  IF( iuplo.EQ.1 ) THEN
316  ioff = ( izero-1 )*lda
317  DO 20 i = 1, izero - 1
318  a( ioff+i ) = zero
319  20 continue
320  ioff = ioff + izero
321  DO 30 i = izero, n
322  a( ioff ) = zero
323  ioff = ioff + lda
324  30 continue
325  ELSE
326  ioff = izero
327  DO 40 i = 1, izero - 1
328  a( ioff ) = zero
329  ioff = ioff + lda
330  40 continue
331  ioff = ioff - izero
332  DO 50 i = izero, n
333  a( ioff+i ) = zero
334  50 continue
335  END IF
336  ELSE
337  ioff = 0
338  IF( iuplo.EQ.1 ) THEN
339 *
340 * Set the first IZERO rows and columns to zero.
341 *
342  DO 70 j = 1, n
343  i2 = min( j, izero )
344  DO 60 i = 1, i2
345  a( ioff+i ) = zero
346  60 continue
347  ioff = ioff + lda
348  70 continue
349  ELSE
350 *
351 * Set the last IZERO rows and columns to zero.
352 *
353  DO 90 j = 1, n
354  i1 = max( j, izero )
355  DO 80 i = i1, n
356  a( ioff+i ) = zero
357  80 continue
358  ioff = ioff + lda
359  90 continue
360  END IF
361  END IF
362  ELSE
363  izero = 0
364  END IF
365 *
366 * Set the imaginary part of the diagonals.
367 *
368  CALL claipd( n, a, lda+1, 0 )
369 *
370  DO 150 ifact = 1, nfact
371 *
372 * Do first for FACT = 'F', then for other values.
373 *
374  fact = facts( ifact )
375 *
376 * Compute the condition number for comparison with
377 * the value returned by CHESVX.
378 *
379  IF( zerot ) THEN
380  IF( ifact.EQ.1 )
381  $ go to 150
382  rcondc = zero
383 *
384  ELSE IF( ifact.EQ.1 ) THEN
385 *
386 * Compute the 1-norm of A.
387 *
388  anorm = clanhe( '1', uplo, n, a, lda, rwork )
389 *
390 * Factor the matrix A.
391 *
392  CALL clacpy( uplo, n, n, a, lda, afac, lda )
393  CALL chetrf( uplo, n, afac, lda, iwork, work,
394  $ lwork, info )
395 *
396 * Compute inv(A) and take its norm.
397 *
398  CALL clacpy( uplo, n, n, afac, lda, ainv, lda )
399  lwork = (n+nb+1)*(nb+3)
400  CALL chetri2( uplo, n, ainv, lda, iwork, work,
401  $ lwork, info )
402  ainvnm = clanhe( '1', uplo, n, ainv, lda, rwork )
403 *
404 * Compute the 1-norm condition number of A.
405 *
406  IF( anorm.LE.zero .OR. ainvnm.LE.zero ) THEN
407  rcondc = one
408  ELSE
409  rcondc = ( one / anorm ) / ainvnm
410  END IF
411  END IF
412 *
413 * Form an exact solution and set the right hand side.
414 *
415  srnamt = 'CLARHS'
416  CALL clarhs( path, xtype, uplo, ' ', n, n, kl, ku,
417  $ nrhs, a, lda, xact, lda, b, lda, iseed,
418  $ info )
419  xtype = 'C'
420 *
421 * --- Test CHESV ---
422 *
423  IF( ifact.EQ.2 ) THEN
424  CALL clacpy( uplo, n, n, a, lda, afac, lda )
425  CALL clacpy( 'Full', n, nrhs, b, lda, x, lda )
426 *
427 * Factor the matrix and solve the system using CHESV.
428 *
429  srnamt = 'CHESV '
430  CALL chesv( uplo, n, nrhs, afac, lda, iwork, x,
431  $ lda, work, lwork, info )
432 *
433 * Adjust the expected value of INFO to account for
434 * pivoting.
435 *
436  k = izero
437  IF( k.GT.0 ) THEN
438  100 continue
439  IF( iwork( k ).LT.0 ) THEN
440  IF( iwork( k ).NE.-k ) THEN
441  k = -iwork( k )
442  go to 100
443  END IF
444  ELSE IF( iwork( k ).NE.k ) THEN
445  k = iwork( k )
446  go to 100
447  END IF
448  END IF
449 *
450 * Check error code from CHESV .
451 *
452  IF( info.NE.k ) THEN
453  CALL alaerh( path, 'CHESV ', info, k, uplo, n,
454  $ n, -1, -1, nrhs, imat, nfail,
455  $ nerrs, nout )
456  go to 120
457  ELSE IF( info.NE.0 ) THEN
458  go to 120
459  END IF
460 *
461 * Reconstruct matrix from factors and compute
462 * residual.
463 *
464  CALL chet01( uplo, n, a, lda, afac, lda, iwork,
465  $ ainv, lda, rwork, result( 1 ) )
466 *
467 * Compute residual of the computed solution.
468 *
469  CALL clacpy( 'Full', n, nrhs, b, lda, work, lda )
470  CALL cpot02( uplo, n, nrhs, a, lda, x, lda, work,
471  $ lda, rwork, result( 2 ) )
472 *
473 * Check solution from generated exact solution.
474 *
475  CALL cget04( n, nrhs, x, lda, xact, lda, rcondc,
476  $ result( 3 ) )
477  nt = 3
478 *
479 * Print information about the tests that did not pass
480 * the threshold.
481 *
482  DO 110 k = 1, nt
483  IF( result( k ).GE.thresh ) THEN
484  IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
485  $ CALL aladhd( nout, path )
486  WRITE( nout, fmt = 9999 )'CHESV ', uplo, n,
487  $ imat, k, result( k )
488  nfail = nfail + 1
489  END IF
490  110 continue
491  nrun = nrun + nt
492  120 continue
493  END IF
494 *
495 * --- Test CHESVX ---
496 *
497  IF( ifact.EQ.2 )
498  $ CALL claset( uplo, n, n, cmplx( zero ),
499  $ cmplx( zero ), afac, lda )
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 CHESVX.
505 *
506  srnamt = 'CHESVX'
507  CALL chesvx( fact, uplo, n, nrhs, a, lda, afac, lda,
508  $ iwork, b, lda, x, lda, rcond, rwork,
509  $ rwork( nrhs+1 ), work, lwork,
510  $ rwork( 2*nrhs+1 ), info )
511 *
512 * Adjust the expected value of INFO to account for
513 * pivoting.
514 *
515  k = izero
516  IF( k.GT.0 ) THEN
517  130 continue
518  IF( iwork( k ).LT.0 ) THEN
519  IF( iwork( k ).NE.-k ) THEN
520  k = -iwork( k )
521  go to 130
522  END IF
523  ELSE IF( iwork( k ).NE.k ) THEN
524  k = iwork( k )
525  go to 130
526  END IF
527  END IF
528 *
529 * Check the error code from CHESVX.
530 *
531  IF( info.NE.k ) THEN
532  CALL alaerh( path, 'CHESVX', info, k, fact // uplo,
533  $ n, n, -1, -1, nrhs, imat, nfail,
534  $ nerrs, nout )
535  go to 150
536  END IF
537 *
538  IF( info.EQ.0 ) THEN
539  IF( ifact.GE.2 ) THEN
540 *
541 * Reconstruct matrix from factors and compute
542 * residual.
543 *
544  CALL chet01( uplo, n, a, lda, afac, lda, iwork,
545  $ ainv, lda, rwork( 2*nrhs+1 ),
546  $ result( 1 ) )
547  k1 = 1
548  ELSE
549  k1 = 2
550  END IF
551 *
552 * Compute residual of the computed solution.
553 *
554  CALL clacpy( 'Full', n, nrhs, b, lda, work, lda )
555  CALL cpot02( uplo, n, nrhs, a, lda, x, lda, work,
556  $ lda, rwork( 2*nrhs+1 ), result( 2 ) )
557 *
558 * Check solution from generated exact solution.
559 *
560  CALL cget04( n, nrhs, x, lda, xact, lda, rcondc,
561  $ result( 3 ) )
562 *
563 * Check the error bounds from iterative refinement.
564 *
565  CALL cpot05( uplo, n, nrhs, a, lda, b, lda, x, lda,
566  $ xact, lda, rwork, rwork( nrhs+1 ),
567  $ result( 4 ) )
568  ELSE
569  k1 = 6
570  END IF
571 *
572 * Compare RCOND from CHESVX 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 140 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  WRITE( nout, fmt = 9998 )'CHESVX', fact, uplo,
585  $ n, imat, k, result( k )
586  nfail = nfail + 1
587  END IF
588  140 continue
589  nrun = nrun + 7 - k1
590 *
591  150 continue
592 *
593  160 continue
594  170 continue
595  180 continue
596 *
597 * Print a summary of the results.
598 *
599  CALL alasvm( path, nout, nfail, nrun, nerrs )
600 *
601  9999 format( 1x, a, ', UPLO=''', a1, ''', N =', i5, ', type ', i2,
602  $ ', test ', i2, ', ratio =', g12.5 )
603  9998 format( 1x, a, ', FACT=''', a1, ''', UPLO=''', a1, ''', N =', i5,
604  $ ', type ', i2, ', test ', i2, ', ratio =', g12.5 )
605  return
606 *
607 * End of CDRVHE
608 *
609  END