3
4
5
6
7
8
9
10 CHARACTER UPLO
11 INTEGER IA, INFO, JA, LRWORK, LWORK, N
12
13
14 INTEGER DESCA( * )
15 REAL D( * ), E( * ), RWORK( * )
16 COMPLEX A( * ), TAU( * ), WORK( * )
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
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266 INTEGER BLOCK_CYCLIC_2D, DLEN_, DTYPE_, CTXT_, M_, N_,
267 $ MB_, NB_, RSRC_, CSRC_, LLD_
268 parameter( block_cyclic_2d = 1, dlen_ = 9, dtype_ = 1,
269 $ ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
270 $ rsrc_ = 7, csrc_ = 8, lld_ = 9 )
271 REAL ONE
272 parameter( one = 1.0e+0 )
273 COMPLEX CONE
274 parameter( cone = ( 1.0e+0, 0.0e+0 ) )
275
276
277 LOGICAL LQUERY, UPPER
278 CHARACTER COLCTOP, ROWCTOP
279 INTEGER ANB, CTXTB, I, IACOL, IAROW, ICOFFA, ICTXT,
280 $ IINFO, INDB, INDRD, INDRE, INDTAU, INDW, IPW,
281 $ IROFFA, J, JB, JX, K, KK, LLRWORK, LLWORK,
282 $ LRWMIN, LWMIN, MINSZ, MYCOL, MYCOLB, MYROW,
283 $ MYROWB, NB, NP, NPCOL, NPCOLB, NPROW, NPROWB,
284 $ NPS, NQ, ONEPMIN, ONEPRMIN, SQNPC, TTLRWMIN,
285 $ TTLWMIN
286
287
288 INTEGER DESCB( DLEN_ ), DESCW( DLEN_ ), IDUM1( 3 ),
289 $ IDUM2( 3 )
290
291
292 EXTERNAL blacs_get, blacs_gridexit, blacs_gridinfo,
297
298
299 LOGICAL LSAME
300 INTEGER INDXG2L, INDXG2P, NUMROC, PJLAENV
302
303
304 INTRINSIC cmplx, ichar, int,
max,
min, mod, real, sqrt
305
306
307
308
309 IF( block_cyclic_2d*csrc_*ctxt_*dlen_*dtype_*lld_*mb_*m_*nb_*n_*
310 $ rsrc_.LT.0 )RETURN
311
312
313 ictxt = desca( ctxt_ )
314 CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
315
316
317
318 info = 0
319 IF( nprow.EQ.-1 ) THEN
320 info = -( 600+ctxt_ )
321 ELSE
322 CALL chk1mat( n, 2, n, 2, ia, ja, desca, 6, info )
323 upper =
lsame( uplo,
'U' )
324 IF( info.EQ.0 ) THEN
325 nb = desca( nb_ )
326 iroffa = mod( ia-1, desca( mb_ ) )
327 icoffa = mod( ja-1, desca( nb_ ) )
328 iarow =
indxg2p( ia, nb, myrow, desca( rsrc_ ), nprow )
329 iacol =
indxg2p( ja, nb, mycol, desca( csrc_ ), npcol )
330 np =
numroc( n, nb, myrow, iarow, nprow )
331 nq =
max( 1,
numroc( n+ja-1, nb, mycol, desca( csrc_ ),
332 $ npcol ) )
333 lwmin =
max( ( np+1 )*nb, 3*nb )
334 anb =
pjlaenv( ictxt, 3,
'PCHETTRD',
'L', 0, 0, 0, 0 )
335 minsz =
pjlaenv( ictxt, 5,
'PCHETTRD',
'L', 0, 0, 0, 0 )
336 sqnpc = int( sqrt( real( nprow*npcol ) ) )
337 nps =
max(
numroc( n, 1, 0, 0, sqnpc ), 2*anb )
338 ttlwmin = 2*( anb+1 )*( 4*nps+2 ) + ( nps+2 )*nps
339 lrwmin = 1
340 ttlrwmin = 2*nps
341
342 work( 1 ) =
cmplx( real( ttlwmin ) )
343 rwork( 1 ) = real( ttlrwmin )
344 lquery = ( lwork.EQ.-1 .OR. lrwork.EQ.-1 )
345 IF( .NOT.upper .AND. .NOT.
lsame( uplo,
'L' ) )
THEN
346 info = -1
347
348
349
350
351 ELSE IF( iroffa.NE.icoffa .OR. icoffa.NE.0 ) THEN
352 info = -5
353 ELSE IF( desca( mb_ ).NE.desca( nb_ ) ) THEN
354 info = -( 600+nb_ )
355 ELSE IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
356 info = -11
357 ELSE IF( lrwork.LT.lrwmin .AND. .NOT.lquery ) THEN
358 info = -13
359 END IF
360 END IF
361 IF( upper ) THEN
362 idum1( 1 ) = ichar( 'U' )
363 ELSE
364 idum1( 1 ) = ichar( 'L' )
365 END IF
366 idum2( 1 ) = 1
367 IF( lwork.EQ.-1 ) THEN
368 idum1( 2 ) = -1
369 ELSE
370 idum1( 2 ) = 1
371 END IF
372 idum2( 2 ) = 11
373 IF( lrwork.EQ.-1 ) THEN
374 idum1( 3 ) = -1
375 ELSE
376 idum1( 3 ) = 1
377 END IF
378 idum2( 3 ) = 13
379 CALL pchk1mat( n, 2, n, 2, ia, ja, desca, 6, 3, idum1, idum2,
380 $ info )
381 END IF
382
383 IF( info.NE.0 ) THEN
384 CALL pxerbla( ictxt,
'PCHENTRD', -info )
385 RETURN
386 ELSE IF( lquery ) THEN
387 RETURN
388 END IF
389
390
391
392 IF( n.EQ.0 )
393 $ RETURN
394
395
396 onepmin = n*n + 3*n + 1
397 llwork = lwork
398 CALL igamn2d( ictxt, 'A', ' ', 1, 1, llwork, 1, 1, -1, -1, -1,
399 $ -1 )
400
401 oneprmin = 2*n
402 llrwork = lrwork
403 CALL igamn2d( ictxt, 'A', ' ', 1, 1, llrwork, 1, 1, -1, -1, -1,
404 $ -1 )
405
406
407
408
409
410 nprowb = 0
411 IF( ( n.LT.minsz .OR. sqnpc.EQ.1 ) .AND. llwork.GE.onepmin .AND.
412 $ llrwork.GE.oneprmin .AND. .NOT.upper ) THEN
413 nprowb = 1
414 nps = n
415 ELSE
416 IF( llwork.GE.ttlwmin .AND. llrwork.GE.ttlrwmin .AND. .NOT.
417 $ upper ) THEN
418 nprowb = sqnpc
419 END IF
420 END IF
421
422 IF( nprowb.GE.1 ) THEN
423 npcolb = nprowb
424 sqnpc = nprowb
425 indb = 1
426 indrd = 1
427 indre = indrd + nps
428 indtau = indb + nps*nps
429 indw = indtau + nps
430 llwork = llwork - indw + 1
431
432 CALL blacs_get( ictxt, 10, ctxtb )
433 CALL blacs_gridinit( ctxtb, 'Row major', sqnpc, sqnpc )
434 CALL blacs_gridinfo( ctxtb, nprowb, npcolb, myrowb, mycolb )
435 CALL descset( descb, n, n, 1, 1, 0, 0, ctxtb, nps )
436
437 CALL pctrmr2d( uplo, 'N', n, n, a, ia, ja, desca, work( indb ),
438 $ 1, 1, descb, ictxt )
439
440
441
442
443 IF( nprowb.GT.0 ) THEN
444
445 IF( nprowb.EQ.1 ) THEN
446 CALL chetrd( uplo, n, work( indb ), nps, rwork( indrd ),
447 $ rwork( indre ), work( indtau ),
448 $ work( indw ), llwork, info )
449 ELSE
450
451 CALL pchettrd(
'L', n, work( indb ), 1, 1, descb,
452 $ rwork( indrd ), rwork( indre ),
453 $ work( indtau ), work( indw ), llwork,
454 $ info )
455
456 END IF
457 END IF
458
459
460
461
462 CALL pslamr1d( n-1, rwork( indre ), 1, 1, descb, e, 1, ja,
463 $ desca )
464
465 CALL pslamr1d( n, rwork( indrd ), 1, 1, descb, d, 1, ja,
466 $ desca )
467
468 CALL pclamr1d( n, work( indtau ), 1, 1, descb, tau, 1, ja,
469 $ desca )
470
471 CALL pctrmr2d( uplo, 'N', n, n, work( indb ), 1, 1, descb, a,
472 $ ia, ja, desca, ictxt )
473
474 IF( myrowb.GE.0 )
475 $ CALL blacs_gridexit( ctxtb )
476
477 ELSE
478
479 CALL pb_topget( ictxt, 'Combine', 'Columnwise', colctop )
480 CALL pb_topget( ictxt, 'Combine', 'Rowwise', rowctop )
481 CALL pb_topset( ictxt, 'Combine', 'Columnwise', '1-tree' )
482 CALL pb_topset( ictxt, 'Combine', 'Rowwise', '1-tree' )
483
484 ipw = np*nb + 1
485
486 IF( upper ) THEN
487
488
489
490 kk = mod( ja+n-1, nb )
491 IF( kk.EQ.0 )
492 $ kk = nb
493 CALL descset( descw, n, nb, nb, nb, iarow,
494 $
indxg2p( ja+n-kk, nb, mycol, desca( csrc_ ),
495 $ npcol ), ictxt,
max( 1, np ) )
496
497 DO 10 k = n - kk + 1, nb + 1, -nb
498 jb =
min( n-k+1, nb )
499 i = ia + k - 1
500 j = ja + k - 1
501
502
503
504
505
506 CALL pclatrd( uplo, k+jb-1, jb, a, ia, ja, desca, d, e,
507 $ tau, work, 1, 1, descw, work( ipw ) )
508
509
510
511
512
513 CALL pcher2k( uplo, 'No transpose', k-1, jb, -cone, a,
514 $ ia, j, desca, work, 1, 1, descw, one, a,
515 $ ia, ja, desca )
516
517
518
519 jx =
min(
indxg2l( j, nb, 0, iacol, npcol ), nq )
521
522 descw( csrc_ ) = mod( descw( csrc_ )+npcol-1, npcol )
523
524 10 CONTINUE
525
526
527
528 CALL pchetd2( uplo,
min( n, nb ), a, ia, ja, desca, d, e,
529 $ tau, work, lwork, iinfo )
530
531 ELSE
532
533
534
535 kk = mod( ja+n-1, nb )
536 IF( kk.EQ.0 )
537 $ kk = nb
538 CALL descset( descw, n, nb, nb, nb, iarow, iacol, ictxt,
540
541 DO 20 k = 1, n - nb, nb
542 i = ia + k - 1
543 j = ja + k - 1
544
545
546
547
548
549 CALL pclatrd( uplo, n-k+1, nb, a, i, j, desca, d, e, tau,
550 $ work, k, 1, descw, work( ipw ) )
551
552
553
554
555
556 CALL pcher2k( uplo, 'No transpose', n-k-nb+1, nb, -cone,
557 $ a, i+nb, j, desca, work, k+nb, 1, descw,
558 $ one, a, i+nb, j+nb, desca )
559
560
561
562 jx =
min(
indxg2l( j+nb-1, nb, 0, iacol, npcol ), nq )
563 CALL pcelset( a, i+nb, j+nb-1, desca,
cmplx( e( jx ) ) )
564
565 descw( csrc_ ) = mod( descw( csrc_ )+1, npcol )
566
567 20 CONTINUE
568
569
570
571 CALL pchetd2( uplo, kk, a, ia+k-1, ja+k-1, desca, d, e, tau,
572 $ work, lwork, iinfo )
573 END IF
574
575 CALL pb_topset( ictxt, 'Combine', 'Columnwise', colctop )
576 CALL pb_topset( ictxt, 'Combine', 'Rowwise', rowctop )
577
578 END IF
579
580 work( 1 ) =
cmplx( real( ttlwmin ) )
581 rwork( 1 ) = real( ttlrwmin )
582
583 RETURN
584
585
586
subroutine chk1mat(ma, mapos0, na, napos0, ia, ja, desca, descapos0, info)
subroutine descset(desc, m, n, mb, nb, irsrc, icsrc, ictxt, lld)
integer function indxg2l(indxglob, nb, iproc, isrcproc, nprocs)
integer function indxg2p(indxglob, nb, iproc, isrcproc, nprocs)
integer function numroc(n, nb, iproc, isrcproc, nprocs)
subroutine pcelset(a, ia, ja, desca, alpha)
subroutine pchetd2(uplo, n, a, ia, ja, desca, d, e, tau, work, lwork, info)
subroutine pchettrd(uplo, n, a, ia, ja, desca, d, e, tau, work, lwork, info)
subroutine pchk1mat(ma, mapos0, na, napos0, ia, ja, desca, descapos0, nextra, ex, expos, info)
subroutine pclamr1d(n, a, ia, ja, desca, b, ib, jb, descb)
subroutine pclatrd(uplo, n, nb, a, ia, ja, desca, d, e, tau, w, iw, jw, descw, work)
integer function pjlaenv(ictxt, ispec, name, opts, n1, n2, n3, n4)
subroutine pslamr1d(n, a, ia, ja, desca, b, ib, jb, descb)
subroutine pxerbla(ictxt, srname, info)