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