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
dchkpo.f
Go to the documentation of this file.
1  SUBROUTINE dchkpo( DOTYPE, NN, NVAL, NNB, NBVAL, NNS, NSVAL,
2  $ thresh, tsterr, nmax, a, afac, ainv, b, x,
3  $ xact, work, rwork, iwork, nout )
4 *
5  include 'plasmaf.h'
6 *
7 *
8 * -- LAPACK test routine (version 3.1) --
9 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
10 * November 2006
11 *
12 * .. Scalar Arguments ..
13  LOGICAL tsterr
14  INTEGER nmax, nn, nnb, nns, nout
15  DOUBLE PRECISION thresh
16 * ..
17 * .. Array Arguments ..
18  LOGICAL dotype( * )
19  INTEGER iwork( * ), nbval( * ), nsval( * ), nval( * )
20  DOUBLE PRECISION a( * ), afac( * ), ainv( * ), b( * ),
21  $ rwork( * ), work( * ), x( * ), xact( * )
22 * ..
23 *
24 * Purpose
25 * =======
26 *
27 * DCHKPO tests DPOTRF, -TRI, -TRS, -RFS, and -CON
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 * NNB (input) INTEGER
44 * The number of values of NB contained in the vector NBVAL.
45 *
46 * NBVAL (input) INTEGER array, dimension (NBVAL)
47 * The values of the blocksize NB.
48 *
49 * NNS (input) INTEGER
50 * The number of values of NRHS contained in the vector NSVAL.
51 *
52 * NSVAL (input) INTEGER array, dimension (NNS)
53 * The values of the number of right hand sides NRHS.
54 *
55 * THRESH (input) DOUBLE PRECISION
56 * The threshold value for the test ratios. A result is
57 * included in the output file if RESULT >= THRESH. To have
58 * every test ratio printed, use THRESH = 0.
59 *
60 * TSTERR (input) LOGICAL
61 * Flag that indicates whether error exits are to be tested.
62 *
63 * NMAX (input) INTEGER
64 * The maximum value permitted for N, used in dimensioning the
65 * work arrays.
66 *
67 * A (workspace) DOUBLE PRECISION array, dimension (NMAX*NMAX)
68 *
69 * AFAC (workspace) DOUBLE PRECISION array, dimension (NMAX*NMAX)
70 *
71 * AINV (workspace) DOUBLE PRECISION array, dimension (NMAX*NMAX)
72 *
73 * B (workspace) DOUBLE PRECISION array, dimension (NMAX*NSMAX)
74 * where NSMAX is the largest entry in NSVAL.
75 *
76 * X (workspace) DOUBLE PRECISION array, dimension (NMAX*NSMAX)
77 *
78 * XACT (workspace) DOUBLE PRECISION array, dimension (NMAX*NSMAX)
79 *
80 * WORK (workspace) DOUBLE PRECISION array, dimension
81 * (NMAX*max(3,NSMAX))
82 *
83 * RWORK (workspace) DOUBLE PRECISION array, dimension
84 * (max(NMAX,2*NSMAX))
85 *
86 * IWORK (workspace) INTEGER array, dimension (NMAX)
87 *
88 * NOUT (input) INTEGER
89 * The unit number for output.
90 *
91 * =====================================================================
92 *
93 * .. Parameters ..
94  DOUBLE PRECISION zero
95  parameter( zero = 0.0d+0 )
96  INTEGER ntypes
97  parameter( ntypes = 9 )
98  INTEGER ntests
99  parameter( ntests = 8 )
100 * ..
101 * .. Local Scalars ..
102  LOGICAL zerot
103  CHARACTER dist, type, uplo, xtype
104  CHARACTER*3 path
105  INTEGER i, imat, in, inb, info, ioff, irhs, iuplo,
106  $ izero, k, kl, ku, lda, mode, n, nb, nerrs,
107  $ nfail, nimat, nrhs, nrun
108  INTEGER plasma_uplo
109  DOUBLE PRECISION anorm, cndnum, rcond, rcondc
110 * ..
111 * .. Local Arrays ..
112  CHARACTER uplos( 2 )
113  INTEGER plasma_uplos( 2 )
114  INTEGER iseed( 4 ), iseedy( 4 )
115  DOUBLE PRECISION result( ntests )
116 * ..
117 * .. External Functions ..
118  DOUBLE PRECISION dget06, dlansy
119  EXTERNAL dget06, dlansy
120 * ..
121 * .. External Subroutines ..
122  EXTERNAL alaerh, alahd, alasum, derrpo, dget04, dlacpy,
123  $ dlarhs, dlatb4, dlatms, dpocon, dporfs, dpot01,
124  $ dpot02, dpot03, dpot05,dpotrf, dpotri, dpotrs,
125  $ xlaenv
126 * ..
127 * .. Scalars in Common ..
128  LOGICAL lerr, ok
129  CHARACTER*32 srnamt
130  INTEGER infot, nunit
131 * ..
132 * .. Common blocks ..
133  common / infoc / infot, nunit, ok, lerr
134  common / srnamc / srnamt
135 * ..
136 * .. Intrinsic Functions ..
137  INTRINSIC max
138 * ..
139 * .. Data statements ..
140  DATA iseedy / 1988, 1989, 1990, 1991 /
141  DATA uplos / 'U', 'L' /
142  DATA plasma_uplos / plasmaupper, plasmalower /
143 * ..
144 * .. Executable Statements ..
145 *
146 * Initialize constants and the random number seed.
147 *
148  path( 1: 1 ) = 'Double precision'
149  path( 2: 3 ) = 'PO'
150  nrun = 0
151  nfail = 0
152  nerrs = 0
153  DO 10 i = 1, 4
154  iseed( i ) = iseedy( i )
155  10 continue
156 *
157 * Test the error exits
158 *
159  IF( tsterr )
160  $ CALL derrpo( path, nout )
161  infot = 0
162  CALL xlaenv( 2, 2 )
163 *
164 * Do for each value of N in NVAL
165 *
166  DO 120 in = 1, nn
167  n = nval( in )
168  lda = max( n, 1 )
169  xtype = 'N'
170  nimat = ntypes
171  IF( n.LE.0 )
172  $ nimat = 1
173 *
174  izero = 0
175  DO 110 imat = 1, nimat
176 *
177 * Do the tests only if DOTYPE( IMAT ) is true.
178 *
179  IF( .NOT.dotype( imat ) )
180  $ go to 110
181 *
182 * Skip types 3, 4, or 5 if the matrix size is too small.
183 *
184  zerot = imat.GE.3 .AND. imat.LE.5
185  IF( zerot .AND. n.LT.imat-2 )
186  $ go to 110
187 *
188 * Do first for UPLO = 'U', then for UPLO = 'L'
189 *
190  DO 100 iuplo = 1, 2
191  uplo = uplos( iuplo )
192  plasma_uplo = plasma_uplos( iuplo )
193 *
194 * Set up parameters with DLATB4 and generate a test matrix
195 * with DLATMS.
196 *
197  CALL dlatb4( path, imat, n, n, type, kl, ku, anorm, mode,
198  $ cndnum, dist )
199 *
200  srnamt = 'DLATMS'
201  CALL dlatms( n, n, dist, iseed, type, rwork, mode,
202  $ cndnum, anorm, kl, ku, uplo, a, lda, work,
203  $ info )
204 *
205 * Check error code from DLATMS.
206 *
207  IF( info.NE.0 ) THEN
208  CALL alaerh( path, 'DLATMS', info, 0, uplo, n, n, -1,
209  $ -1, -1, imat, nfail, nerrs, nout )
210  go to 100
211  END IF
212 *
213 * For types 3-5, zero one row and column of the matrix to
214 * test that INFO is returned correctly.
215 *
216  IF( zerot ) THEN
217  IF( imat.EQ.3 ) THEN
218  izero = 1
219  ELSE IF( imat.EQ.4 ) THEN
220  izero = n
221  ELSE
222  izero = n / 2 + 1
223  END IF
224  ioff = ( izero-1 )*lda
225 *
226 * Set row and column IZERO of A to 0.
227 *
228  IF( iuplo.EQ.1 ) THEN
229  DO 20 i = 1, izero - 1
230  a( ioff+i ) = zero
231  20 continue
232  ioff = ioff + izero
233  DO 30 i = izero, n
234  a( ioff ) = zero
235  ioff = ioff + lda
236  30 continue
237  ELSE
238  ioff = izero
239  DO 40 i = 1, izero - 1
240  a( ioff ) = zero
241  ioff = ioff + lda
242  40 continue
243  ioff = ioff - izero
244  DO 50 i = izero, n
245  a( ioff+i ) = zero
246  50 continue
247  END IF
248  ELSE
249  izero = 0
250  END IF
251 *
252 * Do for each value of NB in NBVAL
253 *
254  DO 90 inb = 1, nnb
255  nb = nbval( inb )
256  CALL xlaenv( 1, nb )
257  IF ( (n / 25) .GT. nb ) THEN
258  goto 90
259  END IF
260  CALL plasma_set( plasma_tile_size, nb, info )
261 *
262 * Compute the L*L' or U'*U factorization of the matrix.
263 *
264  CALL dlacpy( uplo, n, n, a, lda, afac, lda )
265  srnamt = 'DPOTRF'
266  CALL plasma_dpotrf( plasma_uplo, n, afac, lda, info )
267 *
268 * Check error code from DPOTRF.
269 *
270  IF( info.NE.izero ) THEN
271  CALL alaerh( path, 'DPOTRF', info, izero, uplo, n,
272  $ n, -1, -1, nb, imat, nfail, nerrs,
273  $ nout )
274  go to 90
275  END IF
276 *
277 * Skip the tests if INFO is not 0.
278 *
279  IF( info.NE.0 )
280  $ go to 90
281 *
282 *+ TEST 1
283 * Reconstruct matrix from factors and compute residual.
284 *
285  CALL dlacpy( uplo, n, n, afac, lda, ainv, lda )
286  CALL dpot01( uplo, n, a, lda, ainv, lda, rwork,
287  $ result( 1 ) )
288 *
289 *+ TEST 2
290 * Form the inverse and compute the residual.
291 *
292  CALL dlacpy( uplo, n, n, afac, lda, ainv, lda )
293  srnamt = 'DPOTRI'
294  CALL plasma_dpotri( plasma_uplo, n, ainv, lda,
295  $ info )
296 *
297 * Check error code from DPOTRI.
298 *
299  IF( info.NE.0 )
300  $ CALL alaerh( path, 'DPOTRI', info, 0, uplo, n, n,
301  $ -1, -1, -1, imat, nfail, nerrs, nout )
302 *
303  CALL dpot03( uplo, n, a, lda, ainv, lda, work, lda,
304  $ rwork, rcondc, result( 2 ) )
305 *
306 * Print information about the tests that did not pass
307 * the threshold.
308 *
309  DO 60 k = 1, 2
310  IF( result( k ).GE.thresh ) THEN
311  IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
312  $ CALL alahd( nout, path )
313  WRITE( nout, fmt = 9999 )uplo, n, nb, imat, k,
314  $ result( k )
315  nfail = nfail + 1
316  END IF
317  60 continue
318  nrun = nrun + 2
319 *
320 * Skip the rest of the tests unless this is the first
321 * blocksize.
322 *
323  IF( inb.NE.1 )
324  $ go to 90
325 *
326  DO 80 irhs = 1, nns
327  nrhs = nsval( irhs )
328 *
329 *+ TEST 3
330 * Solve and compute residual for A * X = B .
331 *
332  srnamt = 'DLARHS'
333  CALL dlarhs( path, xtype, uplo, ' ', n, n, kl, ku,
334  $ nrhs, a, lda, xact, lda, b, lda,
335  $ iseed, info )
336  CALL dlacpy( 'Full', n, nrhs, b, lda, x, lda )
337 *
338  srnamt = 'DPOTRS'
339  CALL plasma_dpotrs( plasma_uplo, n, nrhs, afac,
340  $ lda, x, lda, info )
341 *
342 *
343 * Check error code from DPOTRS.
344 *
345  IF( info.NE.0 )
346  $ CALL alaerh( path, 'DPOTRS', info, 0, uplo, n,
347  $ n, -1, -1, nrhs, imat, nfail,
348  $ nerrs, nout )
349 *
350  CALL dlacpy( 'Full', n, nrhs, b, lda, work, lda )
351  CALL dpot02( uplo, n, nrhs, a, lda, x, lda, work,
352  $ lda, rwork, result( 3 ) )
353 *+ TEST 4
354 * Check solution from generated exact solution.
355 *
356  CALL dget04( n, nrhs, x, lda, xact, lda, rcondc,
357  $ result( 4 ) )
358 *
359 *+ TESTS 5, 6, and 7
360 * Use iterative refinement to improve the solution.
361 *
362  srnamt = 'DPORFS'
363  CALL dporfs( uplo, n, nrhs, a, lda,
364  $ afac, lda, b, lda, x, lda,
365  $ rwork, rwork( nrhs+1 ), work, iwork,
366  $ info )
367 *
368 * Check error code from DPORFS.
369 *
370  IF( info.NE.0 )
371  $ CALL alaerh( path, 'DPORFS', info, 0, uplo, n,
372  $ n, -1, -1, nrhs, imat, nfail,
373  $ nerrs, nout )
374 *
375  CALL dget04( n, nrhs, x, lda, xact, lda, rcondc,
376  $ result( 5 ) )
377  CALL dpot05( uplo, n, nrhs, a, lda, b, lda, x, lda,
378  $ xact, lda, rwork, rwork( nrhs+1 ),
379  $ result( 6 ) )
380 *
381 * Print information about the tests that did not pass
382 * the threshold.
383 *
384  DO 70 k = 3, 7
385  IF( result( k ).GE.thresh ) THEN
386  IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
387  $ CALL alahd( nout, path )
388  WRITE( nout, fmt = 9998 )uplo, n, nrhs,
389  $ imat, k, result( k )
390  nfail = nfail + 1
391  END IF
392  70 continue
393  nrun = nrun + 5
394  80 continue
395 *
396 *+ TEST 8
397 * Get an estimate of RCOND = 1/CNDNUM.
398 *
399  anorm = dlansy( '1', uplo, n, a, lda, rwork )
400  srnamt = 'DPOCON'
401  CALL dpocon( uplo, n, afac, lda, anorm, rcond, work,
402  $ iwork, info )
403 *
404 * Check error code from DPOCON.
405 *
406  IF( info.NE.0 )
407  $ CALL alaerh( path, 'DPOCON', info, 0, uplo, n, n,
408  $ -1, -1, -1, imat, nfail, nerrs, nout )
409 *
410  result( 8 ) = dget06( rcond, rcondc )
411 *
412 * Print the test ratio if it is .GE. THRESH.
413 *
414  IF( result( 8 ).GE.thresh ) THEN
415  IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
416  $ CALL alahd( nout, path )
417  WRITE( nout, fmt = 9997 )uplo, n, imat, 8,
418  $ result( 8 )
419  nfail = nfail + 1
420  END IF
421  nrun = nrun + 1
422  90 continue
423  100 continue
424  110 continue
425  120 continue
426 *
427 * Print a summary of the results.
428 *
429  CALL alasum( path, nout, nfail, nrun, nerrs )
430 *
431  9999 format( ' UPLO = ''', a1, ''', N =', i5, ', NB =', i4, ', type ',
432  $ i2, ', test ', i2, ', ratio =', g12.5 )
433  9998 format( ' UPLO = ''', a1, ''', N =', i5, ', NRHS=', i3, ', type ',
434  $ i2, ', test(', i2, ') =', g12.5 )
435  9997 format( ' UPLO = ''', a1, ''', N =', i5, ',', 10x, ' type ', i2,
436  $ ', test(', i2, ') =', g12.5 )
437  return
438 *
439 * End of DCHKPO
440 *
441  END