160
161
162
163
164
165 IMPLICIT NONE
166
167
168 CHARACTER UPLO
169 INTEGER N, LDA, LTB, LWORK, INFO
170
171
172 INTEGER IPIV( * ), IPIV2( * )
173 COMPLEX*16 A( LDA, * ), TB( * ), WORK( * )
174
175
176
177
178 COMPLEX*16 ZERO, ONE
179 parameter( zero = ( 0.0e+0, 0.0e+0 ),
180 $ one = ( 1.0e+0, 0.0e+0 ) )
181
182
183 LOGICAL UPPER, TQUERY, WQUERY
184 INTEGER I, J, K, I1, I2, TD
185 INTEGER LDTB, NB, KB, JB, NT, IINFO
186 COMPLEX*16 PIV
187
188
189 LOGICAL LSAME
190 INTEGER ILAENV
192
193
197
198
199 INTRINSIC dconjg, min, max
200
201
202
203
204
205 info = 0
206 upper =
lsame( uplo,
'U' )
207 wquery = ( lwork.EQ.-1 )
208 tquery = ( ltb.EQ.-1 )
209 IF( .NOT.upper .AND. .NOT.
lsame( uplo,
'L' ) )
THEN
210 info = -1
211 ELSE IF( n.LT.0 ) THEN
212 info = -2
213 ELSE IF( lda.LT.max( 1, n ) ) THEN
214 info = -4
215 ELSE IF ( ltb .LT. 4*n .AND. .NOT.tquery ) THEN
216 info = -6
217 ELSE IF ( lwork .LT. n .AND. .NOT.wquery ) THEN
218 info = -10
219 END IF
220
221 IF( info.NE.0 ) THEN
222 CALL xerbla(
'ZHETRF_AA_2STAGE', -info )
223 RETURN
224 END IF
225
226
227
228 nb =
ilaenv( 1,
'ZHETRF_AA_2STAGE', uplo, n, -1, -1, -1 )
229 IF( info.EQ.0 ) THEN
230 IF( tquery ) THEN
231 tb( 1 ) = (3*nb+1)*n
232 END IF
233 IF( wquery ) THEN
234 work( 1 ) = n*nb
235 END IF
236 END IF
237 IF( tquery .OR. wquery ) THEN
238 RETURN
239 END IF
240
241
242
243 IF ( n.EQ.0 ) THEN
244 RETURN
245 ENDIF
246
247
248
249 ldtb = ltb/n
250 IF( ldtb .LT. 3*nb+1 ) THEN
251 nb = (ldtb-1)/3
252 END IF
253 IF( lwork .LT. nb*n ) THEN
254 nb = lwork/n
255 END IF
256
257
258
259 nt = (n+nb-1)/nb
260 td = 2*nb
261 kb = min(nb, n)
262
263
264
265 DO j = 1, kb
266 ipiv( j ) = j
267 END DO
268
269
270
271 tb( 1 ) = nb
272
273 IF( upper ) THEN
274
275
276
277
278
279 DO j = 0, nt-1
280
281
282
283 kb = min(nb, n-j*nb)
284 DO i = 1, j-1
285 IF( i.EQ.1 ) THEN
286
287 IF( i .EQ. (j-1) ) THEN
288 jb = nb+kb
289 ELSE
290 jb = 2*nb
291 END IF
292 CALL zgemm(
'NoTranspose',
'NoTranspose',
293 $ nb, kb, jb,
294 $ one, tb( td+1 + (i*nb)*ldtb ), ldtb-1,
295 $ a( (i-1)*nb+1, j*nb+1 ), lda,
296 $ zero, work( i*nb+1 ), n )
297 ELSE
298
299 IF( i .EQ. (j-1) ) THEN
300 jb = 2*nb+kb
301 ELSE
302 jb = 3*nb
303 END IF
304 CALL zgemm(
'NoTranspose',
'NoTranspose',
305 $ nb, kb, jb,
306 $ one, tb( td+nb+1 + ((i-1)*nb)*ldtb ),
307 $ ldtb-1,
308 $ a( (i-2)*nb+1, j*nb+1 ), lda,
309 $ zero, work( i*nb+1 ), n )
310 END IF
311 END DO
312
313
314
315 CALL zlacpy(
'Upper', kb, kb, a( j*nb+1, j*nb+1 ), lda,
316 $ tb( td+1 + (j*nb)*ldtb ), ldtb-1 )
317 IF( j.GT.1 ) THEN
318
319 CALL zgemm(
'Conjugate transpose',
'NoTranspose',
320 $ kb, kb, (j-1)*nb,
321 $ -one, a( 1, j*nb+1 ), lda,
322 $ work( nb+1 ), n,
323 $ one, tb( td+1 + (j*nb)*ldtb ), ldtb-1 )
324
325 CALL zgemm(
'Conjugate transpose',
'NoTranspose',
326 $ kb, nb, kb,
327 $ one, a( (j-1)*nb+1, j*nb+1 ), lda,
328 $ tb( td+nb+1 + ((j-1)*nb)*ldtb ), ldtb-1,
329 $ zero, work( 1 ), n )
330 CALL zgemm(
'NoTranspose',
'NoTranspose',
331 $ kb, kb, nb,
332 $ -one, work( 1 ), n,
333 $ a( (j-2)*nb+1, j*nb+1 ), lda,
334 $ one, tb( td+1 + (j*nb)*ldtb ), ldtb-1 )
335 END IF
336 IF( j.GT.0 ) THEN
337 CALL zhegst( 1,
'Upper', kb,
338 $ tb( td+1 + (j*nb)*ldtb ), ldtb-1,
339 $ a( (j-1)*nb+1, j*nb+1 ), lda, iinfo )
340 END IF
341
342
343
344 DO i = 1, kb
345 tb( td+1 + (j*nb+i-1)*ldtb )
346 $ = real( tb( td+1 + (j*nb+i-1)*ldtb ) )
347 DO k = i+1, kb
348 tb( td+(k-i)+1 + (j*nb+i-1)*ldtb )
349 $ = dconjg( tb( td-(k-(i+1)) + (j*nb+k-1)*ldtb ) )
350 END DO
351 END DO
352
353 IF( j.LT.nt-1 ) THEN
354 IF( j.GT.0 ) THEN
355
356
357
358 IF( j.EQ.1 ) THEN
359 CALL zgemm(
'NoTranspose',
'NoTranspose',
360 $ kb, kb, kb,
361 $ one, tb( td+1 + (j*nb)*ldtb ), ldtb-1,
362 $ a( (j-1)*nb+1, j*nb+1 ), lda,
363 $ zero, work( j*nb+1 ), n )
364 ELSE
365 CALL zgemm(
'NoTranspose',
'NoTranspose',
366 $ kb, kb, nb+kb,
367 $ one, tb( td+nb+1 + ((j-1)*nb)*ldtb ),
368 $ ldtb-1,
369 $ a( (j-2)*nb+1, j*nb+1 ), lda,
370 $ zero, work( j*nb+1 ), n )
371 END IF
372
373
374
375 CALL zgemm(
'Conjugate transpose',
'NoTranspose',
376 $ nb, n-(j+1)*nb, j*nb,
377 $ -one, work( nb+1 ), n,
378 $ a( 1, (j+1)*nb+1 ), lda,
379 $ one, a( j*nb+1, (j+1)*nb+1 ), lda )
380 END IF
381
382
383
384 DO k = 1, nb
385 CALL zcopy( n-(j+1)*nb,
386 $ a( j*nb+k, (j+1)*nb+1 ), lda,
387 $ work( 1+(k-1)*n ), 1 )
388 END DO
389
390
391
392 CALL zgetrf( n-(j+1)*nb, nb,
393 $ work, n,
394 $ ipiv( (j+1)*nb+1 ), iinfo )
395
396
397
398
399
400
401 DO k = 1, nb
402
403
404
405 CALL zcopy( n-k-(j+1)*nb,
406 $ work( k+1+(k-1)*n ), 1,
407 $ a( j*nb+k, (j+1)*nb+k+1 ), lda )
408
409
410
411 CALL zlacgv( k, work( 1+(k-1)*n ), 1 )
412 END DO
413
414
415
416 kb = min(nb, n-(j+1)*nb)
417 CALL zlaset(
'Full', kb, nb, zero, zero,
418 $ tb( td+nb+1 + (j*nb)*ldtb) , ldtb-1 )
419 CALL zlacpy(
'Upper', kb, nb,
420 $ work, n,
421 $ tb( td+nb+1 + (j*nb)*ldtb ), ldtb-1 )
422 IF( j.GT.0 ) THEN
423 CALL ztrsm(
'R',
'U',
'N',
'U', kb, nb, one,
424 $ a( (j-1)*nb+1, j*nb+1 ), lda,
425 $ tb( td+nb+1 + (j*nb)*ldtb ), ldtb-1 )
426 END IF
427
428
429
430
431 DO k = 1, nb
432 DO i = 1, kb
433 tb( td-nb+k-i+1 + (j*nb+nb+i-1)*ldtb )
434 $ = dconjg( tb( td+nb+i-k+1 + (j*nb+k-1)*ldtb ) )
435 END DO
436 END DO
437 CALL zlaset(
'Lower', kb, nb, zero, one,
438 $ a( j*nb+1, (j+1)*nb+1), lda )
439
440
441
442 DO k = 1, kb
443
444 ipiv( (j+1)*nb+k ) = ipiv( (j+1)*nb+k ) + (j+1)*nb
445
446 i1 = (j+1)*nb+k
447 i2 = ipiv( (j+1)*nb+k )
448 IF( i1.NE.i2 ) THEN
449
450 CALL zswap( k-1, a( (j+1)*nb+1, i1 ), 1,
451 $ a( (j+1)*nb+1, i2 ), 1 )
452
453 IF( i2.GT.(i1+1) ) THEN
454 CALL zswap( i2-i1-1, a( i1, i1+1 ), lda,
455 $ a( i1+1, i2 ), 1 )
456 CALL zlacgv( i2-i1-1, a( i1+1, i2 ), 1 )
457 END IF
458 CALL zlacgv( i2-i1, a( i1, i1+1 ), lda )
459
460 IF( i2.LT.n )
461 $
CALL zswap( n-i2, a( i1, i2+1 ), lda,
462 $ a( i2, i2+1 ), lda )
463
464 piv = a( i1, i1 )
465 a( i1, i1 ) = a( i2, i2 )
466 a( i2, i2 ) = piv
467
468 IF( j.GT.0 ) THEN
469 CALL zswap( j*nb, a( 1, i1 ), 1,
470 $ a( 1, i2 ), 1 )
471 END IF
472 ENDIF
473 END DO
474 END IF
475 END DO
476 ELSE
477
478
479
480
481
482 DO j = 0, nt-1
483
484
485
486 kb = min(nb, n-j*nb)
487 DO i = 1, j-1
488 IF( i.EQ.1 ) THEN
489
490 IF( i .EQ. (j-1) ) THEN
491 jb = nb+kb
492 ELSE
493 jb = 2*nb
494 END IF
495 CALL zgemm(
'NoTranspose',
'Conjugate transpose',
496 $ nb, kb, jb,
497 $ one, tb( td+1 + (i*nb)*ldtb ), ldtb-1,
498 $ a( j*nb+1, (i-1)*nb+1 ), lda,
499 $ zero, work( i*nb+1 ), n )
500 ELSE
501
502 IF( i .EQ. (j-1) ) THEN
503 jb = 2*nb+kb
504 ELSE
505 jb = 3*nb
506 END IF
507 CALL zgemm(
'NoTranspose',
'Conjugate transpose',
508 $ nb, kb, jb,
509 $ one, tb( td+nb+1 + ((i-1)*nb)*ldtb ),
510 $ ldtb-1,
511 $ a( j*nb+1, (i-2)*nb+1 ), lda,
512 $ zero, work( i*nb+1 ), n )
513 END IF
514 END DO
515
516
517
518 CALL zlacpy(
'Lower', kb, kb, a( j*nb+1, j*nb+1 ), lda,
519 $ tb( td+1 + (j*nb)*ldtb ), ldtb-1 )
520 IF( j.GT.1 ) THEN
521
522 CALL zgemm(
'NoTranspose',
'NoTranspose',
523 $ kb, kb, (j-1)*nb,
524 $ -one, a( j*nb+1, 1 ), lda,
525 $ work( nb+1 ), n,
526 $ one, tb( td+1 + (j*nb)*ldtb ), ldtb-1 )
527
528 CALL zgemm(
'NoTranspose',
'NoTranspose',
529 $ kb, nb, kb,
530 $ one, a( j*nb+1, (j-1)*nb+1 ), lda,
531 $ tb( td+nb+1 + ((j-1)*nb)*ldtb ), ldtb-1,
532 $ zero, work( 1 ), n )
533 CALL zgemm(
'NoTranspose',
'Conjugate transpose',
534 $ kb, kb, nb,
535 $ -one, work( 1 ), n,
536 $ a( j*nb+1, (j-2)*nb+1 ), lda,
537 $ one, tb( td+1 + (j*nb)*ldtb ), ldtb-1 )
538 END IF
539 IF( j.GT.0 ) THEN
540 CALL zhegst( 1,
'Lower', kb,
541 $ tb( td+1 + (j*nb)*ldtb ), ldtb-1,
542 $ a( j*nb+1, (j-1)*nb+1 ), lda, iinfo )
543 END IF
544
545
546
547 DO i = 1, kb
548 tb( td+1 + (j*nb+i-1)*ldtb )
549 $ = real( tb( td+1 + (j*nb+i-1)*ldtb ) )
550 DO k = i+1, kb
551 tb( td-(k-(i+1)) + (j*nb+k-1)*ldtb )
552 $ = dconjg( tb( td+(k-i)+1 + (j*nb+i-1)*ldtb ) )
553 END DO
554 END DO
555
556 IF( j.LT.nt-1 ) THEN
557 IF( j.GT.0 ) THEN
558
559
560
561 IF( j.EQ.1 ) THEN
562 CALL zgemm(
'NoTranspose',
'Conjugate transpose',
563 $ kb, kb, kb,
564 $ one, tb( td+1 + (j*nb)*ldtb ), ldtb-1,
565 $ a( j*nb+1, (j-1)*nb+1 ), lda,
566 $ zero, work( j*nb+1 ), n )
567 ELSE
568 CALL zgemm(
'NoTranspose',
'Conjugate transpose',
569 $ kb, kb, nb+kb,
570 $ one, tb( td+nb+1 + ((j-1)*nb)*ldtb ),
571 $ ldtb-1,
572 $ a( j*nb+1, (j-2)*nb+1 ), lda,
573 $ zero, work( j*nb+1 ), n )
574 END IF
575
576
577
578 CALL zgemm(
'NoTranspose',
'NoTranspose',
579 $ n-(j+1)*nb, nb, j*nb,
580 $ -one, a( (j+1)*nb+1, 1 ), lda,
581 $ work( nb+1 ), n,
582 $ one, a( (j+1)*nb+1, j*nb+1 ), lda )
583 END IF
584
585
586
587 CALL zgetrf( n-(j+1)*nb, nb,
588 $ a( (j+1)*nb+1, j*nb+1 ), lda,
589 $ ipiv( (j+1)*nb+1 ), iinfo )
590
591
592
593
594
595
596 kb = min(nb, n-(j+1)*nb)
597 CALL zlaset(
'Full', kb, nb, zero, zero,
598 $ tb( td+nb+1 + (j*nb)*ldtb) , ldtb-1 )
599 CALL zlacpy(
'Upper', kb, nb,
600 $ a( (j+1)*nb+1, j*nb+1 ), lda,
601 $ tb( td+nb+1 + (j*nb)*ldtb ), ldtb-1 )
602 IF( j.GT.0 ) THEN
603 CALL ztrsm(
'R',
'L',
'C',
'U', kb, nb, one,
604 $ a( j*nb+1, (j-1)*nb+1 ), lda,
605 $ tb( td+nb+1 + (j*nb)*ldtb ), ldtb-1 )
606 END IF
607
608
609
610
611 DO k = 1, nb
612 DO i = 1, kb
613 tb( td-nb+k-i+1 + (j*nb+nb+i-1)*ldtb )
614 $ = dconjg( tb( td+nb+i-k+1 + (j*nb+k-1)*ldtb ) )
615 END DO
616 END DO
617 CALL zlaset(
'Upper', kb, nb, zero, one,
618 $ a( (j+1)*nb+1, j*nb+1), lda )
619
620
621
622 DO k = 1, kb
623
624 ipiv( (j+1)*nb+k ) = ipiv( (j+1)*nb+k ) + (j+1)*nb
625
626 i1 = (j+1)*nb+k
627 i2 = ipiv( (j+1)*nb+k )
628 IF( i1.NE.i2 ) THEN
629
630 CALL zswap( k-1, a( i1, (j+1)*nb+1 ), lda,
631 $ a( i2, (j+1)*nb+1 ), lda )
632
633 IF( i2.GT.(i1+1) ) THEN
634 CALL zswap( i2-i1-1, a( i1+1, i1 ), 1,
635 $ a( i2, i1+1 ), lda )
636 CALL zlacgv( i2-i1-1, a( i2, i1+1 ), lda )
637 END IF
638 CALL zlacgv( i2-i1, a( i1+1, i1 ), 1 )
639
640 IF( i2.LT.n )
641 $
CALL zswap( n-i2, a( i2+1, i1 ), 1,
642 $ a( i2+1, i2 ), 1 )
643
644 piv = a( i1, i1 )
645 a( i1, i1 ) = a( i2, i2 )
646 a( i2, i2 ) = piv
647
648 IF( j.GT.0 ) THEN
649 CALL zswap( j*nb, a( i1, 1 ), lda,
650 $ a( i2, 1 ), lda )
651 END IF
652 ENDIF
653 END DO
654
655
656
657
658
659 END IF
660 END DO
661 END IF
662
663
664 CALL zgbtrf( n, n, nb, nb, tb, ldtb, ipiv2, info )
665
666 RETURN
667
668
669
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV
subroutine xerbla(SRNAME, INFO)
XERBLA
logical function lsame(CA, CB)
LSAME
subroutine zswap(N, ZX, INCX, ZY, INCY)
ZSWAP
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
subroutine ztrsm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
ZTRSM
subroutine zgbtrf(M, N, KL, KU, AB, LDAB, IPIV, INFO)
ZGBTRF
subroutine zhegst(ITYPE, UPLO, N, A, LDA, B, LDB, INFO)
ZHEGST
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
subroutine zlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
subroutine zlacgv(N, X, INCX)
ZLACGV conjugates a complex vector.
subroutine zgetrf(M, N, A, LDA, IPIV, INFO)
ZGETRF VARIANT: Crout Level 3 BLAS version of the algorithm.