230 SUBROUTINE sgghd3( 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 REAL A( lda, * ), B( ldb, * ), Q( ldq, * ),
246 $ z( ldz, * ), work( * )
253 parameter ( zero = 0.0e+0, one = 1.0e+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 REAL C, C1, C2, S, S1, S2, TEMP, TEMP1, TEMP2, TEMP3
266 EXTERNAL ilaenv, lsame
279 nb = ilaenv( 1,
'SGGHD3',
' ', n, ilo, ihi, -1 )
280 lwkopt = max( 6*n*nb, 1 )
281 work( 1 ) =
REAL( 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(
'SGGHD3', -info )
312 ELSE IF( lquery )
THEN
319 $
CALL slaset(
'All', n, n, zero, one, q, ldq )
321 $
CALL slaset(
'All', n, n, zero, one, z, ldz )
326 $
CALL slaset(
'Lower', n-1, n-1, zero, zero, b(2, 1), ldb )
338 nbmin = ilaenv( 2,
'SGGHD3',
' ', n, ilo, ihi, -1 )
339 IF( nb.GT.1 .AND. nb.LT.nh )
THEN
343 nx = max( nb, ilaenv( 3,
'SGGHD3',
' ', n, ilo, ihi, -1 ) )
348 IF( lwork.LT.lwkopt )
THEN
354 nbmin = max( 2, ilaenv( 2,
'SGGHD3',
' ', 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,
'SGGHD3',
' ', 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 slaset(
'All', nblst, nblst, zero, one, work, nblst )
389 pw = nblst * nblst + 1
391 CALL slaset(
'All', 2*nnb, 2*nnb, zero, one,
392 $ work( pw ), 2*nnb )
398 DO j = jcol, jcol+nnb-1
405 CALL slartg( 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 slartg( temp, b( jj+1, jj ), c, s,
477 CALL srot( 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 srot( 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 sgemv(
'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 strmv(
'Lower',
'Transpose',
'Non-unit',
547 $ nblst-len, work( len*nblst + 1 ), nblst,
548 $ work( pw+len ), 1 )
549 CALL sgemv(
'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 strmv(
'Upper',
'Transpose',
'Non-unit', len,
586 $ work( ppwo + nnb ), 2*nnb, work( pw ),
588 CALL strmv(
'Lower',
'Transpose',
'Non-unit', nnb,
589 $ work( ppwo + 2*len*nnb ),
590 $ 2*nnb, work( pw + len ), 1 )
591 CALL sgemv(
'Transpose', nnb, len, one,
592 $ work( ppwo ), 2*nnb, a( jrow, j+1 ), 1,
593 $ one, work( pw ), 1 )
594 CALL sgemv(
'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 sgemm(
'Transpose',
'No Transpose', nblst,
613 $ cola, nblst, one, work, nblst,
614 $ a( j, jcol+nnb ), lda, zero, work( pw ),
616 CALL slacpy(
'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 sorm22(
'Left',
'Transpose', 2*nnb, cola, nnb,
633 $ nnb, work( ppwo ), 2*nnb,
634 $ a( j, jcol+nnb ), lda, work( pw ),
640 CALL sgemm(
'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 slacpy(
'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 sgemm(
'No Transpose',
'No Transpose', nh,
662 $ nblst, nblst, one, q( topq, j ), ldq,
663 $ work, nblst, zero, work( pw ), nh )
664 CALL slacpy(
'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 sorm22(
'Right',
'No Transpose', nh, 2*nnb,
678 $ nnb, nnb, work( ppwo ), 2*nnb,
679 $ q( topq, j ), ldq, work( pw ),
685 CALL sgemm(
'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 slacpy(
'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 slaset(
'All', nblst, nblst, zero, one, work,
705 pw = nblst * nblst + 1
707 CALL slaset(
'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 slaset(
'Lower', ihi - jcol - 1, nnb, zero, zero,
756 $ a( jcol + 2, jcol ), lda )
757 CALL slaset(
'Lower', ihi - jcol - 1, nnb, zero, zero,
758 $ b( jcol + 2, jcol ), ldb )
765 CALL sgemm(
'No Transpose',
'No Transpose', top,
766 $ nblst, nblst, one, a( 1, j ), lda,
767 $ work, nblst, zero, work( pw ), top )
768 CALL slacpy(
'All', top, nblst, work( pw ), top,
770 ppwo = nblst*nblst + 1
772 DO j = j0, jcol+1, -nnb
777 CALL sorm22(
'Right',
'No Transpose', top, 2*nnb,
778 $ nnb, nnb, work( ppwo ), 2*nnb,
779 $ a( 1, j ), lda, work( pw ),
785 CALL sgemm(
'No Transpose',
'No Transpose', top,
786 $ 2*nnb, 2*nnb, one, a( 1, j ), lda,
787 $ work( ppwo ), 2*nnb, zero,
789 CALL slacpy(
'All', top, 2*nnb, work( pw ), top,
792 ppwo = ppwo + 4*nnb*nnb
796 CALL sgemm(
'No Transpose',
'No Transpose', top,
797 $ nblst, nblst, one, b( 1, j ), ldb,
798 $ work, nblst, zero, work( pw ), top )
799 CALL slacpy(
'All', top, nblst, work( pw ), top,
801 ppwo = nblst*nblst + 1
803 DO j = j0, jcol+1, -nnb
808 CALL sorm22(
'Right',
'No Transpose', top, 2*nnb,
809 $ nnb, nnb, work( ppwo ), 2*nnb,
810 $ b( 1, j ), ldb, work( pw ),
816 CALL sgemm(
'No Transpose',
'No Transpose', top,
817 $ 2*nnb, 2*nnb, one, b( 1, j ), ldb,
818 $ work( ppwo ), 2*nnb, zero,
820 CALL slacpy(
'All', top, 2*nnb, work( pw ), top,
823 ppwo = ppwo + 4*nnb*nnb
832 topq = max( 2, j - jcol + 1 )
838 CALL sgemm(
'No Transpose',
'No Transpose', nh,
839 $ nblst, nblst, one, z( topq, j ), ldz,
840 $ work, nblst, zero, work( pw ), nh )
841 CALL slacpy(
'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 sorm22(
'Right',
'No Transpose', nh, 2*nnb,
855 $ nnb, nnb, work( ppwo ), 2*nnb,
856 $ z( topq, j ), ldz, work( pw ),
862 CALL sgemm(
'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 slacpy(
'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 sgghrd( compq2, compz2, n, jcol, ihi, a, lda, b, ldb, q,
889 $ ldq, z, ldz, ierr )
890 work( 1 ) =
REAL( lwkopt )
subroutine slartg(F, G, CS, SN, R)
SLARTG generates a plane rotation with real cosine and real sine.
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
subroutine sgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
SGEMV
subroutine srot(N, SX, INCX, SY, INCY, C, S)
SROT
subroutine sgghd3(COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, Q, LDQ, Z, LDZ, WORK, LWORK, INFO)
SGGHD3
subroutine strmv(UPLO, TRANS, DIAG, N, A, LDA, X, INCX)
STRMV
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 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 sgghrd(COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, Q, LDQ, Z, LDZ, INFO)
SGGHRD
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.