Skip to content

Commit 4bbb9fe

Browse files
martin-frbgmartin-frbg
authored
Fix workspace size (Reference-LAPACK PR 774)
1 parent e486254 commit 4bbb9fe

File tree

2 files changed

+50
-31
lines changed

2 files changed

+50
-31
lines changed

‎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 )

‎lapack-netlib/SRC/ztgsen.f‎

Lines changed: 24 additions & 15 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 ZTGSEN + dependencies
109
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ztgsen.f">
1110
*> [TGZ]</a>
1211
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ztgsen.f">
1312
*> [ZIP]</a>
1413
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ztgsen.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
@@ -290,7 +288,7 @@
290288
*> \author Univ. of Colorado Denver
291289
*> \author NAG Ltd.
292290
*
293-
*> \ingroup complex16OTHERcomputational
291+
*> \ingroup tgsen
294292
*
295293
*> \par Further Details:
296294
* =====================
@@ -427,9 +425,11 @@
427425
*> 1996.
428426
*>
429427
* =====================================================================
430-
SUBROUTINE ZTGSEN( IJOB, WANTQ, WANTZ, SELECT, N, A, LDA, B, LDB,
428+
SUBROUTINE ZTGSEN( 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, --
@@ -468,7 +468,8 @@ SUBROUTINE ZTGSEN( IJOB, WANTQ, WANTZ, SELECT, N, A, LDA, B, LDB,
468468
INTEGER ISAVE( 3 )
469469
* ..
470470
* .. External Subroutines ..
471-
EXTERNAL XERBLA, ZLACN2, ZLACPY, ZLASSQ, ZSCAL, ZTGEXC,
471+
EXTERNAL XERBLA, ZLACN2, ZLACPY, ZLASSQ, ZSCAL,
472+
$ ZTGEXC,
472473
$ ZTGSYL
473474
* ..
474475
* .. Intrinsic Functions ..
@@ -530,7 +531,7 @@ SUBROUTINE ZTGSEN( IJOB, WANTQ, WANTZ, SELECT, N, A, LDA, B, LDB,
530531
END IF
531532
*
532533
IF( IJOB.EQ.1 .OR. IJOB.EQ.2 .OR. IJOB.EQ.4 ) THEN
533-
LWMIN = MAX( 1, 2*M*( N-M ) )
534+
LWMIN = MAX( 1, 2*M*( N-M ) + 1 )
534535
LIWMIN = MAX( 1, N+2 )
535536
ELSE IF( IJOB.EQ.3 .OR. IJOB.EQ.5 ) THEN
536537
LWMIN = MAX( 1, 4*M*( N-M ) )
@@ -592,7 +593,8 @@ SUBROUTINE ZTGSEN( IJOB, WANTQ, WANTZ, SELECT, N, A, LDA, B, LDB,
592593
* and Z that will swap adjacent diagonal blocks in (A, B).
593594
*
594595
IF( K.NE.KS )
595-
$ CALL ZTGEXC( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
596+
$ CALL ZTGEXC( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
597+
$ Z,
596598
$ LDZ, K, KS, IERR )
597599
*
598600
IF( IERR.GT.0 ) THEN
@@ -622,7 +624,8 @@ SUBROUTINE ZTGSEN( IJOB, WANTQ, WANTZ, SELECT, N, A, LDA, B, LDB,
622624
N2 = N - M
623625
I = N1 + 1
624626
CALL ZLACPY( 'Full', N1, N2, A( 1, I ), LDA, WORK, N1 )
625-
CALL ZLACPY( 'Full', N1, N2, B( 1, I ), LDB, WORK( N1*N2+1 ),
627+
CALL ZLACPY( 'Full', N1, N2, B( 1, I ), LDB,
628+
$ WORK( N1*N2+1 ),
626629
$ N1 )
627630
IJB = 0
628631
CALL ZTGSYL( 'N', IJB, N1, N2, A, LDA, A( I, I ), LDA, WORK,
@@ -664,14 +667,16 @@ SUBROUTINE ZTGSEN( IJOB, WANTQ, WANTZ, SELECT, N, A, LDA, B, LDB,
664667
*
665668
* Frobenius norm-based Difu estimate.
666669
*
667-
CALL ZTGSYL( 'N', IJB, N1, N2, A, LDA, A( I, I ), LDA, WORK,
670+
CALL ZTGSYL( 'N', IJB, N1, N2, A, LDA, A( I, I ), LDA,
671+
$ WORK,
668672
$ N1, B, LDB, B( I, I ), LDB, WORK( N1*N2+1 ),
669673
$ N1, DSCALE, DIF( 1 ), WORK( N1*N2*2+1 ),
670674
$ LWORK-2*N1*N2, IWORK, IERR )
671675
*
672676
* Frobenius norm-based Difl estimate.
673677
*
674-
CALL ZTGSYL( 'N', IJB, N2, N1, A( I, I ), LDA, A, LDA, WORK,
678+
CALL ZTGSYL( 'N', IJB, N2, N1, A( I, I ), LDA, A, LDA,
679+
$ WORK,
675680
$ N2, B( I, I ), LDB, B, LDB, WORK( N1*N2+1 ),
676681
$ N2, DSCALE, DIF( 2 ), WORK( N1*N2*2+1 ),
677682
$ LWORK-2*N1*N2, IWORK, IERR )
@@ -699,7 +704,8 @@ SUBROUTINE ZTGSEN( IJOB, WANTQ, WANTZ, SELECT, N, A, LDA, B, LDB,
699704
*
700705
* Solve generalized Sylvester equation
701706
*
702-
CALL ZTGSYL( 'N', IJB, N1, N2, A, LDA, A( I, I ), LDA,
707+
CALL ZTGSYL( 'N', IJB, N1, N2, A, LDA, A( I, I ),
708+
$ LDA,
703709
$ WORK, N1, B, LDB, B( I, I ), LDB,
704710
$ WORK( N1*N2+1 ), N1, DSCALE, DIF( 1 ),
705711
$ WORK( N1*N2*2+1 ), LWORK-2*N1*N2, IWORK,
@@ -708,7 +714,8 @@ SUBROUTINE ZTGSEN( IJOB, WANTQ, WANTZ, SELECT, N, A, LDA, B, LDB,
708714
*
709715
* Solve the transposed variant.
710716
*
711-
CALL ZTGSYL( 'C', IJB, N1, N2, A, LDA, A( I, I ), LDA,
717+
CALL ZTGSYL( 'C', IJB, N1, N2, A, LDA, A( I, I ),
718+
$ LDA,
712719
$ WORK, N1, B, LDB, B( I, I ), LDB,
713720
$ WORK( N1*N2+1 ), N1, DSCALE, DIF( 1 ),
714721
$ WORK( N1*N2*2+1 ), LWORK-2*N1*N2, IWORK,
@@ -728,7 +735,8 @@ SUBROUTINE ZTGSEN( IJOB, WANTQ, WANTZ, SELECT, N, A, LDA, B, LDB,
728735
*
729736
* Solve generalized Sylvester equation
730737
*
731-
CALL ZTGSYL( 'N', IJB, N2, N1, A( I, I ), LDA, A, LDA,
738+
CALL ZTGSYL( 'N', IJB, N2, N1, A( I, I ), LDA, A,
739+
$ LDA,
732740
$ WORK, N2, B( I, I ), LDB, B, LDB,
733741
$ WORK( N1*N2+1 ), N2, DSCALE, DIF( 2 ),
734742
$ WORK( N1*N2*2+1 ), LWORK-2*N1*N2, IWORK,
@@ -737,7 +745,8 @@ SUBROUTINE ZTGSEN( IJOB, WANTQ, WANTZ, SELECT, N, A, LDA, B, LDB,
737745
*
738746
* Solve the transposed variant.
739747
*
740-
CALL ZTGSYL( 'C', IJB, N2, N1, A( I, I ), LDA, A, LDA,
748+
CALL ZTGSYL( 'C', IJB, N2, N1, A( I, I ), LDA, A,
749+
$ LDA,
741750
$ WORK, N2, B, LDB, B( I, I ), LDB,
742751
$ WORK( N1*N2+1 ), N2, DSCALE, DIF( 2 ),
743752
$ WORK( N1*N2*2+1 ), LWORK-2*N1*N2, IWORK,

0 commit comments

Comments
 (0)