MAGMA  1.2.0
MatrixAlgebraonGPUandMulticoreArchitectures
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Groups
ddrvac.f
Go to the documentation of this file.
1  SUBROUTINE ddrvac( DOTYPE, NM, MVAL, NNS, NSVAL, THRESH, NMAX,
2  $ a, afac, b, x, work,
3  $ rwork, swork, nout )
4 *
5 * -- LAPACK test routine (version 3.1.2) --
6 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
7 * April 2007
8 *
9 * .. Scalar Arguments ..
10  INTEGER nmax, nm, nns, nout
11  DOUBLE PRECISION thresh
12 * ..
13 * .. Array Arguments ..
14  LOGICAL dotype( * )
15  INTEGER mval( * ), nsval( * )
16  REAL swork(*)
17  DOUBLE PRECISION a( * ), afac( * ), b( * ),
18  $ rwork( * ), work( * ), x( * )
19 * ..
20 *
21 * Purpose
22 * =======
23 *
24 * DDRVAC tests DSPOSV.
25 *
26 * Arguments
27 * =========
28 *
29 * DOTYPE (input) LOGICAL array, dimension (NTYPES)
30 * The matrix types to be used for testing. Matrices of type j
31 * (for 1 <= j <= NTYPES) are used for testing if DOTYPE(j) =
32 * .TRUE.; if DOTYPE(j) = .FALSE., then type j is not used.
33 *
34 * NM (input) INTEGER
35 * The number of values of N contained in the vector MVAL.
36 *
37 * MVAL (input) INTEGER array, dimension (NM)
38 * The values of the matrix dimension N.
39 *
40 * NNS (input) INTEGER
41 * The number of values of NRHS contained in the vector NSVAL.
42 *
43 * NSVAL (input) INTEGER array, dimension (NNS)
44 * The values of the number of right hand sides NRHS.
45 *
46 * THRESH (input) DOUBLE PRECISION
47 * The threshold value for the test ratios. A result is
48 * included in the output file if RESULT >= THRESH. To have
49 * every test ratio printed, use THRESH = 0.
50 *
51 * NMAX (input) INTEGER
52 * The maximum value permitted for N, used in dimensioning the
53 * work arrays.
54 *
55 * A (workspace) DOUBLE PRECISION array, dimension (NMAX*NMAX)
56 *
57 * AFAC (workspace) DOUBLE PRECISION array, dimension (NMAX*NMAX)
58 *
59 * B (workspace) DOUBLE PRECISION array, dimension (NMAX*NSMAX)
60 *
61 * X (workspace) DOUBLE PRECISION array, dimension (NMAX*NSMAX)
62 *
63 * WORK (workspace) DOUBLE PRECISION array, dimension
64 * (NMAX*max(3,NSMAX))
65 *
66 * RWORK (workspace) DOUBLE PRECISION array, dimension
67 * (max(2*NMAX,2*NSMAX+NWORK))
68 *
69 * SWORK (workspace) REAL array, dimension
70 * (NMAX*(NSMAX+NMAX))
71 *
72 * NOUT (input) INTEGER
73 * The unit number for output.
74 *
75 * =====================================================================
76 *
77 * .. Parameters ..
78  DOUBLE PRECISION one, zero
79  parameter( one = 1.0d+0, zero = 0.0d+0 )
80  INTEGER ntypes
81  parameter( ntypes = 9 )
82  INTEGER ntests
83  parameter( ntests = 1 )
84 * ..
85 * .. Local Scalars ..
86  LOGICAL zerot
87  CHARACTER dist, type, uplo, xtype
88  CHARACTER*3 path
89  INTEGER i, im, imat, info, ioff, irhs, iuplo,
90  $ izero, kl, ku, lda, mode, n,
91  $ nerrs, nfail, nimat, nrhs, nrun
92  DOUBLE PRECISION anorm, cndnum
93 * ..
94 * .. Local Arrays ..
95  CHARACTER uplos( 2 )
96  INTEGER iseed( 4 ), iseedy( 4 )
97  DOUBLE PRECISION result( ntests )
98 * ..
99 * .. Local Variables ..
100  INTEGER iter, kase
101 * ..
102 * .. External Functions ..
103  LOGICAL lsame
104  EXTERNAL lsame
105 * ..
106 * .. External Subroutines ..
107  EXTERNAL alaerh, dlacpy,
108  $ dlarhs, dlaset, dlatb4, dlatms,
109  $ dpot06, dsposv
110 * ..
111 * .. Intrinsic Functions ..
112  INTRINSIC dble, max, sqrt
113 * ..
114 * .. Scalars in Common ..
115  LOGICAL lerr, ok
116  CHARACTER*32 srnamt
117  INTEGER infot, nunit
118 * ..
119 * .. Common blocks ..
120  common / infoc / infot, nunit, ok, lerr
121  common / srnamc / srnamt
122 * ..
123 * .. Data statements ..
124  DATA iseedy / 1988, 1989, 1990, 1991 /
125  DATA uplos / 'U', 'L' /
126 * ..
127 * .. Executable Statements ..
128 *
129 * Initialize constants and the random number seed.
130 *
131  kase = 0
132  path( 1: 1 ) = 'Double precision'
133  path( 2: 3 ) = 'PO'
134  nrun = 0
135  nfail = 0
136  nerrs = 0
137  DO 10 i = 1, 4
138  iseed( i ) = iseedy( i )
139  10 continue
140 *
141  infot = 0
142 *
143 * Do for each value of N in MVAL
144 *
145  DO 120 im = 1, nm
146  n = mval( im )
147  lda = max( n, 1 )
148  nimat = ntypes
149  IF( n.LE.0 )
150  $ nimat = 1
151 *
152  DO 110 imat = 1, nimat
153 *
154 * Do the tests only if DOTYPE( IMAT ) is true.
155 *
156  IF( .NOT.dotype( imat ) )
157  $ go to 110
158 *
159 * Skip types 3, 4, or 5 if the matrix size is too small.
160 *
161  zerot = imat.GE.3 .AND. imat.LE.5
162  IF( zerot .AND. n.LT.imat-2 )
163  $ go to 110
164 *
165 * Do first for UPLO = 'U', then for UPLO = 'L'
166 *
167  DO 100 iuplo = 1, 2
168  uplo = uplos( iuplo )
169 *
170 * Set up parameters with DLATB4 and generate a test matrix
171 * with DLATMS.
172 *
173  CALL dlatb4( path, imat, n, n, type, kl, ku, anorm, mode,
174  $ cndnum, dist )
175 *
176  srnamt = 'DLATMS'
177  CALL dlatms( n, n, dist, iseed, type, rwork, mode,
178  $ cndnum, anorm, kl, ku, uplo, a, lda, work,
179  $ info )
180 *
181 * Check error code from DLATMS.
182 *
183  IF( info.NE.0 ) THEN
184  CALL alaerh( path, 'DLATMS', info, 0, uplo, n, n, -1,
185  $ -1, -1, imat, nfail, nerrs, nout )
186  go to 100
187  END IF
188 *
189 * For types 3-5, zero one row and column of the matrix to
190 * test that INFO is returned correctly.
191 *
192  IF( zerot ) THEN
193  IF( imat.EQ.3 ) THEN
194  izero = 1
195  ELSE IF( imat.EQ.4 ) THEN
196  izero = n
197  ELSE
198  izero = n / 2 + 1
199  END IF
200  ioff = ( izero-1 )*lda
201 *
202 * Set row and column IZERO of A to 0.
203 *
204  IF( iuplo.EQ.1 ) THEN
205  DO 20 i = 1, izero - 1
206  a( ioff+i ) = zero
207  20 continue
208  ioff = ioff + izero
209  DO 30 i = izero, n
210  a( ioff ) = zero
211  ioff = ioff + lda
212  30 continue
213  ELSE
214  ioff = izero
215  DO 40 i = 1, izero - 1
216  a( ioff ) = zero
217  ioff = ioff + lda
218  40 continue
219  ioff = ioff - izero
220  DO 50 i = izero, n
221  a( ioff+i ) = zero
222  50 continue
223  END IF
224  ELSE
225  izero = 0
226  END IF
227 *
228  DO 60 irhs = 1, nns
229  nrhs = nsval( irhs )
230  xtype = 'N'
231 *
232 * Form an exact solution and set the right hand side.
233 *
234  srnamt = 'DLARHS'
235  CALL dlarhs( path, xtype, uplo, ' ', n, n, kl, ku,
236  $ nrhs, a, lda, x, lda, b, lda,
237  $ iseed, info )
238 *
239 * Compute the L*L' or U'*U factorization of the
240 * matrix and solve the system.
241 *
242  srnamt = 'DSPOSV '
243  kase = kase + 1
244 *
245  CALL dlacpy( 'All', n, n, a, lda, afac, lda)
246 *
247  CALL dsposv( uplo, n, nrhs, afac, lda, b, lda, x, lda,
248  $ work, swork, iter, info )
249 
250  IF (iter.LT.0) THEN
251  CALL dlacpy( 'All', n, n, a, lda, afac, lda )
252  ENDIF
253 *
254 * Check error code from DSPOSV .
255 *
256  IF( info.NE.izero ) THEN
257 *
258  IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
259  $ CALL alahd( nout, path )
260  nerrs = nerrs + 1
261 *
262  IF( info.NE.izero .AND. izero.NE.0 ) THEN
263  WRITE( nout, fmt = 9988 )'DSPOSV',info,izero,n,
264  $ imat
265  ELSE
266  WRITE( nout, fmt = 9975 )'DSPOSV',info,n,imat
267  END IF
268  END IF
269 *
270 * Skip the remaining test if the matrix is singular.
271 *
272  IF( info.NE.0 )
273  $ go to 110
274 *
275 * Check the quality of the solution
276 *
277  CALL dlacpy( 'All', n, nrhs, b, lda, work, lda )
278 *
279  CALL dpot06( uplo, n, nrhs, a, lda, x, lda, work,
280  $ lda, rwork, result( 1 ) )
281 *
282 * Check if the test passes the tesing.
283 * Print information about the tests that did not
284 * pass the testing.
285 *
286 * If iterative refinement has been used and claimed to
287 * be successful (ITER>0), we want
288 * NORM1(B - A*X)/(NORM1(A)*NORM1(X)*EPS*SRQT(N)) < 1
289 *
290 * If double precision has been used (ITER<0), we want
291 * NORM1(B - A*X)/(NORM1(A)*NORM1(X)*EPS) < THRES
292 * (Cf. the linear solver testing routines)
293 *
294  IF ((thresh.LE.0.0e+00)
295  $ .OR.((iter.GE.0).AND.(n.GT.0)
296  $ .AND.(result(1).GE.sqrt(dble(n))))
297  $ .OR.((iter.LT.0).AND.(result(1).GE.thresh))) THEN
298 *
299  IF( nfail.EQ.0 .AND. nerrs.EQ.0 ) THEN
300  WRITE( nout, fmt = 8999 )'DPO'
301  WRITE( nout, fmt = '( '' Matrix types:'' )' )
302  WRITE( nout, fmt = 8979 )
303  WRITE( nout, fmt = '( '' Test ratios:'' )' )
304  WRITE( nout, fmt = 8960 )1
305  WRITE( nout, fmt = '( '' Messages:'' )' )
306  END IF
307 *
308  WRITE( nout, fmt = 9998 )uplo, n, nrhs, imat, 1,
309  $ result( 1 )
310 *
311  nfail = nfail + 1
312 *
313  END IF
314 *
315  nrun = nrun + 1
316 *
317  60 continue
318  100 continue
319  110 continue
320  120 continue
321 *
322  130 continue
323 *
324 * Print a summary of the results.
325 *
326  IF( nfail.GT.0 ) THEN
327  WRITE( nout, fmt = 9996 )'DSPOSV', nfail, nrun
328  ELSE
329  WRITE( nout, fmt = 9995 )'DSPOSV', nrun
330  END IF
331  IF( nerrs.GT.0 ) THEN
332  WRITE( nout, fmt = 9994 )nerrs
333  END IF
334 *
335  9998 format( ' UPLO=''', a1, ''', N =', i5, ', NRHS=', i3, ', type ',
336  $ i2, ', test(', i2, ') =', g12.5 )
337  9996 format( 1x, a6, ': ', i6, ' out of ', i6,
338  $ ' tests failed to pass the threshold' )
339  9995 format( /1x, 'All tests for ', a6,
340  $ ' routines passed the threshold (', i6, ' tests run)' )
341  9994 format( 6x, i6, ' error messages recorded' )
342 *
343 * SUBNAM, INFO, INFOE, N, IMAT
344 *
345  9988 format( ' *** ', a6, ' returned with INFO =', i5, ' instead of ',
346  $ i5, / ' ==> N =', i5, ', type ',
347  $ i2 )
348 *
349 * SUBNAM, INFO, N, IMAT
350 *
351  9975 format( ' *** Error code from ', a6, '=', i5, ' for M=', i5,
352  $ ', type ', i2 )
353  8999 format( / 1x, a3, ': positive definite dense matrices' )
354  8979 format( 4x, '1. Diagonal', 24x, '7. Last n/2 columns zero', / 4x,
355  $ '2. Upper triangular', 16x,
356  $ '8. Random, CNDNUM = sqrt(0.1/EPS)', / 4x,
357  $ '3. Lower triangular', 16x, '9. Random, CNDNUM = 0.1/EPS',
358  $ / 4x, '4. Random, CNDNUM = 2', 13x,
359  $ '10. Scaled near underflow', / 4x, '5. First column zero',
360  $ 14x, '11. Scaled near overflow', / 4x,
361  $ '6. Last column zero' )
362  8960 format( 3x, i2, ': norm_1( B - A * X ) / ',
363  $ '( norm_1(A) * norm_1(X) * EPS * SQRT(N) ) > 1 if ITERREF',
364  $ / 4x, 'or norm_1( B - A * X ) / ',
365  $ '( norm_1(A) * norm_1(X) * EPS ) > THRES if DPOTRF' )
366 
367  return
368 *
369 * End of DDRVAC
370 *
371  END