3
4
5
6
7
8
9
10 INTEGER LWORK, M, N, NB, NOUT, NPCOL, NPROW
11 REAL THRESH
12
13
14 INTEGER ISEED( 4 ), RESULT( 9 )
15 REAL 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 REAL ZERO, ONE
216 parameter( zero = 0.0e0, one = 1.0e0 )
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, PTRS, PTRSC, PTRU, PTRUC,
224 $ PTRVT, PTRVTC, PTRWORK, SETHET, SIZE, SIZEQ,
225 $ WPSGESVD, WPSLAGGE, WPSSVDCHK, WPSSVDCMP
226 REAL 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, sgamn2d, sgamx2d, slabad, sscal,
234 $ igamn2d, igamx2d, igebr2d, igebs2d,
pselset,
237
238
239 INTEGER NUMROC
240 REAL PSLAMCH
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 pslagge( m, n, work( ptrd ), work( ptra ), ia, ja, desca,
325 $ iseed, SIZE, work( ptrwork ), -1, dinfo )
326 wpslagge = int( work( ptrwork ) )
327
328 CALL psgesvd(
'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 wpsgesvd = int( work( ptrwork ) )
333
334 CALL pssvdchk( 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 wpssvdchk = int( work( ptrwork ) )
339
340 CALL pssvdcmp( 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 wpssvdcmp = int( work( ptrwork ) )
345
346
347
348 lwmin = 1 + 2*lda*nq + 3*SIZE +
349 $
max( wpslagge, ldu*sizeq+ldvt*nq+
max( ldu*sizeq,
350 $ ldvt*nq )+wpsgesvd+
max( wpssvdchk, wpssvdcmp ) )
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_ ),
'PSSVDTST', -info )
365 RETURN
366 END IF
367
369 unfl =
pslamch( context,
'Safe min' )
370 ovfl = one / unfl
371 CALL slabad( 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 pslaset(
'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 pslaset(
'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 pslaset(
'All', m, n, zero, zero, work( ptra ),
434 $ ia, ja, desca )
435
436 DO 40 i = 1, SIZE
437 CALL pselset( 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 pslagge( 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 sscal( SIZE, rtovfl, work( ptrd ), 1 )
454
455 CALL pslagge( 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 sscal( SIZE, rtunfl, work( ptrd ), 1 )
464 CALL pslagge( 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 pslacpy(
'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
512 CALL blacs_barrier( context, 'All' )
514
515 CALL psgesvd( jobu, jobvt, m, n, work( ptrac ), ia, ja,
516 $ desca, work( ptrs ), work( ptru ), iu, ju,
517 $ descu, work( ptrvt ), ivt, jvt, descvt,
518 $ work( ptrwork ), wpsgesvd, info )
519
521 CALL slcombine( context,
'All',
'>',
'W', 1, 1, wtime )
522 CALL slcombine( context,
'All',
'>',
'C', 1, 1, ctime )
523
524
525
526
527 itmp( 1 ) = info
528 itmp( 2 ) = info
529
530 CALL igamn2d( desca( ctxt_ ), 'a', ' ', 1, 1, itmp, 1, 1,
531 $ 1, -1, -1, 0 )
532 CALL igamx2d( desca( ctxt_ ), 'a', ' ', 1, 1, itmp( 2 ),
533 $ 1, 1, 1, -1, -1, 0 )
534
535 IF( itmp( 1 ).NE.itmp( 2 ) ) THEN
536 IF( myrow.EQ.0 .AND. mycol.EQ.0 ) THEN
537 WRITE( nout, fmt = * )
538 $ 'Different processes return different INFO'
539 GO TO 120
540 END IF
541 END IF
542
543
544
545
546
547 IF( info.EQ.( size+1 ) ) THEN
548 hetero = 'P'
549 sethet = 1
550 END IF
551
552
553
554 IF( info.EQ.zero ) THEN
555
556 DO 50 i = 1, SIZE
557 work( i+ptrwork ) = work( i+ptrs-1 )
558 work( i+size+ptrwork ) = work( i+ptrs-1 )
559 50 CONTINUE
560
561 CALL sgamn2d( desca( ctxt_ ), 'a', ' ', SIZE, 1,
562 $ work( 1+ptrwork ), SIZE, 1, 1, -1, -1,
563 $ 0 )
564 CALL sgamx2d( desca( ctxt_ ), 'a', ' ', SIZE, 1,
565 $ work( 1+size+ptrwork ), SIZE, 1, 1, -1,
566 $ -1, 0 )
567
568 DO 60 i = 1, SIZE
569 IF( abs( work( i+ptrwork )-work( size+i+
570 $ ptrwork ) ).GT.zero ) THEN
571 WRITE( nout, fmt = * )'I= ', i, ' MIN=',
572 $ work( i+ptrwork ), ' MAX=',
573 $ work( size+i+ptrwork )
574 hetero = 'T'
575 sethet = 1
576 GO TO 70
577 END IF
578
579 60 CONTINUE
580 70 CONTINUE
581
582 END IF
583
584 IF( sethet.NE.1 )
585 $ hetero = 'N'
586
587
588
589 CALL pslacpy(
'A', m, n, work( ptra ), ia, ja, desca,
590 $ work( ptrac ), ia, ja, desca )
591
592
593
594
595
596 CALL pslacpy(
'A', m,
SIZE, work( ptru ), iu, ju, descu,
597 $ work( ptruc ), iu, ju, descu )
598
599
600
601 CALL pssvdchk( m, n, work( ptrac ), ia, ja, desca,
602 $ work( ptruc ), iu, ju, descu,
603 $ work( ptrvt ), ivt, jvt, descvt,
604 $ work( ptrs ), thresh, work( ptrwork ),
605 $ llwork, result, chk, mtm )
606
607 ELSE
608
609
610
611 CALL psgesvd( jobu, jobvt, m, n, work( ptrac ), ia, ja,
612 $ desca, work( ptrsc ), work( ptruc ), iu,
613 $ ju, descu, work( ptrvtc ), ivt, jvt,
614 $ descvt, work( ptrwork ), wpsgesvd, info )
615
616 CALL pssvdcmp( m, n, jobtype, work( ptrs ),
617 $ work( ptrsc ), work( ptru ),
618 $ work( ptruc ), iu, ju, descu,
619 $ work( ptrvt ), work( ptrvtc ), ivt, jvt,
620 $ descvt, thresh, result, delta,
621 $ work( ptrwork ), llwork )
622
623 END IF
624
625 80 CONTINUE
626
627 IF( myrow.EQ.0 .AND. mycol.EQ.0 ) THEN
628 DO 90 i = 1, 9
629 IF( result( i ).EQ.1 ) THEN
630 pass = 1
631 WRITE( nout, fmt = * )'Test I = ', i, 'has failed'
632 WRITE( nout, fmt = * )' '
633 END IF
634 90 CONTINUE
635 IF( pass.EQ.0 ) THEN
636 WRITE( nout, fmt = 9999 )'Passed', wtime( 1 ),
637 $ ctime( 1 ), m, n, nprow, npcol, nb, itype, chk, mtm,
638 $ delta, hetero
639 END IF
640 END IF
641 100 CONTINUE
642 CALL blacs_gridexit( context )
643 110 CONTINUE
644
645 9999 FORMAT( a6, 2e10.3, 2i6, 2i4, i5, i6, 3f6.2, 4x, a1 )
646 120 CONTINUE
647
648
649
650 RETURN
subroutine descinit(desc, m, n, mb, nb, irsrc, icsrc, ictxt, lld, info)
integer function numroc(n, nb, iproc, isrcproc, nprocs)
real function pslamch(ictxt, cmach)
subroutine pslaset(uplo, m, n, alpha, beta, a, ia, ja, desca)
subroutine pselset(a, ia, ja, desca, alpha)
subroutine psgesvd(jobu, jobvt, m, n, a, ia, ja, desca, s, u, iu, ju, descu, vt, ivt, jvt, descvt, work, lwork, info)
subroutine pslacpy(uplo, m, n, a, ia, ja, desca, b, ib, jb, descb)
subroutine pslagge(m, n, d, a, ia, ja, desca, iseed, order, work, lwork, info)
subroutine pssvdchk(m, n, a, ia, ja, desca, u, iu, ju, descu, vt, ivt, jvt, descvt, s, thresh, work, lwork, result, chk, mtm)
subroutine pssvdcmp(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)