228 SUBROUTINE sgghd3( COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, Q,
229 $ LDQ, Z, LDZ, WORK, LWORK, INFO )
238 CHARACTER COMPQ, COMPZ
239 INTEGER IHI, ILO, INFO, LDA, LDB, LDQ, LDZ, N, LWORK
242 REAL A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
243 $ z( ldz, * ), work( * )
250 parameter( zero = 0.0e+0, one = 1.0e+0 )
253 LOGICAL BLK22, INITQ, INITZ, LQUERY, WANTQ, WANTZ
254 CHARACTER*1 COMPQ2, COMPZ2
255 INTEGER COLA, I, IERR, J, J0, JCOL, JJ, JROW, K,
256 $ kacc22, len, lwkopt, n2nb, nb, nblst, nbmin,
257 $ nh, nnb, nx, ppw, ppwo, pw, top, topq
258 REAL C, C1, C2, S, S1, S2, TEMP, TEMP1, TEMP2, TEMP3
263 EXTERNAL ilaenv, lsame
277 nb = ilaenv( 1,
'SGGHD3',
' ', n, ilo, ihi, -1 )
278 lwkopt = max( 6*n*nb, 1 )
279 work( 1 ) = real( lwkopt )
280 initq = lsame( compq,
'I' )
281 wantq = initq .OR. lsame( compq,
'V' )
282 initz = lsame( compz,
'I' )
283 wantz = initz .OR. lsame( compz,
'V' )
284 lquery = ( lwork.EQ.-1 )
286 IF( .NOT.lsame( compq,
'N' ) .AND. .NOT.wantq )
THEN
288 ELSE IF( .NOT.lsame( compz,
'N' ) .AND. .NOT.wantz )
THEN
290 ELSE IF( n.LT.0 )
THEN
292 ELSE IF( ilo.LT.1 )
THEN
294 ELSE IF( ihi.GT.n .OR. ihi.LT.ilo-1 )
THEN
296 ELSE IF( lda.LT.max( 1, n ) )
THEN
298 ELSE IF( ldb.LT.max( 1, n ) )
THEN
300 ELSE IF( ( wantq .AND. ldq.LT.n ) .OR. ldq.LT.1 )
THEN
302 ELSE IF( ( wantz .AND. ldz.LT.n ) .OR. ldz.LT.1 )
THEN
304 ELSE IF( lwork.LT.1 .AND. .NOT.lquery )
THEN
308 CALL xerbla(
'SGGHD3', -info )
310 ELSE IF( lquery )
THEN
317 $
CALL slaset(
'All', n, n, zero, one, q, ldq )
319 $
CALL slaset(
'All', n, n, zero, one, z, ldz )
324 $
CALL slaset(
'Lower', n-1, n-1, zero, zero, b(2, 1), ldb )
336 nbmin = ilaenv( 2,
'SGGHD3',
' ', n, ilo, ihi, -1 )
337 IF( nb.GT.1 .AND. nb.LT.nh )
THEN
341 nx = max( nb, ilaenv( 3,
'SGGHD3',
' ', n, ilo, ihi, -1 ) )
346 IF( lwork.LT.lwkopt )
THEN
352 nbmin = max( 2, ilaenv( 2,
'SGGHD3',
' ', n, ilo, ihi,
354 IF( lwork.GE.6*n*nbmin )
THEN
363 IF( nb.LT.nbmin .OR. nb.GE.nh )
THEN
373 kacc22 = ilaenv( 16,
'SGGHD3',
' ', n, ilo, ihi, -1 )
375 DO jcol = ilo, ihi-2, nb
376 nnb = min( nb, ihi-jcol-1 )
384 n2nb = ( ihi-jcol-1 ) / nnb - 1
385 nblst = ihi - jcol - n2nb*nnb
386 CALL slaset(
'All', nblst, nblst, zero, one, work, nblst )
387 pw = nblst * nblst + 1
389 CALL slaset(
'All', 2*nnb, 2*nnb, zero, one,
390 $ work( pw ), 2*nnb )
396 DO j = jcol, jcol+nnb-1
403 CALL slartg( temp, a( i, j ), c, s, a( i-1, j ) )
410 ppw = ( nblst + 1 )*( nblst - 2 ) - j + jcol + 1
412 jrow = j + n2nb*nnb + 2
416 DO jj = ppw, ppw+len-1
417 temp = work( jj + nblst )
418 work( jj + nblst ) = c*temp - s*work( jj )
419 work( jj ) = s*temp + c*work( jj )
422 ppw = ppw - nblst - 1
425 ppwo = nblst*nblst + ( nnb+j-jcol-1 )*2*nnb + nnb
427 DO jrow = j0, j+2, -nnb
430 DO i = jrow+nnb-1, jrow, -1
433 DO jj = ppw, ppw+len-1
434 temp = work( jj + 2*nnb )
435 work( jj + 2*nnb ) = c*temp - s*work( jj )
436 work( jj ) = s*temp + c*work( jj )
439 ppw = ppw - 2*nnb - 1
441 ppwo = ppwo + 4*nnb*nnb
460 DO i = min( jj+1, ihi ), j+2, -1
464 b( i, jj ) = c*temp - s*b( i-1, jj )
465 b( i-1, jj ) = s*temp + c*b( i-1, jj )
471 temp = b( jj+1, jj+1 )
472 CALL slartg( temp, b( jj+1, jj ), c, s,
475 CALL srot( jj-top, b( top+1, jj+1 ), 1,
476 $ b( top+1, jj ), 1, c, s )
489 jj = mod( ihi-j-1, 3 )
490 DO i = ihi-j-3, jj+1, -3
500 temp1 = a( k, j+i+1 )
501 temp2 = a( k, j+i+2 )
502 temp3 = a( k, j+i+3 )
503 a( k, j+i+3 ) = c2*temp3 + s2*temp2
504 temp2 = -s2*temp3 + c2*temp2
505 a( k, j+i+2 ) = c1*temp2 + s1*temp1
506 temp1 = -s1*temp2 + c1*temp1
507 a( k, j+i+1 ) = c*temp1 + s*temp
508 a( k, j+i ) = -s*temp1 + c*temp
514 CALL srot( ihi-top, a( top+1, j+i+1 ), 1,
515 $ a( top+1, j+i ), 1, a( j+1+i, j ),
522 IF ( j .LT. jcol + nnb - 1 )
THEN
535 jrow = ihi - nblst + 1
536 CALL sgemv(
'Transpose', nblst, len, one, work,
537 $ nblst, a( jrow, j+1 ), 1, zero,
540 DO i = jrow, jrow+nblst-len-1
541 work( ppw ) = a( i, j+1 )
544 CALL strmv(
'Lower',
'Transpose',
'Non-unit',
545 $ nblst-len, work( len*nblst + 1 ), nblst,
546 $ work( pw+len ), 1 )
547 CALL sgemv(
'Transpose', len, nblst-len, one,
548 $ work( (len+1)*nblst - len + 1 ), nblst,
549 $ a( jrow+nblst-len, j+1 ), 1, one,
550 $ work( pw+len ), 1 )
552 DO i = jrow, jrow+nblst-1
553 a( i, j+1 ) = work( ppw )
570 ppwo = 1 + nblst*nblst
572 DO jrow = j0, jcol+1, -nnb
574 DO i = jrow, jrow+nnb-1
575 work( ppw ) = a( i, j+1 )
579 DO i = jrow+nnb, jrow+nnb+len-1
580 work( ppw ) = a( i, j+1 )
583 CALL strmv(
'Upper',
'Transpose',
'Non-unit', len,
584 $ work( ppwo + nnb ), 2*nnb, work( pw ),
586 CALL strmv(
'Lower',
'Transpose',
'Non-unit', nnb,
587 $ work( ppwo + 2*len*nnb ),
588 $ 2*nnb, work( pw + len ), 1 )
589 CALL sgemv(
'Transpose', nnb, len, one,
590 $ work( ppwo ), 2*nnb, a( jrow, j+1 ), 1,
591 $ one, work( pw ), 1 )
592 CALL sgemv(
'Transpose', len, nnb, one,
593 $ work( ppwo + 2*len*nnb + nnb ), 2*nnb,
594 $ a( jrow+nnb, j+1 ), 1, one,
595 $ work( pw+len ), 1 )
597 DO i = jrow, jrow+len+nnb-1
598 a( i, j+1 ) = work( ppw )
601 ppwo = ppwo + 4*nnb*nnb
608 cola = n - jcol - nnb + 1
610 CALL sgemm(
'Transpose',
'No Transpose', nblst,
611 $ cola, nblst, one, work, nblst,
612 $ a( j, jcol+nnb ), lda, zero, work( pw ),
614 CALL slacpy(
'All', nblst, cola, work( pw ), nblst,
615 $ a( j, jcol+nnb ), lda )
616 ppwo = nblst*nblst + 1
618 DO j = j0, jcol+1, -nnb
630 CALL sorm22(
'Left',
'Transpose', 2*nnb, cola, nnb,
631 $ nnb, work( ppwo ), 2*nnb,
632 $ a( j, jcol+nnb ), lda, work( pw ),
638 CALL sgemm(
'Transpose',
'No Transpose', 2*nnb,
639 $ cola, 2*nnb, one, work( ppwo ), 2*nnb,
640 $ a( j, jcol+nnb ), lda, zero, work( pw ),
642 CALL slacpy(
'All', 2*nnb, cola, work( pw ), 2*nnb,
643 $ a( j, jcol+nnb ), lda )
645 ppwo = ppwo + 4*nnb*nnb
653 topq = max( 2, j - jcol + 1 )
659 CALL sgemm(
'No Transpose',
'No Transpose', nh,
660 $ nblst, nblst, one, q( topq, j ), ldq,
661 $ work, nblst, zero, work( pw ), nh )
662 CALL slacpy(
'All', nh, nblst, work( pw ), nh,
663 $ q( topq, j ), ldq )
664 ppwo = nblst*nblst + 1
666 DO j = j0, jcol+1, -nnb
668 topq = max( 2, j - jcol + 1 )
675 CALL sorm22(
'Right',
'No Transpose', nh, 2*nnb,
676 $ nnb, nnb, work( ppwo ), 2*nnb,
677 $ q( topq, j ), ldq, work( pw ),
683 CALL sgemm(
'No Transpose',
'No Transpose', nh,
684 $ 2*nnb, 2*nnb, one, q( topq, j ), ldq,
685 $ work( ppwo ), 2*nnb, zero, work( pw ),
687 CALL slacpy(
'All', nh, 2*nnb, work( pw ), nh,
688 $ q( topq, j ), ldq )
690 ppwo = ppwo + 4*nnb*nnb
696 IF ( wantz .OR. top.GT.0 )
THEN
701 CALL slaset(
'All', nblst, nblst, zero, one, work,
703 pw = nblst * nblst + 1
705 CALL slaset(
'All', 2*nnb, 2*nnb, zero, one,
706 $ work( pw ), 2*nnb )
712 DO j = jcol, jcol+nnb-1
713 ppw = ( nblst + 1 )*( nblst - 2 ) - j + jcol + 1
715 jrow = j + n2nb*nnb + 2
721 DO jj = ppw, ppw+len-1
722 temp = work( jj + nblst )
723 work( jj + nblst ) = c*temp - s*work( jj )
724 work( jj ) = s*temp + c*work( jj )
727 ppw = ppw - nblst - 1
730 ppwo = nblst*nblst + ( nnb+j-jcol-1 )*2*nnb + nnb
732 DO jrow = j0, j+2, -nnb
735 DO i = jrow+nnb-1, jrow, -1
740 DO jj = ppw, ppw+len-1
741 temp = work( jj + 2*nnb )
742 work( jj + 2*nnb ) = c*temp - s*work( jj )
743 work( jj ) = s*temp + c*work( jj )
746 ppw = ppw - 2*nnb - 1
748 ppwo = ppwo + 4*nnb*nnb
753 CALL slaset(
'Lower', ihi - jcol - 1, nnb, zero, zero,
754 $ a( jcol + 2, jcol ), lda )
755 CALL slaset(
'Lower', ihi - jcol - 1, nnb, zero, zero,
756 $ b( jcol + 2, jcol ), ldb )
763 CALL sgemm(
'No Transpose',
'No Transpose', top,
764 $ nblst, nblst, one, a( 1, j ), lda,
765 $ work, nblst, zero, work( pw ), top )
766 CALL slacpy(
'All', top, nblst, work( pw ), top,
768 ppwo = nblst*nblst + 1
770 DO j = j0, jcol+1, -nnb
775 CALL sorm22(
'Right',
'No Transpose', top, 2*nnb,
776 $ nnb, nnb, work( ppwo ), 2*nnb,
777 $ a( 1, j ), lda, work( pw ),
783 CALL sgemm(
'No Transpose',
'No Transpose', top,
784 $ 2*nnb, 2*nnb, one, a( 1, j ), lda,
785 $ work( ppwo ), 2*nnb, zero,
787 CALL slacpy(
'All', top, 2*nnb, work( pw ), top,
790 ppwo = ppwo + 4*nnb*nnb
794 CALL sgemm(
'No Transpose',
'No Transpose', top,
795 $ nblst, nblst, one, b( 1, j ), ldb,
796 $ work, nblst, zero, work( pw ), top )
797 CALL slacpy(
'All', top, nblst, work( pw ), top,
799 ppwo = nblst*nblst + 1
801 DO j = j0, jcol+1, -nnb
806 CALL sorm22(
'Right',
'No Transpose', top, 2*nnb,
807 $ nnb, nnb, work( ppwo ), 2*nnb,
808 $ b( 1, j ), ldb, work( pw ),
814 CALL sgemm(
'No Transpose',
'No Transpose', top,
815 $ 2*nnb, 2*nnb, one, b( 1, j ), ldb,
816 $ work( ppwo ), 2*nnb, zero,
818 CALL slacpy(
'All', top, 2*nnb, work( pw ), top,
821 ppwo = ppwo + 4*nnb*nnb
830 topq = max( 2, j - jcol + 1 )
836 CALL sgemm(
'No Transpose',
'No Transpose', nh,
837 $ nblst, nblst, one, z( topq, j ), ldz,
838 $ work, nblst, zero, work( pw ), nh )
839 CALL slacpy(
'All', nh, nblst, work( pw ), nh,
840 $ z( topq, j ), ldz )
841 ppwo = nblst*nblst + 1
843 DO j = j0, jcol+1, -nnb
845 topq = max( 2, j - jcol + 1 )
852 CALL sorm22(
'Right',
'No Transpose', nh, 2*nnb,
853 $ nnb, nnb, work( ppwo ), 2*nnb,
854 $ z( topq, j ), ldz, work( pw ),
860 CALL sgemm(
'No Transpose',
'No Transpose', nh,
861 $ 2*nnb, 2*nnb, one, z( topq, j ), ldz,
862 $ work( ppwo ), 2*nnb, zero, work( pw ),
864 CALL slacpy(
'All', nh, 2*nnb, work( pw ), nh,
865 $ z( topq, j ), ldz )
867 ppwo = ppwo + 4*nnb*nnb
878 IF ( jcol.NE.ilo )
THEN
886 $
CALL sgghrd( compq2, compz2, n, jcol, ihi, a, lda, b, ldb, q,
887 $ ldq, z, ldz, ierr )
888 work( 1 ) = real( lwkopt )
subroutine slaset(UPLO, M, N, ALPHA, BETA, A, LDA)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
subroutine slartg(f, g, c, s, r)
SLARTG generates a plane rotation with real cosine and real sine.
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine sorm22(SIDE, TRANS, M, N, N1, N2, Q, LDQ, C, LDC, WORK, LWORK, INFO)
SORM22 multiplies a general matrix by a banded orthogonal matrix.
subroutine sgghd3(COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, Q, LDQ, Z, LDZ, WORK, LWORK, INFO)
SGGHD3
subroutine sgghrd(COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, Q, LDQ, Z, LDZ, INFO)
SGGHRD
subroutine srot(N, SX, INCX, SY, INCY, C, S)
SROT
subroutine strmv(UPLO, TRANS, DIAG, N, A, LDA, X, INCX)
STRMV
subroutine sgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
SGEMV
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM