230 SUBROUTINE dgghd3( COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, Q,
231 $ ldq, z, ldz, work, lwork, info )
241 CHARACTER COMPQ, COMPZ
242 INTEGER IHI, ILO, INFO, LDA, LDB, LDQ, LDZ, N, LWORK
245 DOUBLE PRECISION A( lda, * ), B( ldb, * ), Q( ldq, * ),
246 $ z( ldz, * ), work( * )
252 DOUBLE PRECISION ZERO, ONE
253 parameter ( zero = 0.0d+0, one = 1.0d+0 )
256 LOGICAL BLK22, INITQ, INITZ, LQUERY, WANTQ, WANTZ
257 CHARACTER*1 COMPQ2, COMPZ2
258 INTEGER COLA, I, IERR, J, J0, JCOL, JJ, JROW, K,
259 $ kacc22, len, lwkopt, n2nb, nb, nblst, nbmin,
260 $ nh, nnb, nx, ppw, ppwo, pw, top, topq
261 DOUBLE PRECISION C, C1, C2, S, S1, S2, TEMP, TEMP1, TEMP2, TEMP3
266 EXTERNAL ilaenv, lsame
279 nb = ilaenv( 1,
'DGGHD3',
' ', n, ilo, ihi, -1 )
280 lwkopt = max( 6*n*nb, 1 )
281 work( 1 ) = dble( lwkopt )
282 initq = lsame( compq,
'I' )
283 wantq = initq .OR. lsame( compq,
'V' )
284 initz = lsame( compz,
'I' )
285 wantz = initz .OR. lsame( compz,
'V' )
286 lquery = ( lwork.EQ.-1 )
288 IF( .NOT.lsame( compq,
'N' ) .AND. .NOT.wantq )
THEN
290 ELSE IF( .NOT.lsame( compz,
'N' ) .AND. .NOT.wantz )
THEN
292 ELSE IF( n.LT.0 )
THEN
294 ELSE IF( ilo.LT.1 )
THEN
296 ELSE IF( ihi.GT.n .OR. ihi.LT.ilo-1 )
THEN
298 ELSE IF( lda.LT.max( 1, n ) )
THEN
300 ELSE IF( ldb.LT.max( 1, n ) )
THEN
302 ELSE IF( ( wantq .AND. ldq.LT.n ) .OR. ldq.LT.1 )
THEN
304 ELSE IF( ( wantz .AND. ldz.LT.n ) .OR. ldz.LT.1 )
THEN
306 ELSE IF( lwork.LT.1 .AND. .NOT.lquery )
THEN
310 CALL xerbla(
'DGGHD3', -info )
312 ELSE IF( lquery )
THEN
319 $
CALL dlaset(
'All', n, n, zero, one, q, ldq )
321 $
CALL dlaset(
'All', n, n, zero, one, z, ldz )
326 $
CALL dlaset(
'Lower', n-1, n-1, zero, zero, b(2, 1), ldb )
338 nbmin = ilaenv( 2,
'DGGHD3',
' ', n, ilo, ihi, -1 )
339 IF( nb.GT.1 .AND. nb.LT.nh )
THEN
343 nx = max( nb, ilaenv( 3,
'DGGHD3',
' ', n, ilo, ihi, -1 ) )
348 IF( lwork.LT.lwkopt )
THEN
354 nbmin = max( 2, ilaenv( 2,
'DGGHD3',
' ', n, ilo, ihi,
356 IF( lwork.GE.6*n*nbmin )
THEN
365 IF( nb.LT.nbmin .OR. nb.GE.nh )
THEN
375 kacc22 = ilaenv( 16,
'DGGHD3',
' ', n, ilo, ihi, -1 )
377 DO jcol = ilo, ihi-2, nb
378 nnb = min( nb, ihi-jcol-1 )
386 n2nb = ( ihi-jcol-1 ) / nnb - 1
387 nblst = ihi - jcol - n2nb*nnb
388 CALL dlaset(
'All', nblst, nblst, zero, one, work, nblst )
389 pw = nblst * nblst + 1
391 CALL dlaset(
'All', 2*nnb, 2*nnb, zero, one,
392 $ work( pw ), 2*nnb )
398 DO j = jcol, jcol+nnb-1
405 CALL dlartg( temp, a( i, j ), c, s, a( i-1, j ) )
412 ppw = ( nblst + 1 )*( nblst - 2 ) - j + jcol + 1
414 jrow = j + n2nb*nnb + 2
418 DO jj = ppw, ppw+len-1
419 temp = work( jj + nblst )
420 work( jj + nblst ) = c*temp - s*work( jj )
421 work( jj ) = s*temp + c*work( jj )
424 ppw = ppw - nblst - 1
427 ppwo = nblst*nblst + ( nnb+j-jcol-1 )*2*nnb + nnb
429 DO jrow = j0, j+2, -nnb
432 DO i = jrow+nnb-1, jrow, -1
435 DO jj = ppw, ppw+len-1
436 temp = work( jj + 2*nnb )
437 work( jj + 2*nnb ) = c*temp - s*work( jj )
438 work( jj ) = s*temp + c*work( jj )
441 ppw = ppw - 2*nnb - 1
443 ppwo = ppwo + 4*nnb*nnb
462 DO i = min( jj+1, ihi ), j+2, -1
466 b( i, jj ) = c*temp - s*b( i-1, jj )
467 b( i-1, jj ) = s*temp + c*b( i-1, jj )
473 temp = b( jj+1, jj+1 )
474 CALL dlartg( temp, b( jj+1, jj ), c, s,
477 CALL drot( jj-top, b( top+1, jj+1 ), 1,
478 $ b( top+1, jj ), 1, c, s )
491 jj = mod( ihi-j-1, 3 )
492 DO i = ihi-j-3, jj+1, -3
502 temp1 = a( k, j+i+1 )
503 temp2 = a( k, j+i+2 )
504 temp3 = a( k, j+i+3 )
505 a( k, j+i+3 ) = c2*temp3 + s2*temp2
506 temp2 = -s2*temp3 + c2*temp2
507 a( k, j+i+2 ) = c1*temp2 + s1*temp1
508 temp1 = -s1*temp2 + c1*temp1
509 a( k, j+i+1 ) = c*temp1 + s*temp
510 a( k, j+i ) = -s*temp1 + c*temp
516 CALL drot( ihi-top, a( top+1, j+i+1 ), 1,
517 $ a( top+1, j+i ), 1, a( j+1+i, j ),
524 IF ( j .LT. jcol + nnb - 1 )
THEN
537 jrow = ihi - nblst + 1
538 CALL dgemv(
'Transpose', nblst, len, one, work,
539 $ nblst, a( jrow, j+1 ), 1, zero,
542 DO i = jrow, jrow+nblst-len-1
543 work( ppw ) = a( i, j+1 )
546 CALL dtrmv(
'Lower',
'Transpose',
'Non-unit',
547 $ nblst-len, work( len*nblst + 1 ), nblst,
548 $ work( pw+len ), 1 )
549 CALL dgemv(
'Transpose', len, nblst-len, one,
550 $ work( (len+1)*nblst - len + 1 ), nblst,
551 $ a( jrow+nblst-len, j+1 ), 1, one,
552 $ work( pw+len ), 1 )
554 DO i = jrow, jrow+nblst-1
555 a( i, j+1 ) = work( ppw )
572 ppwo = 1 + nblst*nblst
574 DO jrow = j0, jcol+1, -nnb
576 DO i = jrow, jrow+nnb-1
577 work( ppw ) = a( i, j+1 )
581 DO i = jrow+nnb, jrow+nnb+len-1
582 work( ppw ) = a( i, j+1 )
585 CALL dtrmv(
'Upper',
'Transpose',
'Non-unit', len,
586 $ work( ppwo + nnb ), 2*nnb, work( pw ),
588 CALL dtrmv(
'Lower',
'Transpose',
'Non-unit', nnb,
589 $ work( ppwo + 2*len*nnb ),
590 $ 2*nnb, work( pw + len ), 1 )
591 CALL dgemv(
'Transpose', nnb, len, one,
592 $ work( ppwo ), 2*nnb, a( jrow, j+1 ), 1,
593 $ one, work( pw ), 1 )
594 CALL dgemv(
'Transpose', len, nnb, one,
595 $ work( ppwo + 2*len*nnb + nnb ), 2*nnb,
596 $ a( jrow+nnb, j+1 ), 1, one,
597 $ work( pw+len ), 1 )
599 DO i = jrow, jrow+len+nnb-1
600 a( i, j+1 ) = work( ppw )
603 ppwo = ppwo + 4*nnb*nnb
610 cola = n - jcol - nnb + 1
612 CALL dgemm(
'Transpose',
'No Transpose', nblst,
613 $ cola, nblst, one, work, nblst,
614 $ a( j, jcol+nnb ), lda, zero, work( pw ),
616 CALL dlacpy(
'All', nblst, cola, work( pw ), nblst,
617 $ a( j, jcol+nnb ), lda )
618 ppwo = nblst*nblst + 1
620 DO j = j0, jcol+1, -nnb
632 CALL dorm22(
'Left',
'Transpose', 2*nnb, cola, nnb,
633 $ nnb, work( ppwo ), 2*nnb,
634 $ a( j, jcol+nnb ), lda, work( pw ),
640 CALL dgemm(
'Transpose',
'No Transpose', 2*nnb,
641 $ cola, 2*nnb, one, work( ppwo ), 2*nnb,
642 $ a( j, jcol+nnb ), lda, zero, work( pw ),
644 CALL dlacpy(
'All', 2*nnb, cola, work( pw ), 2*nnb,
645 $ a( j, jcol+nnb ), lda )
647 ppwo = ppwo + 4*nnb*nnb
655 topq = max( 2, j - jcol + 1 )
661 CALL dgemm(
'No Transpose',
'No Transpose', nh,
662 $ nblst, nblst, one, q( topq, j ), ldq,
663 $ work, nblst, zero, work( pw ), nh )
664 CALL dlacpy(
'All', nh, nblst, work( pw ), nh,
665 $ q( topq, j ), ldq )
666 ppwo = nblst*nblst + 1
668 DO j = j0, jcol+1, -nnb
670 topq = max( 2, j - jcol + 1 )
677 CALL dorm22(
'Right',
'No Transpose', nh, 2*nnb,
678 $ nnb, nnb, work( ppwo ), 2*nnb,
679 $ q( topq, j ), ldq, work( pw ),
685 CALL dgemm(
'No Transpose',
'No Transpose', nh,
686 $ 2*nnb, 2*nnb, one, q( topq, j ), ldq,
687 $ work( ppwo ), 2*nnb, zero, work( pw ),
689 CALL dlacpy(
'All', nh, 2*nnb, work( pw ), nh,
690 $ q( topq, j ), ldq )
692 ppwo = ppwo + 4*nnb*nnb
698 IF ( wantz .OR. top.GT.0 )
THEN
703 CALL dlaset(
'All', nblst, nblst, zero, one, work,
705 pw = nblst * nblst + 1
707 CALL dlaset(
'All', 2*nnb, 2*nnb, zero, one,
708 $ work( pw ), 2*nnb )
714 DO j = jcol, jcol+nnb-1
715 ppw = ( nblst + 1 )*( nblst - 2 ) - j + jcol + 1
717 jrow = j + n2nb*nnb + 2
723 DO jj = ppw, ppw+len-1
724 temp = work( jj + nblst )
725 work( jj + nblst ) = c*temp - s*work( jj )
726 work( jj ) = s*temp + c*work( jj )
729 ppw = ppw - nblst - 1
732 ppwo = nblst*nblst + ( nnb+j-jcol-1 )*2*nnb + nnb
734 DO jrow = j0, j+2, -nnb
737 DO i = jrow+nnb-1, jrow, -1
742 DO jj = ppw, ppw+len-1
743 temp = work( jj + 2*nnb )
744 work( jj + 2*nnb ) = c*temp - s*work( jj )
745 work( jj ) = s*temp + c*work( jj )
748 ppw = ppw - 2*nnb - 1
750 ppwo = ppwo + 4*nnb*nnb
755 CALL dlaset(
'Lower', ihi - jcol - 1, nnb, zero, zero,
756 $ a( jcol + 2, jcol ), lda )
757 CALL dlaset(
'Lower', ihi - jcol - 1, nnb, zero, zero,
758 $ b( jcol + 2, jcol ), ldb )
765 CALL dgemm(
'No Transpose',
'No Transpose', top,
766 $ nblst, nblst, one, a( 1, j ), lda,
767 $ work, nblst, zero, work( pw ), top )
768 CALL dlacpy(
'All', top, nblst, work( pw ), top,
770 ppwo = nblst*nblst + 1
772 DO j = j0, jcol+1, -nnb
777 CALL dorm22(
'Right',
'No Transpose', top, 2*nnb,
778 $ nnb, nnb, work( ppwo ), 2*nnb,
779 $ a( 1, j ), lda, work( pw ),
785 CALL dgemm(
'No Transpose',
'No Transpose', top,
786 $ 2*nnb, 2*nnb, one, a( 1, j ), lda,
787 $ work( ppwo ), 2*nnb, zero,
789 CALL dlacpy(
'All', top, 2*nnb, work( pw ), top,
792 ppwo = ppwo + 4*nnb*nnb
796 CALL dgemm(
'No Transpose',
'No Transpose', top,
797 $ nblst, nblst, one, b( 1, j ), ldb,
798 $ work, nblst, zero, work( pw ), top )
799 CALL dlacpy(
'All', top, nblst, work( pw ), top,
801 ppwo = nblst*nblst + 1
803 DO j = j0, jcol+1, -nnb
808 CALL dorm22(
'Right',
'No Transpose', top, 2*nnb,
809 $ nnb, nnb, work( ppwo ), 2*nnb,
810 $ b( 1, j ), ldb, work( pw ),
816 CALL dgemm(
'No Transpose',
'No Transpose', top,
817 $ 2*nnb, 2*nnb, one, b( 1, j ), ldb,
818 $ work( ppwo ), 2*nnb, zero,
820 CALL dlacpy(
'All', top, 2*nnb, work( pw ), top,
823 ppwo = ppwo + 4*nnb*nnb
832 topq = max( 2, j - jcol + 1 )
838 CALL dgemm(
'No Transpose',
'No Transpose', nh,
839 $ nblst, nblst, one, z( topq, j ), ldz,
840 $ work, nblst, zero, work( pw ), nh )
841 CALL dlacpy(
'All', nh, nblst, work( pw ), nh,
842 $ z( topq, j ), ldz )
843 ppwo = nblst*nblst + 1
845 DO j = j0, jcol+1, -nnb
847 topq = max( 2, j - jcol + 1 )
854 CALL dorm22(
'Right',
'No Transpose', nh, 2*nnb,
855 $ nnb, nnb, work( ppwo ), 2*nnb,
856 $ z( topq, j ), ldz, work( pw ),
862 CALL dgemm(
'No Transpose',
'No Transpose', nh,
863 $ 2*nnb, 2*nnb, one, z( topq, j ), ldz,
864 $ work( ppwo ), 2*nnb, zero, work( pw ),
866 CALL dlacpy(
'All', nh, 2*nnb, work( pw ), nh,
867 $ z( topq, j ), ldz )
869 ppwo = ppwo + 4*nnb*nnb
880 IF ( jcol.NE.ilo )
THEN
888 $
CALL dgghrd( compq2, compz2, n, jcol, ihi, a, lda, b, ldb, q,
889 $ ldq, z, ldz, ierr )
890 work( 1 ) = dble( lwkopt )
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
subroutine dorm22(SIDE, TRANS, M, N, N1, N2, Q, LDQ, C, LDC, WORK, LWORK, INFO)
DORM22 multiplies a general matrix by a banded orthogonal matrix.
subroutine dgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
DGEMV
subroutine drot(N, DX, INCX, DY, INCY, C, S)
DROT
subroutine dlartg(F, G, CS, SN, R)
DLARTG generates a plane rotation with real cosine and real sine.
subroutine dlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine dgghd3(COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, Q, LDQ, Z, LDZ, WORK, LWORK, INFO)
DGGHD3
subroutine dtrmv(UPLO, TRANS, DIAG, N, A, LDA, X, INCX)
DTRMV
subroutine dgghrd(COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, Q, LDQ, Z, LDZ, INFO)
DGGHRD