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