3
4
5
6
7
8
9
10 CHARACTER JOBZ, UPLO
11 INTEGER IA, INFO, IZ, JA, JZ, LRWORK, LWORK, N
12
13
14 INTEGER DESCA( * ), DESCZ( * )
15 DOUBLE PRECISION RWORK( * ), W( * )
16 COMPLEX*16 A( * ), WORK( * ), Z( * )
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
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238 INTEGER BLOCK_CYCLIC_2D, DLEN_, DTYPE_, CTXT_, M_, N_,
239 $ MB_, NB_, RSRC_, CSRC_, LLD_
240 parameter( block_cyclic_2d = 1, dlen_ = 9, dtype_ = 1,
241 $ ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
242 $ rsrc_ = 7, csrc_ = 8, lld_ = 9 )
243 DOUBLE PRECISION ZERO, ONE
244 parameter( zero = 0.0d+0, one = 1.0d+0 )
245 COMPLEX*16 CZERO, CONE
246 parameter( czero = ( 0.0d+0, 0.0d+0 ),
247 $ cone = ( 1.0d+0, 0.0d+0 ) )
248 INTEGER ITHVAL
249 parameter( ithval = 10 )
250
251
252 LOGICAL LOWER, WANTZ
253 INTEGER CONTEXTC, CSRC_A, I, IACOL, IAROW, ICOFFA,
254 $ IINFO, INDD, INDE, INDRD, INDRE, INDRWORK,
255 $ INDTAU, INDWORK, INDWORK2, IROFFA, IROFFZ,
256 $ ISCALE, IZROW, J, K, LDC, LLRWORK, LLWORK,
257 $ LRMIN, LRWMIN, LWMIN, MB_A, MB_Z, MYCOL,
258 $ MYPCOLC, MYPROWC, MYROW, NB, NB_A, NB_Z, NP0,
259 $ NPCOL, NPCOLC, NPROCS, NPROW, NPROWC, NQ0, NRC,
260 $ RSIZEZSTEQR2, RSRC_A, RSRC_Z, SIZEPZHETRD,
261 $ SIZEPZUNMTR, SIZEZSTEQR2
262 DOUBLE PRECISION ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
263 $ SMLNUM
264
265
266 INTEGER DESCQR( 10 ), IDUM1( 3 ), IDUM2( 3 )
267
268
269 LOGICAL LSAME
270 INTEGER INDXG2P, NUMROC, SL_GRIDRESHAPE
271 DOUBLE PRECISION PDLAMCH, PZLANHE
274
275
276 EXTERNAL blacs_gridexit, blacs_gridinfo,
chk1mat, dcopy,
280
281
282 INTRINSIC abs, dble, dcmplx, ichar, int,
max,
min, mod,
283 $ sqrt
284
285
286
287 IF( block_cyclic_2d*csrc_*ctxt_*dlen_*dtype_*lld_*mb_*m_*nb_*n_*
288 $ rsrc_.LT.0 )RETURN
289
290
291
292 IF( n.EQ.0 )
293 $ RETURN
294
295
296
297 CALL blacs_gridinfo( desca( ctxt_ ), nprow, npcol, myrow, mycol )
298 info = 0
299
300
301
302 indtau = 1
303 indd = 1
304 inde = 1
305 indwork = 1
306 indwork2 = 1
307
308 indre = 1
309 indrd = 1
310 indrwork = 1
311
312 wantz =
lsame( jobz,
'V' )
313 IF( nprow.EQ.-1 ) THEN
314 info = -( 700+ctxt_ )
315 ELSE IF( wantz ) THEN
316 IF( desca( ctxt_ ).NE.descz( ctxt_ ) ) THEN
317 info = -( 1200+ctxt_ )
318 END IF
319 END IF
320 IF( info.EQ.0 ) THEN
321 CALL chk1mat( n, 3, n, 3, ia, ja, desca, 7, info )
322 IF( wantz )
323 $
CALL chk1mat( n, 3, n, 3, iz, jz, descz, 12, info )
324
325 IF( info.EQ.0 ) THEN
326
327
328
329 safmin =
pdlamch( desca( ctxt_ ),
'Safe minimum' )
330 eps =
pdlamch( desca( ctxt_ ),
'Precision' )
331 smlnum = safmin / eps
332 bignum = one / smlnum
333 rmin = sqrt( smlnum )
334 rmax =
min( sqrt( bignum ), one / sqrt( sqrt( safmin ) ) )
335
336 nprocs = nprow*npcol
337 nb_a = desca( nb_ )
338 mb_a = desca( mb_ )
339 nb = nb_a
340 lower =
lsame( uplo,
'L' )
341
342 rsrc_a = desca( rsrc_ )
343 csrc_a = desca( csrc_ )
344 iroffa = mod( ia-1, mb_a )
345 icoffa = mod( ja-1, nb_a )
346 iarow =
indxg2p( 1, nb_a, myrow, rsrc_a, nprow )
347 iacol =
indxg2p( 1, mb_a, mycol, csrc_a, npcol )
348 np0 =
numroc( n+iroffa, nb, myrow, iarow, nprow )
349 nq0 =
numroc( n+icoffa, nb, mycol, iacol, npcol )
350 IF( wantz ) THEN
351 nb_z = descz( nb_ )
352 mb_z = descz( mb_ )
353 rsrc_z = descz( rsrc_ )
354 iroffz = mod( iz-1, mb_a )
355 izrow =
indxg2p( 1, nb_a, myrow, rsrc_z, nprow )
356 ELSE
357 iroffz = 0
358 izrow = 0
359 END IF
360
361
362
363 CALL pzhetrd( uplo, n, a, ia, ja, desca, rwork( indd ),
364 $ rwork( inde ), work( indtau ),
365 $ work( indwork ), -1, iinfo )
366 sizepzhetrd = int( abs( work( 1 ) ) )
367
368
369
370 IF( wantz ) THEN
371 CALL pzunmtr(
'L', uplo,
'N', n, n, a, ia, ja, desca,
372 $ work( indtau ), z, iz, jz, descz,
373 $ work( indwork ), -1, iinfo )
374 sizepzunmtr = int( abs( work( 1 ) ) )
375 ELSE
376 sizepzunmtr = 0
377 END IF
378
379
380
381 IF( wantz ) THEN
382 rsizezsteqr2 =
max( 1, 2*n-2 )
383 ELSE
384 rsizezsteqr2 = 0
385 END IF
386
387
388
389
390
391
392
393 ldc = 0
394 IF( wantz ) THEN
396 $ nprocs, 1 )
397 CALL blacs_gridinfo( contextc, nprowc, npcolc, myprowc,
398 $ mypcolc )
399 nrc =
numroc( n, nb_a, myprowc, 0, nprocs )
401 CALL descinit( descqr, n, n, nb, nb, 0, 0, contextc, ldc,
402 $ info )
403 END IF
404
405
406
407 IF( wantz ) THEN
408 sizezsteqr2 = n*ldc
409 ELSE
410 sizezsteqr2 = 0
411 END IF
412
413
414
415 indtau = 1
416 indd = indtau + n
417 inde = indd + n
418 indwork = inde + n
419 indwork2 = indwork + n*ldc
420 llwork = lwork - indwork + 1
421
422
423
424 indre = 1
425 indrd = indre + n
426 indrwork = indrd + n
427 llrwork = lrwork - indrwork + 1
428
429
430
431 lrwmin = 2*n + rsizezsteqr2
432 lwmin = 3*n +
max( sizepzhetrd, sizepzunmtr, sizezsteqr2 )
433
434 END IF
435 IF( info.EQ.0 ) THEN
436 IF( .NOT.( wantz .OR.
lsame( jobz,
'N' ) ) )
THEN
437 info = -1
438 ELSE IF( .NOT.( lower .OR.
lsame( uplo,
'U' ) ) )
THEN
439 info = -2
440 ELSE IF( lwork.LT.lwmin .AND. lwork.NE.-1 ) THEN
441 info = -14
442 ELSE IF( lrwork.LT.lrwmin .AND. lrwork.NE.-1 ) THEN
443 info = -16
444 ELSE IF( iroffa.NE.0 ) THEN
445 info = -5
446 ELSE IF( desca( mb_ ).NE.desca( nb_ ) ) THEN
447 info = -( 700+nb_ )
448 END IF
449 IF( wantz ) THEN
450 IF( iroffa.NE.iroffz ) THEN
451 info = -10
452 ELSE IF( iarow.NE.izrow ) THEN
453 info = -10
454 ELSE IF( desca( m_ ).NE.descz( m_ ) ) THEN
455 info = -( 1200+m_ )
456 ELSE IF( desca( n_ ).NE.descz( n_ ) ) THEN
457 info = -( 1200+n_ )
458 ELSE IF( desca( mb_ ).NE.descz( mb_ ) ) THEN
459 info = -( 1200+mb_ )
460 ELSE IF( desca( nb_ ).NE.descz( nb_ ) ) THEN
461 info = -( 1200+nb_ )
462 ELSE IF( desca( rsrc_ ).NE.descz( rsrc_ ) ) THEN
463 info = -( 1200+rsrc_ )
464 ELSE IF( desca( ctxt_ ).NE.descz( ctxt_ ) ) THEN
465 info = -( 1200+ctxt_ )
466 END IF
467 END IF
468 END IF
469 IF( wantz ) THEN
470 idum1( 1 ) = ichar( 'V' )
471 ELSE
472 idum1( 1 ) = ichar( 'N' )
473 END IF
474 idum2( 1 ) = 1
475 IF( lower ) THEN
476 idum1( 2 ) = ichar( 'L' )
477 ELSE
478 idum1( 2 ) = ichar( 'U' )
479 END IF
480 idum2( 2 ) = 2
481 IF( lwork.EQ.-1 ) THEN
482 idum1( 3 ) = -1
483 ELSE
484 idum1( 3 ) = 1
485 END IF
486 idum2( 3 ) = 3
487 IF( wantz ) THEN
488 CALL pchk2mat( n, 3, n, 3, ia, ja, desca, 7, n, 3, n, 3, iz,
489 $ jz, descz, 12, 3, idum1, idum2, info )
490 ELSE
491 CALL pchk1mat( n, 3, n, 3, ia, ja, desca, 7, 3, idum1,
492 $ idum2, info )
493 END IF
494 work( 1 ) = dcmplx( lwmin )
495 rwork( 1 ) = dble( lrwmin )
496 END IF
497
498 IF( info.NE.0 ) THEN
499 CALL pxerbla( desca( ctxt_ ),
'PZHEEV', -info )
500 IF( wantz )
501 $ CALL blacs_gridexit( contextc )
502 RETURN
503 ELSE IF( lwork.EQ.-1 .OR. lrwork.EQ.-1 ) THEN
504 IF( wantz )
505 $ CALL blacs_gridexit( contextc )
506 RETURN
507 END IF
508
509
510
511 iscale = 0
512
513 anrm =
pzlanhe(
'M', uplo, n, a, ia, ja, desca,
514 $ rwork( indrwork ) )
515
516
517 IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
518 iscale = 1
519 sigma = rmin / anrm
520 ELSE IF( anrm.GT.rmax ) THEN
521 iscale = 1
522 sigma = rmax / anrm
523 END IF
524
525 IF( iscale.EQ.1 ) THEN
526 CALL pzlascl( uplo, one, sigma, n, n, a, ia, ja, desca, iinfo )
527 END IF
528
529
530
531 CALL pzhetrd( uplo, n, a, ia, ja, desca, rwork( indrd ),
532 $ rwork( indre ), work( indtau ), work( indwork ),
533 $ llwork, iinfo )
534
535
536
537 DO 10 i = 1, n
538 CALL pzelget(
'A',
' ', work( indd+i-1 ), a, i+ia-1, i+ja-1,
539 $ desca )
540 rwork( indrd+i-1 ) = dble( work( indd+i-1 ) )
541 10 CONTINUE
542 IF(
lsame( uplo,
'U' ) )
THEN
543 DO 20 i = 1, n - 1
544 CALL pzelget(
'A',
' ', work( inde+i-1 ), a, i+ia-1, i+ja,
545 $ desca )
546 rwork( indre+i-1 ) = dble( work( inde+i-1 ) )
547 20 CONTINUE
548 ELSE
549 DO 30 i = 1, n - 1
550 CALL pzelget(
'A',
' ', work( inde+i-1 ), a, i+ia, i+ja-1,
551 $ desca )
552 rwork( indre+i-1 ) = dble( work( inde+i-1 ) )
553 30 CONTINUE
554 END IF
555
556 IF( wantz ) THEN
557
558 CALL pzlaset(
'Full', n, n, czero, cone, work( indwork ), 1, 1,
559 $ descqr )
560
561
562
563
564
565 CALL zsteqr2(
'I', n, rwork( indrd ), rwork( indre ),
566 $ work( indwork ), ldc, nrc, rwork( indrwork ),
567 $ info )
568
569 CALL pzgemr2d( n, n, work( indwork ), 1, 1, descqr, z, ia, ja,
570 $ descz, contextc )
571
572 CALL pzunmtr(
'L', uplo,
'N', n, n, a, ia, ja, desca,
573 $ work( indtau ), z, iz, jz, descz,
574 $ work( indwork ), llwork, iinfo )
575
576 ELSE
577
578 CALL zsteqr2(
'N', n, rwork( indrd ), rwork( indre ),
579 $ work( indwork ), 1, 1, rwork( indrwork ), info )
580 END IF
581
582
583
584 CALL dcopy( n, rwork( indd ), 1, w, 1 )
585
586
587
588 IF( iscale.EQ.1 ) THEN
589 CALL dscal( n, one / sigma, w, 1 )
590 END IF
591
592 work( 1 ) = dble( lwmin )
593
594
595
596 IF( wantz ) THEN
597 CALL blacs_gridexit( contextc )
598 END IF
599
600
601
602
603 IF( n.LE.ithval ) THEN
604 j = n
605 k = 1
606 ELSE
607 j = n / ithval
608 k = ithval
609 END IF
610
611 lrmin = int( rwork( 1 ) )
612 indtau = 0
613 inde = indtau + j
614 DO 40 i = 1, j
615 rwork( i+indtau ) = w( ( i-1 )*k+1 )
616 rwork( i+inde ) = w( ( i-1 )*k+1 )
617 40 CONTINUE
618
619 CALL dgamn2d( desca( ctxt_ ), 'All', ' ', j, 1, rwork( 1+indtau ),
620 $ j, 1, 1, -1, -1, 0 )
621 CALL dgamx2d( desca( ctxt_ ), 'All', ' ', j, 1, rwork( 1+inde ),
622 $ j, 1, 1, -1, -1, 0 )
623
624 DO 50 i = 1, j
625 IF( info.EQ.0 .AND. ( rwork( i+indtau )-rwork( i+inde ).NE.
626 $ zero ) ) THEN
627 info = n + 1
628 END IF
629 50 CONTINUE
630 rwork( 1 ) = lrmin
631
632 RETURN
633
634
635
Int sl_gridreshape(Int *ctxt, Int *pstart, Int *row_major_in, Int *row_major_out, Int *P, Int *Q)
subroutine chk1mat(ma, mapos0, na, napos0, ia, ja, desca, descapos0, info)
subroutine descinit(desc, m, n, mb, nb, irsrc, icsrc, ictxt, lld, info)
integer function indxg2p(indxglob, nb, iproc, isrcproc, nprocs)
integer function numroc(n, nb, iproc, isrcproc, nprocs)
subroutine pchk1mat(ma, mapos0, na, napos0, ia, ja, desca, descapos0, nextra, ex, expos, info)
subroutine pchk2mat(ma, mapos0, na, napos0, ia, ja, desca, descapos0, mb, mbpos0, nb, nbpos0, ib, jb, descb, descbpos0, nextra, ex, expos, info)
double precision function pdlamch(ictxt, cmach)
subroutine pxerbla(ictxt, srname, info)
subroutine pzlaset(uplo, m, n, alpha, beta, a, ia, ja, desca)
subroutine pzelget(scope, top, alpha, a, ia, ja, desca)
subroutine pzhetrd(uplo, n, a, ia, ja, desca, d, e, tau, work, lwork, info)
double precision function pzlanhe(norm, uplo, n, a, ia, ja, desca, work)
subroutine pzlascl(type, cfrom, cto, m, n, a, ia, ja, desca, info)
subroutine pzunmtr(side, uplo, trans, m, n, a, ia, ja, desca, tau, c, ic, jc, descc, work, lwork, info)
subroutine zsteqr2(compz, n, d, e, z, ldz, nr, work, info)