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
ddrvls.f
Go to the documentation of this file.
1  SUBROUTINE ddrvls( DOTYPE, NM, MVAL, NN, NVAL, NNS, NSVAL, NNB,
2  $ nbval, nxval, thresh, tsterr, a, copya, b,
3  $ copyb, c, s, copys, ibval, work, iwork, nout )
4 *
5  include 'plasmaf.h'
6 *
7 * -- LAPACK test routine (version 3.1.1) --
8 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
9 * January 2007
10 *
11 * .. Scalar Arguments ..
12  LOGICAL tsterr
13  INTEGER nm, nn, nnb, nns, nout
14  DOUBLE PRECISION thresh
15 * ..
16 * .. Array Arguments ..
17  LOGICAL dotype( * )
18  INTEGER iwork( * ), mval( * ), nbval( * ), nsval( * ),
19  $ nval( * ), nxval( * ), ibval( * )
20  DOUBLE PRECISION a( * ), b( * ), c( * ), copya( * ), copyb( * ),
21  $ copys( * ), s( * ), work( * )
22 * ..
23 *
24 * Purpose
25 * =======
26 *
27 * DDRVLS tests the least squares driver routines DGELS, DGELSS, DGELSX,
28 * DGELSY and DGELSD.
29 *
30 * Arguments
31 * =========
32 *
33 * DOTYPE (input) LOGICAL array, dimension (NTYPES)
34 * The matrix types to be used for testing. Matrices of type j
35 * (for 1 <= j <= NTYPES) are used for testing if DOTYPE(j) =
36 * .TRUE.; if DOTYPE(j) = .FALSE., then type j is not used.
37 * The matrix of type j is generated as follows:
38 * j=1: A = U*D*V where U and V are random orthogonal matrices
39 * and D has random entries (> 0.1) taken from a uniform
40 * distribution (0,1). A is full rank.
41 * j=2: The same of 1, but A is scaled up.
42 * j=3: The same of 1, but A is scaled down.
43 * j=4: A = U*D*V where U and V are random orthogonal matrices
44 * and D has 3*min(M,N)/4 random entries (> 0.1) taken
45 * from a uniform distribution (0,1) and the remaining
46 * entries set to 0. A is rank-deficient.
47 * j=5: The same of 4, but A is scaled up.
48 * j=6: The same of 5, but A is scaled down.
49 *
50 * NM (input) INTEGER
51 * The number of values of M contained in the vector MVAL.
52 *
53 * MVAL (input) INTEGER array, dimension (NM)
54 * The values of the matrix row dimension M.
55 *
56 * NN (input) INTEGER
57 * The number of values of N contained in the vector NVAL.
58 *
59 * NVAL (input) INTEGER array, dimension (NN)
60 * The values of the matrix column dimension N.
61 *
62 * NNS (input) INTEGER
63 * The number of values of NRHS contained in the vector NSVAL.
64 *
65 * NSVAL (input) INTEGER array, dimension (NNS)
66 * The values of the number of right hand sides NRHS.
67 *
68 * NNB (input) INTEGER
69 * The number of values of NB and NX contained in the
70 * vectors NBVAL and NXVAL. The blocking parameters are used
71 * in pairs (NB,NX).
72 *
73 * NBVAL (input) INTEGER array, dimension (NNB)
74 * The values of the tile size NB.
75 *
76 * IBVAL (input) INTEGER array, dimension (NNB)
77 * The values of the inner block size IB.
78 *
79 * NXVAL (input) INTEGER array, dimension (NNB)
80 * The values of the crossover point NX.
81 *
82 * THRESH (input) DOUBLE PRECISION
83 * The threshold value for the test ratios. A result is
84 * included in the output file if RESULT >= THRESH. To have
85 * every test ratio printed, use THRESH = 0.
86 *
87 * TSTERR (input) LOGICAL
88 * Flag that indicates whether error exits are to be tested.
89 *
90 * A (workspace) DOUBLE PRECISION array, dimension (MMAX*NMAX)
91 * where MMAX is the maximum value of M in MVAL and NMAX is the
92 * maximum value of N in NVAL.
93 *
94 * COPYA (workspace) DOUBLE PRECISION array, dimension (MMAX*NMAX)
95 *
96 * B (workspace) DOUBLE PRECISION array, dimension (MMAX*NSMAX)
97 * where MMAX is the maximum value of M in MVAL and NSMAX is the
98 * maximum value of NRHS in NSVAL.
99 *
100 * COPYB (workspace) DOUBLE PRECISION array, dimension (MMAX*NSMAX)
101 *
102 * C (workspace) DOUBLE PRECISION array, dimension (MMAX*NSMAX)
103 *
104 * S (workspace) DOUBLE PRECISION array, dimension
105 * (min(MMAX,NMAX))
106 *
107 * COPYS (workspace) DOUBLE PRECISION array, dimension
108 * (min(MMAX,NMAX))
109 *
110 * WORK (workspace) DOUBLE PRECISION array,
111 * dimension (MMAX*NMAX + 4*NMAX + MMAX).
112 *
113 * IWORK (workspace) INTEGER array, dimension (15*NMAX)
114 *
115 * NOUT (input) INTEGER
116 * The unit number for output.
117 *
118 * =====================================================================
119 *
120 * .. Parameters ..
121  INTEGER ntests
122  parameter( ntests = 18 )
123  INTEGER smlsiz
124  parameter( smlsiz = 25 )
125  DOUBLE PRECISION one, two, zero
126  parameter( one = 1.0d0, two = 2.0d0, zero = 0.0d0 )
127 * ..
128 * .. Local Scalars ..
129  CHARACTER trans
130  CHARACTER*3 path
131  INTEGER crank, i, im, in, inb, info, ins, irank,
132  $ iscale, itran, itype, j, k, lda, ldb, ldwork,
133  $ lwlsy, lwork, m, mnmin, n, nb, ib, ncols,
134  $ nerrs, nfail, nlvl, nrhs, nrows, nrun, rank,
135  $ plasma_trans
136  INTEGER ht( 2 )
137  DOUBLE PRECISION eps, norma, normb, rcond
138 * ..
139 * .. Local Arrays ..
140  INTEGER iseed( 4 ), iseedy( 4 )
141  DOUBLE PRECISION result( ntests )
142 * ..
143 * .. External Functions ..
144  DOUBLE PRECISION dasum, dlamch, dqrt14, dqrt17
145  EXTERNAL dasum, dlamch, dqrt14, dqrt17
146 * ..
147 * .. External Subroutines ..
148  EXTERNAL alaerh, alahd, alasvm, daxpy, derrls, dgels,
149  $ dgelsd, dgelss, dgelsx, dgelsy, dgemm, dlacpy,
150  $ dlarnv, dlasrt, dqrt13, dqrt15, dqrt16, dscal,
151  $ xlaenv
152 * ..
153 * .. Intrinsic Functions ..
154  INTRINSIC dble, int, log, max, min, sqrt
155 * ..
156 * .. Scalars in Common ..
157  LOGICAL lerr, ok
158  CHARACTER*32 srnamt
159  INTEGER infot, iounit
160 * ..
161 * .. Common blocks ..
162  common / infoc / infot, iounit, ok, lerr
163  common / srnamc / srnamt
164 * ..
165 * .. Data statements ..
166  DATA iseedy / 1988, 1989, 1990, 1991 /
167 * ..
168 * .. Executable Statements ..
169 *
170 * Initialize constants and the random number seed.
171 *
172  path( 1: 1 ) = 'Double precision'
173  path( 2: 3 ) = 'LS'
174  nrun = 0
175  nfail = 0
176  nerrs = 0
177  DO 10 i = 1, 4
178  iseed( i ) = iseedy( i )
179  10 continue
180  eps = dlamch( 'Epsilon' )
181 *
182 * Threshold for rank estimation
183 *
184  rcond = sqrt( eps ) - ( sqrt( eps )-eps ) / 2
185 *
186 * Test the error exits
187 *
188  CALL xlaenv( 2, 2 )
189  CALL xlaenv( 9, smlsiz )
190  IF( tsterr )
191  $ CALL derrls( path, nout )
192 *
193 * Print the header if NM = 0 or NN = 0 and THRESH = 0.
194 *
195  IF( ( nm.EQ.0 .OR. nn.EQ.0 ) .AND. thresh.EQ.zero )
196  $ CALL alahd( nout, path )
197  infot = 0
198  CALL xlaenv( 2, 2 )
199  CALL xlaenv( 9, smlsiz )
200 *
201  DO 150 im = 1, nm
202  m = mval( im )
203  lda = max( 1, m )
204 *
205  DO 140 in = 1, nn
206  n = nval( in )
207  mnmin = min( m, n )
208  ldb = max( 1, m, n )
209 *
210  DO 130 ins = 1, nns
211  nrhs = nsval( ins )
212  nlvl = max( int( log( max( one, dble( mnmin ) ) /
213  $ dble( smlsiz+1 ) ) / log( two ) ) + 1, 0 )
214  lwork = max( 1, ( m+nrhs )*( n+2 ), ( n+nrhs )*( m+2 ),
215  $ m*n+4*mnmin+max( m, n ), 12*mnmin+2*mnmin*smlsiz+
216  $ 8*mnmin*nlvl+mnmin*nrhs+(smlsiz+1)**2 )
217 *
218  DO 120 irank = 1, 2
219  DO 110 iscale = 1, 3
220  itype = ( irank-1 )*3 + iscale
221  IF( .NOT.dotype( itype ) )
222  $ go to 110
223 *
224  IF( irank.EQ.1 ) THEN
225 *
226 * Test DGELS
227 *
228 * Generate a matrix of scaling type ISCALE
229 *
230  CALL dqrt13( iscale, m, n, copya, lda, norma,
231  $ iseed )
232  DO 40 inb = 1, nnb
233  nb = nbval( inb )
234  ib = ibval( inb )
235  CALL xlaenv( 1, nb )
236  CALL xlaenv( 3, nxval( inb ) )
237  IF ( (max(m, n) / 25) .GT. nb ) THEN
238  goto 40
239  END IF
240  CALL plasma_set( plasma_tile_size, nb, info )
241  CALL plasma_set( plasma_inner_block_size, ib,
242  $ info )
243 *
244 * Allocate T
245 *
246  CALL plasma_alloc_workspace_dgels( m, n , ht,
247  $ info )
248 *
249 * DO 30 ITRAN = 1, 2
250  DO 30 itran = 1, 1
251 *
252 * ONLY PLASMANOTRANS supported !
253 *
254  IF( itran.EQ.1 ) THEN
255  trans = 'N'
256  plasma_trans = plasmanotrans
257  nrows = m
258  ncols = n
259  ELSE
260  trans = 'T'
261  plasma_trans = plasmatrans
262  nrows = n
263  ncols = m
264  END IF
265  ldwork = max( 1, ncols )
266 *
267 * Set up a consistent rhs
268 *
269  IF( ncols.GT.0 ) THEN
270  CALL dlarnv( 2, iseed, ncols*nrhs,
271  $ work )
272  CALL dscal( ncols*nrhs,
273  $ one / dble( ncols ), work,
274  $ 1 )
275  END IF
276  CALL dgemm( trans, 'No transpose', nrows,
277  $ nrhs, ncols, one, copya, lda,
278  $ work, ldwork, zero, b, ldb )
279  CALL dlacpy( 'Full', nrows, nrhs, b, ldb,
280  $ copyb, ldb )
281 *
282 * Solve LS or overdetermined system
283 *
284  IF( m.GT.0 .AND. n.GT.0 ) THEN
285  CALL dlacpy( 'Full', m, n, copya, lda,
286  $ a, lda )
287  CALL dlacpy( 'Full', nrows, nrhs,
288  $ copyb, ldb, b, ldb )
289  END IF
290  srnamt = 'DGELS '
291 *
292  CALL plasma_dgels( plasma_trans,
293  $ m, n, nrhs,
294  $ a, lda, ht, b, ldb,
295  $ info )
296 *
297  IF( info.NE.0 )
298  $ CALL alaerh( path, 'DGELS ', info, 0,
299  $ trans, m, n, nrhs, -1, nb,
300  $ itype, nfail, nerrs,
301  $ nout )
302 *
303 * Check correctness of results
304 *
305  ldwork = max( 1, nrows )
306  IF( nrows.GT.0 .AND. nrhs.GT.0 )
307  $ CALL dlacpy( 'Full', nrows, nrhs,
308  $ copyb, ldb, c, ldb )
309  CALL dqrt16( trans, m, n, nrhs, copya,
310  $ lda, b, ldb, c, ldb, work,
311  $ result( 1 ) )
312 *
313  IF( ( itran.EQ.1 .AND. m.GE.n ) .OR.
314  $ ( itran.EQ.2 .AND. m.LT.n ) ) THEN
315 *
316 * Solving LS system
317 *
318  result( 2 ) = dqrt17( trans, 1, m, n,
319  $ nrhs, copya, lda, b, ldb,
320  $ copyb, ldb, c, work,
321  $ lwork )
322  ELSE
323 *
324 * Solving overdetermined system
325 *
326  result( 2 ) = dqrt14( trans, m, n,
327  $ nrhs, copya, lda, b, ldb,
328  $ work, lwork )
329  END IF
330 *
331 * Print information about the tests that
332 * did not pass the threshold.
333 *
334  DO 20 k = 1, 2
335  IF( result( k ).GE.thresh ) THEN
336  IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
337  $ CALL alahd( nout, path )
338  WRITE( nout, fmt = 9999 )trans, m,
339  $ n, nrhs, nb, itype, k,
340  $ result( k )
341  nfail = nfail + 1
342  END IF
343  20 continue
344  nrun = nrun + 2
345  30 continue
346 *
347 * Deallocate T
348 *
349  CALL plasma_dealloc_handle( ht, info )
350  40 continue
351  END IF
352 *
353  110 continue
354  120 continue
355  130 continue
356  140 continue
357  150 continue
358 *
359 * Print a summary of the results.
360 *
361  CALL alasvm( path, nout, nfail, nrun, nerrs )
362 *
363  9999 format( ' TRANS=''', a1, ''', M=', i5, ', N=', i5, ', NRHS=', i4,
364  $ ', NB=', i4, ', type', i2, ', test(', i2, ')=', g12.5 )
365  9998 format( ' M=', i5, ', N=', i5, ', NRHS=', i4, ', NB=', i4,
366  $ ', type', i2, ', test(', i2, ')=', g12.5 )
367  return
368 *
369 * End of DDRVLS
370 *
371  END