3
4
5
6
7
8
9
10 INTEGER LWORK, M, N, NB, NOUT, NPCOL, NPROW
11 DOUBLE PRECISION THRESH
12
13
14 INTEGER ISEED( 4 ), RESULT( 9 )
15 DOUBLE PRECISION WORK( * )
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210 INTEGER BLOCK_CYCLIC_2D, DLEN_, DTYPE_, CTXT_, M_, N_,
211 $ MB_, NB_, RSRC_, CSRC_, LLD_, NTYPES
212 parameter( block_cyclic_2d = 1, dlen_ = 9, dtype_ = 1,
213 $ ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
214 $ rsrc_ = 7, csrc_ = 8, lld_ = 9, ntypes = 6 )
215 DOUBLE PRECISION ZERO, ONE
216 parameter( zero = 0.0d0, one = 1.0d0 )
217
218
219 CHARACTER HETERO, JOBU, JOBVT
220 INTEGER CONTEXT, DINFO, I, IA, IAM, INFO, ITYPE, IU,
221 $ IVT, JA, JOBTYPE, JU, JVT, LDA, LDU, LDVT,
222 $ LLWORK, LWMIN, MYCOL, MYROW, NNODES, NQ, PASS,
223 $ PTRA, PTRAC, PTRD, PTRWORK, PTRS, PTRSC, PTRU,
224 $ PTRUC, PTRVT, PTRVTC, SETHET, SIZE, SIZEQ,
225 $ WPDGESVD, WPDLAGGE, WPDSVDCHK, WPDSVDCMP
226 DOUBLE PRECISION CHK, DELTA, H, MTM, OVFL, RTOVFL, RTUNFL, ULP,
227 $ UNFL
228
229
230 EXTERNAL blacs_barrier, blacs_get, blacs_gridexit,
231 $ blacs_gridinfo, blacs_gridinit, blacs_pinfo,
232 $ blacs_set,
233 $
descinit, dgamn2d, dgamx2d, dlabad, dscal,
234 $ igamn2d, igamx2d, igebr2d, igebs2d,
pdelset,
237
238
239 INTEGER NUMROC
240 DOUBLE PRECISION PDLAMCH
242
243
244 INTEGER DESCA( DLEN_ ), DESCU( DLEN_ ),
245 $ DESCVT( DLEN_ ), ITMP( 2 )
246 DOUBLE PRECISION CTIME( 1 ), WTIME( 1 )
247
248
249 INTRINSIC abs, int,
max,
min, sqrt
250
251
252
253 IF( block_cyclic_2d*csrc_*dtype_*lld_*mb_*m_*nb_*n_*rsrc_.LT.0 )
254 $ RETURN
255
256 CALL blacs_pinfo( iam, nnodes )
257 CALL blacs_get( -1, 0, context )
258 CALL blacs_gridinit( context, 'R', nprow, npcol )
259 CALL blacs_gridinfo( context, nprow, npcol, myrow, mycol )
260
261
262
263 IF( ( myrow.GE.nprow ) .OR. ( myrow.LT.0 ) .OR.
264 $ ( mycol.GE.npcol ) .OR. ( mycol.LT.0 ) )GO TO 110
265 CALL blacs_set( context, 15, 1 )
266 info = 0
267
268
269
270 IF( m.LE.0 ) THEN
271 info = -1
272 ELSE IF( n.LE.0 ) THEN
273 info = -2
274 ELSE IF( nprow.LE.0 ) THEN
275 info = -3
276 ELSE IF( npcol.LE.0 ) THEN
277 info = -4
278 ELSE IF( nb.LE.0 ) THEN
279 info = -5
280 ELSE IF( thresh.LE.0 ) THEN
281 info = -7
282 END IF
283
285
286
287
288 ia = 1
289 ja = 1
290 iu = 1
291 ju = 1
292 ivt = 1
293 jvt = 1
294
295 lda =
numroc( m, nb, myrow, 0, nprow )
297 nq =
numroc( n, nb, mycol, 0, npcol )
298 ldu = lda
299 sizeq =
numroc(
SIZE, nb, mycol, 0, npcol )
300 ldvt =
numroc(
SIZE, nb, myrow, 0, nprow )
301 ldvt =
max( 1, ldvt )
302 CALL descinit( desca, m, n, nb, nb, 0, 0, context, lda, dinfo )
303 CALL descinit( descu, m,
SIZE, nb, nb, 0, 0, context, ldu, dinfo )
304 CALL descinit( descvt,
SIZE, n, nb, nb, 0, 0, context, ldvt,
305 $ dinfo )
306
307
308
309 ptra = 2
310 ptrac = ptra + lda*nq
311 ptrd = ptrac + lda*nq
312 ptrs = ptrd + SIZE
313 ptrsc = ptrs + SIZE
314 ptrwork = ptrsc + SIZE
315
316 ptru = ptrwork
317 ptrvt = ptrwork
318 ptruc = ptrwork
319 ptrvtc = ptrwork
320
321
322
323
324 CALL pdlagge( m, n, work( ptrd ), work( ptra ), ia, ja, desca,
325 $ iseed, SIZE, work( ptrwork ), -1, dinfo )
326 wpdlagge = int( work( ptrwork ) )
327
328 CALL pdgesvd(
'V',
'V', m, n, work( ptra ), ia, ja, desca,
329 $ work( ptrs ), work( ptru ), iu, ju, descu,
330 $ work( ptrvt ), ivt, jvt, descvt,
331 $ work( ptrwork ), -1, dinfo )
332 wpdgesvd = int( work( ptrwork ) )
333
334 CALL pdsvdchk( m, n, work( ptrac ), ia, ja, desca, work( ptruc ),
335 $ iu, ju, descu, work( ptrvt ), ivt, jvt, descvt,
336 $ work( ptrs ), thresh, work( ptrwork ), -1,
337 $ result, chk, mtm )
338 wpdsvdchk = int( work( ptrwork ) )
339
340 CALL pdsvdcmp( m, n, 1, work( ptrs ), work( ptrsc ), work( ptru ),
341 $ work( ptruc ), iu, ju, descu, work( ptrvt ),
342 $ work( ptrvtc ), ivt, jvt, descvt, thresh,
343 $ result, delta, work( ptrwork ), -1 )
344 wpdsvdcmp = int( work( ptrwork ) )
345
346
347
348 lwmin = 1 + 2*lda*nq + 3*SIZE +
349 $
max( wpdlagge, ldu*sizeq+ldvt*nq+
max( ldu*sizeq,
350 $ ldvt*nq )+wpdgesvd+
max( wpdsvdchk, wpdsvdcmp ) )
351 work( 1 ) = lwmin
352
353
354
355 IF( lwork.EQ.-1 )
356 $ GO TO 120
357 IF( info.EQ.0 ) THEN
358 IF( lwork.LT.lwmin ) THEN
359 info = -10
360 END IF
361 END IF
362
363 IF( info.NE.0 ) THEN
364 CALL pxerbla( desca( ctxt_ ),
'PDSVDTST', -info )
365 RETURN
366 END IF
367
369 unfl =
pdlamch( context,
'Safe min' )
370 ovfl = one / unfl
371 CALL dlabad( unfl, ovfl )
372 rtunfl = sqrt( unfl )
373 rtovfl = sqrt( ovfl )
374
375
376
377 IF( myrow.EQ.0 .AND. mycol.EQ.0 ) THEN
378 CALL igebs2d( context, 'a', ' ', 4, 1, iseed, 4 )
379 ELSE
380 CALL igebr2d( context, 'a', ' ', 4, 1, iseed, 4, 0, 0 )
381 END IF
382
383
384
385 DO 100 itype = 1, ntypes
386
387 pass = 0
388 sethet = 0
389 ptrwork = ptrsc + SIZE
390 llwork = lwork - ptrwork + 1
391
392
393
394 IF( itype.EQ.1 ) THEN
395
396
397
398 DO 10 i = 1, SIZE
399 work( ptrd+i-1 ) = zero
400 10 CONTINUE
401
402 CALL pdlaset(
'All', m, n, zero, zero, work( ptra ),
403 $ ia, ja, desca )
404
405 ELSE IF( itype.EQ.2 ) THEN
406
407
408
409 DO 20 i = 1, SIZE
410 work( ptrd+i-1 ) = one
411 20 CONTINUE
412
413 CALL pdlaset(
'All', m, n, zero, one, work( ptra ),
414 $ ia, ja, desca )
415
416 ELSE IF( itype.GT.2 ) THEN
417
418
419
420 IF( size.NE.1 ) THEN
421 h = ( ulp-1 ) / ( size-1 )
422 DO 30 i = 1, SIZE
423 work( ptrd+i-1 ) = 1 + h*( i-1 )
424 30 CONTINUE
425 ELSE
426 work( ptrd ) = 1
427 END IF
428
429 IF( itype.EQ.3 ) THEN
430
431
432
433 CALL pdlaset(
'All', m, n, zero, zero, work( ptra ),
434 $ ia, ja, desca )
435
436 DO 40 i = 1, SIZE
437 CALL pdelset( work( ptra ), i, i, desca,
438 $ work( ptrd+i-1 ) )
439 40 CONTINUE
440
441 ELSE IF( itype.EQ.4 ) THEN
442
443
444
445 CALL pdlagge( m, n, work( ptrd ), work( ptra ), ia, ja,
446 $ desca, iseed, SIZE, work( ptrwork ),
447 $ llwork, info )
448
449 ELSE IF( itype.EQ.5 ) THEN
450
451
452
453 CALL dscal( SIZE, rtovfl, work( ptrd ), 1 )
454
455 CALL pdlagge( m, n, work( ptrd ), work( ptra ), ia, ja,
456 $ desca, iseed, SIZE, work( ptrwork ),
457 $ llwork, info )
458
459 ELSE IF( itype.EQ.6 ) THEN
460
461
462
463 CALL dscal( SIZE, rtunfl, work( ptrd ), 1 )
464 CALL pdlagge( m, n, work( ptrd ), work( ptra ), ia, ja,
465 $ desca, iseed, SIZE, work( ptrwork ),
466 $ llwork, info )
467
468 END IF
469
470 END IF
471
472
473
474
475 DO 80 jobtype = 1, 4
476
477 IF( jobtype.EQ.1 ) THEN
478 jobu = 'V'
479 jobvt = 'V'
480 ptrvt = ptru + ldu*sizeq
481 ptruc = ptrvt + ldvt*nq
482 ptrwork = ptruc + ldu*sizeq
483 llwork = lwork - ptrwork + 1
484 ELSE IF( jobtype.EQ.2 ) THEN
485 jobu = 'V'
486 jobvt = 'N'
487 ELSE IF( jobtype.EQ.3 ) THEN
488 jobu = 'N'
489 jobvt = 'V'
490 ptrvtc = ptruc
491 ptrwork = ptrvtc + ldvt*nq
492 llwork = lwork - ptrwork + 1
493 ELSE IF( jobtype.EQ.4 ) THEN
494 jobu = 'N'
495 jobvt = 'N'
496 ptrwork = ptruc
497 llwork = lwork - ptrwork + 1
498 END IF
499
500
501
502 CALL pdlacpy(
'A', m, n, work( ptra ), ia, ja, desca,
503 $ work( ptrac ), ia, ja, desca )
504
505
506
507
508 IF( jobtype.EQ.1 ) THEN
509
510
511
513 CALL blacs_barrier( context, 'All' )
515
516 CALL pdgesvd( jobu, jobvt, m, n, work( ptrac ), ia, ja,
517 $ desca, work( ptrs ), work( ptru ), iu, ju,
518 $ descu, work( ptrvt ), ivt, jvt, descvt,
519 $ work( ptrwork ), wpdgesvd, info )
520
522 CALL slcombine( context,
'All',
'>',
'W', 1, 1, wtime )
523 CALL slcombine( context,
'All',
'>',
'C', 1, 1, ctime )
524
525
526
527
528 itmp( 1 ) = info
529 itmp( 2 ) = info
530
531 CALL igamn2d( desca( ctxt_ ), 'a', ' ', 1, 1, itmp, 1, 1,
532 $ 1, -1, -1, 0 )
533 CALL igamx2d( desca( ctxt_ ), 'a', ' ', 1, 1, itmp( 2 ),
534 $ 1, 1, 1, -1, -1, 0 )
535
536 IF( itmp( 1 ).NE.itmp( 2 ) ) THEN
537 IF( myrow.EQ.0 .AND. mycol.EQ.0 ) THEN
538 WRITE( nout, fmt = * )
539 $ 'Different processes return different INFO'
540 GO TO 120
541 END IF
542 END IF
543
544
545
546
547
548 IF( info.EQ.( size+1 ) ) THEN
549 hetero = 'P'
550 sethet = 1
551 END IF
552
553
554
555 IF( info.EQ.zero ) THEN
556
557 DO 50 i = 1, SIZE
558 work( i+ptrwork ) = work( i+ptrs-1 )
559 work( i+size+ptrwork ) = work( i+ptrs-1 )
560 50 CONTINUE
561
562 CALL dgamn2d( desca( ctxt_ ), 'a', ' ', SIZE, 1,
563 $ work( 1+ptrwork ), SIZE, 1, 1, -1, -1,
564 $ 0 )
565 CALL dgamx2d( desca( ctxt_ ), 'a', ' ', SIZE, 1,
566 $ work( 1+size+ptrwork ), SIZE, 1, 1, -1,
567 $ -1, 0 )
568
569 DO 60 i = 1, SIZE
570 IF( abs( work( i+ptrwork )-work( size+i+
571 $ ptrwork ) ).GT.zero ) THEN
572 WRITE( nout, fmt = * )'I= ', i, ' MIN=',
573 $ work( i+ptrwork ), ' MAX=',
574 $ work( size+i+ptrwork )
575 hetero = 'T'
576 sethet = 1
577 GO TO 70
578 END IF
579
580 60 CONTINUE
581 70 CONTINUE
582
583 END IF
584
585 IF( sethet.NE.1 )
586 $ hetero = 'N'
587
588
589
590 CALL pdlacpy(
'A', m, n, work( ptra ), ia, ja, desca,
591 $ work( ptrac ), ia, ja, desca )
592
593
594
595
596
597 CALL pdlacpy(
'A', m,
SIZE, work( ptru ), iu, ju, descu,
598 $ work( ptruc ), iu, ju, descu )
599
600
601
602 CALL pdsvdchk( m, n, work( ptrac ), ia, ja, desca,
603 $ work( ptruc ), iu, ju, descu,
604 $ work( ptrvt ), ivt, jvt, descvt,
605 $ work( ptrs ), thresh, work( ptrwork ),
606 $ llwork, result, chk, mtm )
607
608 ELSE
609
610
611
612 CALL pdgesvd( jobu, jobvt, m, n, work( ptrac ), ia, ja,
613 $ desca, work( ptrsc ), work( ptruc ), iu,
614 $ ju, descu, work( ptrvtc ), ivt, jvt,
615 $ descvt, work( ptrwork ), wpdgesvd, info )
616
617 CALL pdsvdcmp( m, n, jobtype, work( ptrs ),
618 $ work( ptrsc ), work( ptru ),
619 $ work( ptruc ), iu, ju, descu,
620 $ work( ptrvt ), work( ptrvtc ), ivt, jvt,
621 $ descvt, thresh, result, delta,
622 $ work( ptrwork ), llwork )
623
624 END IF
625
626 80 CONTINUE
627
628 IF( myrow.EQ.0 .AND. mycol.EQ.0 ) THEN
629 DO 90 i = 1, 9
630 IF( result( i ).EQ.1 ) THEN
631 pass = 1
632 WRITE( nout, fmt = * )'Test I = ', i, 'has failed'
633 WRITE( nout, fmt = * )' '
634 END IF
635 90 CONTINUE
636 IF( pass.EQ.0 ) THEN
637 WRITE( nout, fmt = 9999 )'Passed', wtime( 1 ),
638 $ ctime( 1 ), m, n, nprow, npcol, nb, itype, chk, mtm,
639 $ delta, hetero
640 END IF
641 END IF
642 100 CONTINUE
643 CALL blacs_gridexit( context )
644 110 CONTINUE
645
646 9999 FORMAT( a6, 2e10.3, 2i6, 2i4, i5, i6, 3f6.2, 4x, a1 )
647 120 CONTINUE
648
649
650
651 RETURN
subroutine descinit(desc, m, n, mb, nb, irsrc, icsrc, ictxt, lld, info)
integer function numroc(n, nb, iproc, isrcproc, nprocs)
subroutine pdlaset(uplo, m, n, alpha, beta, a, ia, ja, desca)
double precision function pdlamch(ictxt, cmach)
subroutine pdelset(a, ia, ja, desca, alpha)
subroutine pdgesvd(jobu, jobvt, m, n, a, ia, ja, desca, s, u, iu, ju, descu, vt, ivt, jvt, descvt, work, lwork, info)
subroutine pdlacpy(uplo, m, n, a, ia, ja, desca, b, ib, jb, descb)
subroutine pdlagge(m, n, d, a, ia, ja, desca, iseed, order, work, lwork, info)
subroutine pdsvdchk(m, n, a, ia, ja, desca, u, iu, ju, descu, vt, ivt, jvt, descvt, s, thresh, work, lwork, result, chk, mtm)
subroutine pdsvdcmp(m, n, jobtype, s, sc, u, uc, iu, ju, descu, vt, vtc, ivt, jvt, descvt, thresh, result, delta, work, lwork)
subroutine pxerbla(ictxt, srname, info)
subroutine slcombine(ictxt, scope, op, timetype, n, ibeg, times)