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