Skip to content

Commit d27e98c

Browse files
martin-frbgmartin-frbg
authored
Merge pull request #5734 from martin-frbg/lapack774
Fix workspace size in ?TGSEN (Reference-LAPACK PR 774)
2 parents 429d23f + cc74393 commit d27e98c

File tree

4 files changed

+98
-60
lines changed

4 files changed

+98
-60
lines changed

‎lapack-netlib/SRC/ctgsen.f‎

Lines changed: 23 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -5,15 +5,13 @@
55
* Online html documentation available at
66
* http://www.netlib.org/lapack/explore-html/
77
*
8-
*> \htmlonly
98
*> Download CTGSEN + dependencies
109
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ctgsen.f">
1110
*> [TGZ]</a>
1211
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ctgsen.f">
1312
*> [ZIP]</a>
1413
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ctgsen.f">
1514
*> [TXT]</a>
16-
*> \endhtmlonly
1715
*
1816
* Definition:
1917
* ===========
@@ -242,7 +240,7 @@
242240
*> \verbatim
243241
*> LWORK is INTEGER
244242
*> The dimension of the array WORK. LWORK >= 1
245-
*> If IJOB = 1, 2 or 4, LWORK >= 2*M*(N-M)
243+
*> If IJOB = 1, 2 or 4, LWORK >= 2*M*(N-M) + 1
246244
*> If IJOB = 3 or 5, LWORK >= 4*M*(N-M)
247245
*>
248246
*> If LWORK = -1, then a workspace query is assumed; the routine
@@ -427,9 +425,11 @@
427425
*> 1996.
428426
*>
429427
* =====================================================================
430-
SUBROUTINE CTGSEN( IJOB, WANTQ, WANTZ, SELECT, N, A, LDA, B, LDB,
428+
SUBROUTINE CTGSEN( IJOB, WANTQ, WANTZ, SELECT, N, A, LDA, B,
429+
$ LDB,
431430
$ ALPHA, BETA, Q, LDQ, Z, LDZ, M, PL, PR, DIF,
432431
$ WORK, LWORK, IWORK, LIWORK, INFO )
432+
IMPLICIT NONE
433433
*
434434
* -- LAPACK computational routine --
435435
* -- LAPACK is a software package provided by Univ. of Tennessee, --
@@ -473,7 +473,8 @@ SUBROUTINE CTGSEN( IJOB, WANTQ, WANTZ, SELECT, N, A, LDA, B, LDB,
473473
* ..
474474
* .. External Subroutines ..
475475
REAL SLAMCH
476-
EXTERNAL CLACN2, CLACPY, CLASSQ, CSCAL, CTGEXC, CTGSYL,
476+
EXTERNAL CLACN2, CLACPY, CLASSQ, CSCAL, CTGEXC,
477+
$ CTGSYL,
477478
$ SLAMCH, XERBLA
478479
* ..
479480
* .. Intrinsic Functions ..
@@ -531,7 +532,7 @@ SUBROUTINE CTGSEN( IJOB, WANTQ, WANTZ, SELECT, N, A, LDA, B, LDB,
531532
END IF
532533
*
533534
IF( IJOB.EQ.1 .OR. IJOB.EQ.2 .OR. IJOB.EQ.4 ) THEN
534-
LWMIN = MAX( 1, 2*M*(N-M) )
535+
LWMIN = MAX( 1, 2*M*(N-M) + 1 )
535536
LIWMIN = MAX( 1, N+2 )
536537
ELSE IF( IJOB.EQ.3 .OR. IJOB.EQ.5 ) THEN
537538
LWMIN = MAX( 1, 4*M*(N-M) )
@@ -593,7 +594,8 @@ SUBROUTINE CTGSEN( IJOB, WANTQ, WANTZ, SELECT, N, A, LDA, B, LDB,
593594
* and Z that will swap adjacent diagonal blocks in (A, B).
594595
*
595596
IF( K.NE.KS )
596-
$ CALL CTGEXC( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
597+
$ CALL CTGEXC( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
598+
$ Z,
597599
$ LDZ, K, KS, IERR )
598600
*
599601
IF( IERR.GT.0 ) THEN
@@ -623,7 +625,8 @@ SUBROUTINE CTGSEN( IJOB, WANTQ, WANTZ, SELECT, N, A, LDA, B, LDB,
623625
N2 = N - M
624626
I = N1 + 1
625627
CALL CLACPY( 'Full', N1, N2, A( 1, I ), LDA, WORK, N1 )
626-
CALL CLACPY( 'Full', N1, N2, B( 1, I ), LDB, WORK( N1*N2+1 ),
628+
CALL CLACPY( 'Full', N1, N2, B( 1, I ), LDB,
629+
$ WORK( N1*N2+1 ),
627630
$ N1 )
628631
IJB = 0
629632
CALL CTGSYL( 'N', IJB, N1, N2, A, LDA, A( I, I ), LDA, WORK,
@@ -665,14 +668,16 @@ SUBROUTINE CTGSEN( IJOB, WANTQ, WANTZ, SELECT, N, A, LDA, B, LDB,
665668
*
666669
* Frobenius norm-based Difu estimate.
667670
*
668-
CALL CTGSYL( 'N', IJB, N1, N2, A, LDA, A( I, I ), LDA, WORK,
671+
CALL CTGSYL( 'N', IJB, N1, N2, A, LDA, A( I, I ), LDA,
672+
$ WORK,
669673
$ N1, B, LDB, B( I, I ), LDB, WORK( N1*N2+1 ),
670674
$ N1, DSCALE, DIF( 1 ), WORK( N1*N2*2+1 ),
671675
$ LWORK-2*N1*N2, IWORK, IERR )
672676
*
673677
* Frobenius norm-based Difl estimate.
674678
*
675-
CALL CTGSYL( 'N', IJB, N2, N1, A( I, I ), LDA, A, LDA, WORK,
679+
CALL CTGSYL( 'N', IJB, N2, N1, A( I, I ), LDA, A, LDA,
680+
$ WORK,
676681
$ N2, B( I, I ), LDB, B, LDB, WORK( N1*N2+1 ),
677682
$ N2, DSCALE, DIF( 2 ), WORK( N1*N2*2+1 ),
678683
$ LWORK-2*N1*N2, IWORK, IERR )
@@ -700,7 +705,8 @@ SUBROUTINE CTGSEN( IJOB, WANTQ, WANTZ, SELECT, N, A, LDA, B, LDB,
700705
*
701706
* Solve generalized Sylvester equation
702707
*
703-
CALL CTGSYL( 'N', IJB, N1, N2, A, LDA, A( I, I ), LDA,
708+
CALL CTGSYL( 'N', IJB, N1, N2, A, LDA, A( I, I ),
709+
$ LDA,
704710
$ WORK, N1, B, LDB, B( I, I ), LDB,
705711
$ WORK( N1*N2+1 ), N1, DSCALE, DIF( 1 ),
706712
$ WORK( N1*N2*2+1 ), LWORK-2*N1*N2, IWORK,
@@ -709,7 +715,8 @@ SUBROUTINE CTGSEN( IJOB, WANTQ, WANTZ, SELECT, N, A, LDA, B, LDB,
709715
*
710716
* Solve the transposed variant.
711717
*
712-
CALL CTGSYL( 'C', IJB, N1, N2, A, LDA, A( I, I ), LDA,
718+
CALL CTGSYL( 'C', IJB, N1, N2, A, LDA, A( I, I ),
719+
$ LDA,
713720
$ WORK, N1, B, LDB, B( I, I ), LDB,
714721
$ WORK( N1*N2+1 ), N1, DSCALE, DIF( 1 ),
715722
$ WORK( N1*N2*2+1 ), LWORK-2*N1*N2, IWORK,
@@ -729,7 +736,8 @@ SUBROUTINE CTGSEN( IJOB, WANTQ, WANTZ, SELECT, N, A, LDA, B, LDB,
729736
*
730737
* Solve generalized Sylvester equation
731738
*
732-
CALL CTGSYL( 'N', IJB, N2, N1, A( I, I ), LDA, A, LDA,
739+
CALL CTGSYL( 'N', IJB, N2, N1, A( I, I ), LDA, A,
740+
$ LDA,
733741
$ WORK, N2, B( I, I ), LDB, B, LDB,
734742
$ WORK( N1*N2+1 ), N2, DSCALE, DIF( 2 ),
735743
$ WORK( N1*N2*2+1 ), LWORK-2*N1*N2, IWORK,
@@ -738,7 +746,8 @@ SUBROUTINE CTGSEN( IJOB, WANTQ, WANTZ, SELECT, N, A, LDA, B, LDB,
738746
*
739747
* Solve the transposed variant.
740748
*
741-
CALL CTGSYL( 'C', IJB, N2, N1, A( I, I ), LDA, A, LDA,
749+
CALL CTGSYL( 'C', IJB, N2, N1, A( I, I ), LDA, A,
750+
$ LDA,
742751
$ WORK, N2, B, LDB, B( I, I ), LDB,
743752
$ WORK( N1*N2+1 ), N2, DSCALE, DIF( 2 ),
744753
$ WORK( N1*N2*2+1 ), LWORK-2*N1*N2, IWORK,

‎lapack-netlib/SRC/dtgsen.f‎

Lines changed: 26 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -5,15 +5,13 @@
55
* Online html documentation available at
66
* http://www.netlib.org/lapack/explore-html/
77
*
8-
*> \htmlonly
98
*> Download DTGSEN + dependencies
109
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtgsen.f">
1110
*> [TGZ]</a>
1211
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtgsen.f">
1312
*> [ZIP]</a>
1413
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtgsen.f">
1514
*> [TXT]</a>
16-
*> \endhtmlonly
1715
*
1816
* Definition:
1917
* ===========
@@ -256,7 +254,7 @@
256254
*> \verbatim
257255
*> LWORK is INTEGER
258256
*> The dimension of the array WORK. LWORK >= 4*N+16.
259-
*> If IJOB = 1, 2 or 4, LWORK >= MAX(4*N+16, 2*M*(N-M)).
257+
*> If IJOB = 1, 2 or 4, LWORK >= MAX(4*N+16, 2*M*(N-M) + 1).
260258
*> If IJOB = 3 or 5, LWORK >= MAX(4*N+16, 4*M*(N-M)).
261259
*>
262260
*> If LWORK = -1, then a workspace query is assumed; the routine
@@ -304,7 +302,7 @@
304302
*> \author Univ. of Colorado Denver
305303
*> \author NAG Ltd.
306304
*
307-
*> \ingroup doubleOTHERcomputational
305+
*> \ingroup tgsen
308306
*
309307
*> \par Further Details:
310308
* =====================
@@ -445,9 +443,11 @@
445443
*> \endverbatim
446444
*>
447445
* =====================================================================
448-
SUBROUTINE DTGSEN( IJOB, WANTQ, WANTZ, SELECT, N, A, LDA, B, LDB,
446+
SUBROUTINE DTGSEN( IJOB, WANTQ, WANTZ, SELECT, N, A, LDA, B,
447+
$ LDB,
449448
$ ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDZ, M, PL,
450449
$ PR, DIF, WORK, LWORK, IWORK, LIWORK, INFO )
450+
IMPLICIT NONE
451451
*
452452
* -- LAPACK computational routine --
453453
* -- LAPACK is a software package provided by Univ. of Tennessee, --
@@ -486,7 +486,8 @@ SUBROUTINE DTGSEN( IJOB, WANTQ, WANTZ, SELECT, N, A, LDA, B, LDB,
486486
INTEGER ISAVE( 3 )
487487
* ..
488488
* .. External Subroutines ..
489-
EXTERNAL DLACN2, DLACPY, DLAG2, DLASSQ, DTGEXC, DTGSYL,
489+
EXTERNAL DLACN2, DLACPY, DLAG2, DLASSQ, DTGEXC,
490+
$ DTGSYL,
490491
$ XERBLA
491492
* ..
492493
* .. External Functions ..
@@ -561,7 +562,7 @@ SUBROUTINE DTGSEN( IJOB, WANTQ, WANTZ, SELECT, N, A, LDA, B, LDB,
561562
END IF
562563
*
563564
IF( IJOB.EQ.1 .OR. IJOB.EQ.2 .OR. IJOB.EQ.4 ) THEN
564-
LWMIN = MAX( 1, 4*N+16, 2*M*( N-M ) )
565+
LWMIN = MAX( 1, 4*N+16, 2*M*( N-M ) + 1 )
565566
LIWMIN = MAX( 1, N+6 )
566567
ELSE IF( IJOB.EQ.3 .OR. IJOB.EQ.5 ) THEN
567568
LWMIN = MAX( 1, 4*N+16, 4*M*( N-M ) )
@@ -634,7 +635,8 @@ SUBROUTINE DTGSEN( IJOB, WANTQ, WANTZ, SELECT, N, A, LDA, B, LDB,
634635
*
635636
KK = K
636637
IF( K.NE.KS )
637-
$ CALL DTGEXC( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
638+
$ CALL DTGEXC( WANTQ, WANTZ, N, A, LDA, B, LDB, Q,
639+
$ LDQ,
638640
$ Z, LDZ, KK, KS, WORK, LWORK, IERR )
639641
*
640642
IF( IERR.GT.0 ) THEN
@@ -668,7 +670,8 @@ SUBROUTINE DTGSEN( IJOB, WANTQ, WANTZ, SELECT, N, A, LDA, B, LDB,
668670
I = N1 + 1
669671
IJB = 0
670672
CALL DLACPY( 'Full', N1, N2, A( 1, I ), LDA, WORK, N1 )
671-
CALL DLACPY( 'Full', N1, N2, B( 1, I ), LDB, WORK( N1*N2+1 ),
673+
CALL DLACPY( 'Full', N1, N2, B( 1, I ), LDB,
674+
$ WORK( N1*N2+1 ),
672675
$ N1 )
673676
CALL DTGSYL( 'N', IJB, N1, N2, A, LDA, A( I, I ), LDA, WORK,
674677
$ N1, B, LDB, B( I, I ), LDB, WORK( N1*N2+1 ), N1,
@@ -710,14 +713,16 @@ SUBROUTINE DTGSEN( IJOB, WANTQ, WANTZ, SELECT, N, A, LDA, B, LDB,
710713
*
711714
* Frobenius norm-based Difu-estimate.
712715
*
713-
CALL DTGSYL( 'N', IJB, N1, N2, A, LDA, A( I, I ), LDA, WORK,
716+
CALL DTGSYL( 'N', IJB, N1, N2, A, LDA, A( I, I ), LDA,
717+
$ WORK,
714718
$ N1, B, LDB, B( I, I ), LDB, WORK( N1*N2+1 ),
715719
$ N1, DSCALE, DIF( 1 ), WORK( 2*N1*N2+1 ),
716720
$ LWORK-2*N1*N2, IWORK, IERR )
717721
*
718722
* Frobenius norm-based Difl-estimate.
719723
*
720-
CALL DTGSYL( 'N', IJB, N2, N1, A( I, I ), LDA, A, LDA, WORK,
724+
CALL DTGSYL( 'N', IJB, N2, N1, A( I, I ), LDA, A, LDA,
725+
$ WORK,
721726
$ N2, B( I, I ), LDB, B, LDB, WORK( N1*N2+1 ),
722727
$ N2, DSCALE, DIF( 2 ), WORK( 2*N1*N2+1 ),
723728
$ LWORK-2*N1*N2, IWORK, IERR )
@@ -746,7 +751,8 @@ SUBROUTINE DTGSEN( IJOB, WANTQ, WANTZ, SELECT, N, A, LDA, B, LDB,
746751
*
747752
* Solve generalized Sylvester equation.
748753
*
749-
CALL DTGSYL( 'N', IJB, N1, N2, A, LDA, A( I, I ), LDA,
754+
CALL DTGSYL( 'N', IJB, N1, N2, A, LDA, A( I, I ),
755+
$ LDA,
750756
$ WORK, N1, B, LDB, B( I, I ), LDB,
751757
$ WORK( N1*N2+1 ), N1, DSCALE, DIF( 1 ),
752758
$ WORK( 2*N1*N2+1 ), LWORK-2*N1*N2, IWORK,
@@ -755,7 +761,8 @@ SUBROUTINE DTGSEN( IJOB, WANTQ, WANTZ, SELECT, N, A, LDA, B, LDB,
755761
*
756762
* Solve the transposed variant.
757763
*
758-
CALL DTGSYL( 'T', IJB, N1, N2, A, LDA, A( I, I ), LDA,
764+
CALL DTGSYL( 'T', IJB, N1, N2, A, LDA, A( I, I ),
765+
$ LDA,
759766
$ WORK, N1, B, LDB, B( I, I ), LDB,
760767
$ WORK( N1*N2+1 ), N1, DSCALE, DIF( 1 ),
761768
$ WORK( 2*N1*N2+1 ), LWORK-2*N1*N2, IWORK,
@@ -775,7 +782,8 @@ SUBROUTINE DTGSEN( IJOB, WANTQ, WANTZ, SELECT, N, A, LDA, B, LDB,
775782
*
776783
* Solve generalized Sylvester equation.
777784
*
778-
CALL DTGSYL( 'N', IJB, N2, N1, A( I, I ), LDA, A, LDA,
785+
CALL DTGSYL( 'N', IJB, N2, N1, A( I, I ), LDA, A,
786+
$ LDA,
779787
$ WORK, N2, B( I, I ), LDB, B, LDB,
780788
$ WORK( N1*N2+1 ), N2, DSCALE, DIF( 2 ),
781789
$ WORK( 2*N1*N2+1 ), LWORK-2*N1*N2, IWORK,
@@ -784,7 +792,8 @@ SUBROUTINE DTGSEN( IJOB, WANTQ, WANTZ, SELECT, N, A, LDA, B, LDB,
784792
*
785793
* Solve the transposed variant.
786794
*
787-
CALL DTGSYL( 'T', IJB, N2, N1, A( I, I ), LDA, A, LDA,
795+
CALL DTGSYL( 'T', IJB, N2, N1, A( I, I ), LDA, A,
796+
$ LDA,
788797
$ WORK, N2, B( I, I ), LDB, B, LDB,
789798
$ WORK( N1*N2+1 ), N2, DSCALE, DIF( 2 ),
790799
$ WORK( 2*N1*N2+1 ), LWORK-2*N1*N2, IWORK,
@@ -826,7 +835,8 @@ SUBROUTINE DTGSEN( IJOB, WANTQ, WANTZ, SELECT, N, A, LDA, B, LDB,
826835
WORK( 6 ) = B( K+1, K )
827836
WORK( 7 ) = B( K, K+1 )
828837
WORK( 8 ) = B( K+1, K+1 )
829-
CALL DLAG2( WORK, 2, WORK( 5 ), 2, SMLNUM*EPS, BETA( K ),
838+
CALL DLAG2( WORK, 2, WORK( 5 ), 2, SMLNUM*EPS,
839+
$ BETA( K ),
830840
$ BETA( K+1 ), ALPHAR( K ), ALPHAR( K+1 ),
831841
$ ALPHAI( K ) )
832842
ALPHAI( K+1 ) = -ALPHAI( K )

0 commit comments

Comments
 (0)