Skip to content

Commit fdcae32

Browse files
christoph-conradschristoph-conrads
committed
xGGQRCS: fix accidental memory allocation
1 parent dfe2bdd commit fdcae32

File tree

2 files changed

+74
-70
lines changed

2 files changed

+74
-70
lines changed

‎SRC/dggqrcs.f

+37-35
Original file line numberDiff line numberDiff line change
@@ -350,11 +350,10 @@ RECURSIVE SUBROUTINE DGGQRCS( JOBU1, JOBU2, JOBX, M, N, P, L,
350350
*
351351
* .. Local Scalars ..
352352
LOGICAL WANTU1, WANTU2, WANTX, LQUERY
353-
INTEGER I, J, K, K1, LMAX, Z, LDG, LDX, LDVT, LWKOPT
353+
INTEGER I, J, K, K1, LMAX, IG, IG11, IG21, IG22,
354+
$ IVT, IVT12, LDG, LDX, LDVT, LWKOPT
354355
DOUBLE PRECISION BASE, NAN, NORMA, NORMB, NORMG, TOL, ULP, UNFL,
355356
$ THETA, IOTA, W
356-
* .. Local Arrays ..
357-
DOUBLE PRECISION G( M + P, N ), VT( N, N )
358357
* ..
359358
* .. External Functions ..
360359
LOGICAL LSAME
@@ -438,12 +437,15 @@ RECURSIVE SUBROUTINE DGGQRCS( JOBU1, JOBU2, JOBX, M, N, P, L,
438437
*
439438
SWAPPED = .FALSE.
440439
L = 0
441-
LMAX = MIN( M + P, N )
442-
Z = ( M + P ) * N
443-
G = WORK( 1 )
444440
LDG = M + P
445-
VT = 0
446441
LDVT = N
442+
LMAX = MIN( M + P, N )
443+
IG = 1
444+
IG11 = 1
445+
IG21 = M + 1
446+
IG22 = LDG * M + M + 1
447+
IVT = LDG * N + 2
448+
IVT12 = IVT + LDVT * M
447449
THETA = NAN
448450
IOTA = NAN
449451
W = NAN
@@ -453,21 +455,23 @@ RECURSIVE SUBROUTINE DGGQRCS( JOBU1, JOBU2, JOBX, M, N, P, L,
453455
IF( INFO.EQ.0 ) THEN
454456
LWKOPT = 0
455457
*
456-
CALL DGEQP3( M+P, N, G, LDG, IWORK, ALPHA, WORK, -1, INFO )
458+
CALL DGEQP3( M + P, N, WORK( IG ), LDG, IWORK, ALPHA, WORK, -1,
459+
$ INFO )
457460
LWKOPT = MAX( LWKOPT, INT( WORK( 1 ) ) )
458461
LWKOPT = INT( WORK( 1 ) )
459462
*
460-
CALL DORGQR( M + P, LMAX, LMAX, G, LDG, ALPHA, WORK, -1, INFO )
463+
CALL DORGQR( M + P, LMAX, LMAX, WORK( IG ), LDG, ALPHA, WORK,
464+
$ -1, INFO )
461465
LWKOPT = MAX( LWKOPT, INT( WORK( 1 ) ) )
462466
*
463467
CALL DORCSD2BY1( JOBU1, JOBU2, JOBX, M + P, M, LMAX,
464-
$ G, LDG, G, LDG,
468+
$ WORK( IG ), LDG, WORK( IG ), LDG,
465469
$ ALPHA,
466-
$ U1, LDU1, U2, LDU2, VT, LDVT,
470+
$ U1, LDU1, U2, LDU2, WORK( IVT ), LDVT,
467471
$ WORK, -1, IWORK, INFO )
468472
LWKOPT = MAX( LWKOPT, INT( WORK( 1 ) ) )
469473
* The matrix (A, B) must be stored sequentially for DORGQR
470-
LWKOPT = LWKOPT + Z
474+
LWKOPT = LWKOPT + IVT
471475
* 2-by-1 CSD matrix V1 must be stored
472476
IF( WANTX ) THEN
473477
LWKOPT = LWKOPT + LDVT*N
@@ -484,9 +488,7 @@ RECURSIVE SUBROUTINE DGGQRCS( JOBU1, JOBU2, JOBX, M, N, P, L,
484488
RETURN
485489
ENDIF
486490
* Finish initialization
487-
IF( WANTX ) THEN
488-
VT = WORK( Z + 1 )
489-
ELSE
491+
IF( .NOT.WANTX ) THEN
490492
LDVT = 0
491493
END IF
492494
*
@@ -506,8 +508,8 @@ RECURSIVE SUBROUTINE DGGQRCS( JOBU1, JOBU2, JOBX, M, N, P, L,
506508
*
507509
* Copy matrices A, B into the (M+P) x N matrix G
508510
*
509-
CALL DLACPY( 'A', M, N, A, LDA, G( 1, 1 ), LDG )
510-
CALL DLACPY( 'A', P, N, B, LDB, G( M + 1, 1 ), LDG )
511+
CALL DLACPY( 'A', M, N, A, LDA, WORK( IG11 ), LDG )
512+
CALL DLACPY( 'A', P, N, B, LDB, WORK( IG21 ), LDG )
511513
*
512514
* DEBUG
513515
*
@@ -533,16 +535,16 @@ RECURSIVE SUBROUTINE DGGQRCS( JOBU1, JOBU2, JOBX, M, N, P, L,
533535
*
534536
* Compute the QR factorization with column pivoting GΠ = Q1 R1
535537
*
536-
CALL DGEQP3( M + P, N, G, LDG, IWORK, ALPHA,
537-
$ WORK( Z + 1 ), LWORK - Z, INFO )
538+
CALL DGEQP3( M + P, N, WORK( IG ), LDG, IWORK, ALPHA,
539+
$ WORK( IVT ), LWORK - IVT + 1, INFO )
538540
IF( INFO.NE.0 ) THEN
539541
RETURN
540542
END IF
541543
*
542544
* Determine the rank of G
543545
*
544546
DO I = 1, MIN( M + P, N )
545-
IF( ABS( G( I, I ) ).LE.TOL ) THEN
547+
IF( ABS( WORK( (I-1) * LDG + I ) ).LE.TOL ) THEN
546548
EXIT
547549
END IF
548550
L = L + 1
@@ -566,11 +568,11 @@ RECURSIVE SUBROUTINE DGGQRCS( JOBU1, JOBU2, JOBX, M, N, P, L,
566568
*
567569
IF( WANTX ) THEN
568570
IF( L.LE.M ) THEN
569-
CALL DLACPY( 'U', L, N, G, LDG, A, LDA )
571+
CALL DLACPY( 'U', L, N, WORK( IG ), LDG, A, LDA )
570572
CALL DLASET( 'L', L - 1, N, 0.0D0, 0.0D0, A( 2, 1 ), LDA )
571573
ELSE
572-
CALL DLACPY( 'U', M, N, G, LDG, A, LDA )
573-
CALL DLACPY( 'U', L - M, N - M, G( M+1,M+1 ), LDG, B, LDB )
574+
CALL DLACPY( 'U', M, N, WORK( IG ), LDG, A, LDA )
575+
CALL DLACPY( 'U', L - M, N - M, WORK( IG22 ), LDG, B, LDB )
574576
*
575577
CALL DLASET( 'L', M - 1, N, 0.0D0, 0.0D0, A( 2, 1 ), LDA )
576578
CALL DLASET( 'L', L-M-1, N, 0.0D0, 0.0D0, B( 2, 1 ), LDB )
@@ -579,8 +581,8 @@ RECURSIVE SUBROUTINE DGGQRCS( JOBU1, JOBU2, JOBX, M, N, P, L,
579581
*
580582
* Explicitly form Q1 so that we can compute the CS decomposition
581583
*
582-
CALL DORGQR( M + P, L, L, G, LDG, ALPHA,
583-
$ WORK( Z + 1 ), LWORK - Z, INFO )
584+
CALL DORGQR( M + P, L, L, WORK( IG ), LDG, ALPHA,
585+
$ WORK( IVT ), LWORK - IVT + 1, INFO )
584586
IF ( INFO.NE.0 ) THEN
585587
RETURN
586588
END IF
@@ -595,10 +597,10 @@ RECURSIVE SUBROUTINE DGGQRCS( JOBU1, JOBU2, JOBX, M, N, P, L,
595597
K = MIN( M, P, L, M + P - L )
596598
K1 = MAX( L - P, 0 )
597599
CALL DORCSD2BY1( JOBU1, JOBU2, JOBX, M + P, M, L,
598-
$ G( 1, 1 ), LDG, G( M + 1, 1 ), LDG,
600+
$ WORK( IG11 ), LDG, WORK( IG21 ), LDG,
599601
$ ALPHA,
600-
$ U1, LDU1, U2, LDU2, VT, LDVT,
601-
$ WORK( Z + LDVT*N + 1 ), LWORK - Z - LDVT*N,
602+
$ U1, LDU1, U2, LDU2, WORK( IVT ), LDVT,
603+
$ WORK( IVT + LDVT*N ), LWORK - IVT - LDVT*N + 1,
602604
$ IWORK( N + 1 ), INFO )
603605
IF( INFO.NE.0 ) THEN
604606
RETURN
@@ -614,14 +616,14 @@ RECURSIVE SUBROUTINE DGGQRCS( JOBU1, JOBU2, JOBX, M, N, P, L,
614616
LDX = L
615617
IF ( L.LE.M ) THEN
616618
CALL DGEMM( 'N', 'N', L, N, L,
617-
$ 1.0D0, VT, LDVT, A, LDA,
619+
$ 1.0D0, WORK( IVT ), LDVT, A, LDA,
618620
$ 0.0D0, WORK( 2 ), LDX )
619621
ELSE
620622
CALL DGEMM( 'N', 'N', L, N, M,
621-
$ 1.0D0, VT( 1, 1 ), LDVT, A, LDA,
623+
$ 1.0D0, WORK( IVT ), LDVT, A, LDA,
622624
$ 0.0D0, WORK( 2 ), LDX )
623625
CALL DGEMM( 'N', 'N', L, N - M, L - M,
624-
$ 1.0D0, VT( 1, M + 1 ), LDVT, B, LDB,
626+
$ 1.0D0, WORK( IVT12 ), LDVT, B, LDB,
625627
$ 1.0D0, WORK( L*M + 2 ), LDX )
626628
END IF
627629
* Revert column permutation Π by permuting the columns of X
@@ -640,7 +642,7 @@ RECURSIVE SUBROUTINE DGGQRCS( JOBU1, JOBU2, JOBX, M, N, P, L,
640642
ALPHA( I ) = 0.0D0
641643
BETA( I ) = 1.0D0
642644
IF( WANTX ) THEN
643-
WORK( Z + I + 1 ) = 1.0D0
645+
WORK( IVT + I ) = 1.0D0
644646
END IF
645647
ELSE
646648
* iota comes in the greek alphabet after theta
@@ -651,13 +653,13 @@ RECURSIVE SUBROUTINE DGGQRCS( JOBU1, JOBU2, JOBX, M, N, P, L,
651653
ALPHA( I ) = ( SIN( IOTA ) / TAN( THETA ) ) / W
652654
BETA( I ) = SIN( IOTA )
653655
IF( WANTX ) THEN
654-
WORK( Z + I + 1 ) = SIN( THETA ) / SIN( IOTA )
656+
WORK( IVT + I ) = SIN( THETA ) / SIN( IOTA )
655657
END IF
656658
ELSE
657659
ALPHA( I ) = COS( IOTA )
658660
BETA( I ) = SIN( IOTA )
659661
IF( WANTX ) THEN
660-
WORK( Z + I + 1 ) = COS( THETA ) / COS( IOTA ) / W
662+
WORK( IVT + I ) = COS( THETA ) / COS( IOTA ) / W
661663
END IF
662664
END IF
663665
END IF
@@ -670,7 +672,7 @@ RECURSIVE SUBROUTINE DGGQRCS( JOBU1, JOBU2, JOBX, M, N, P, L,
670672
END DO
671673
DO I = 1, K
672674
WORK( LDX*J + I + K1 + 1 ) =
673-
$ WORK( LDX*J + I + K1 + 1 ) * WORK( Z + I + 1 )
675+
$ WORK( LDX*J + I + K1 + 1 ) * WORK( IVT + I )
674676
END DO
675677
END DO
676678
END IF

‎SRC/sggqrcs.f

+37-35
Original file line numberDiff line numberDiff line change
@@ -350,11 +350,10 @@ RECURSIVE SUBROUTINE SGGQRCS( JOBU1, JOBU2, JOBX, M, N, P, L,
350350
*
351351
* .. Local Scalars ..
352352
LOGICAL WANTU1, WANTU2, WANTX, LQUERY
353-
INTEGER I, J, K, K1, LMAX, Z, LDG, LDX, LDVT, LWKOPT
353+
INTEGER I, J, K, K1, LMAX, IG, IG11, IG21, IG22,
354+
$ IVT, IVT12, LDG, LDX, LDVT, LWKOPT
354355
REAL BASE, NAN, NORMA, NORMB, NORMG, TOL, ULP, UNFL,
355356
$ THETA, IOTA, W
356-
* .. Local Arrays ..
357-
REAL G( M + P, N ), VT( N, N )
358357
* ..
359358
* .. External Functions ..
360359
LOGICAL LSAME
@@ -438,12 +437,15 @@ RECURSIVE SUBROUTINE SGGQRCS( JOBU1, JOBU2, JOBX, M, N, P, L,
438437
*
439438
SWAPPED = .FALSE.
440439
L = 0
441-
LMAX = MIN( M + P, N )
442-
Z = ( M + P ) * N
443-
G = WORK( 1 )
444440
LDG = M + P
445-
VT = 0
446441
LDVT = N
442+
LMAX = MIN( M + P, N )
443+
IG = 1
444+
IG11 = 1
445+
IG21 = M + 1
446+
IG22 = LDG * M + M + 1
447+
IVT = LDG * N + 2
448+
IVT12 = IVT + LDVT * M
447449
THETA = NAN
448450
IOTA = NAN
449451
W = NAN
@@ -453,21 +455,23 @@ RECURSIVE SUBROUTINE SGGQRCS( JOBU1, JOBU2, JOBX, M, N, P, L,
453455
IF( INFO.EQ.0 ) THEN
454456
LWKOPT = 0
455457
*
456-
CALL SGEQP3( M+P, N, G, LDG, IWORK, ALPHA, WORK, -1, INFO )
458+
CALL SGEQP3( M + P, N, WORK( IG ), LDG, IWORK, ALPHA, WORK, -1,
459+
$ INFO )
457460
LWKOPT = MAX( LWKOPT, INT( WORK( 1 ) ) )
458461
LWKOPT = INT( WORK( 1 ) )
459462
*
460-
CALL SORGQR( M + P, LMAX, LMAX, G, LDG, ALPHA, WORK, -1, INFO )
463+
CALL SORGQR( M + P, LMAX, LMAX, WORK( IG ), LDG, ALPHA, WORK,
464+
$ -1, INFO )
461465
LWKOPT = MAX( LWKOPT, INT( WORK( 1 ) ) )
462466
*
463467
CALL SORCSD2BY1( JOBU1, JOBU2, JOBX, M + P, M, LMAX,
464-
$ G, LDG, G, LDG,
468+
$ WORK( IG ), LDG, WORK( IG ), LDG,
465469
$ ALPHA,
466-
$ U1, LDU1, U2, LDU2, VT, LDVT,
470+
$ U1, LDU1, U2, LDU2, WORK( IVT ), LDVT,
467471
$ WORK, -1, IWORK, INFO )
468472
LWKOPT = MAX( LWKOPT, INT( WORK( 1 ) ) )
469473
* The matrix (A, B) must be stored sequentially for SORGQR
470-
LWKOPT = LWKOPT + Z
474+
LWKOPT = LWKOPT + IVT
471475
* 2-by-1 CSD matrix V1 must be stored
472476
IF( WANTX ) THEN
473477
LWKOPT = LWKOPT + LDVT*N
@@ -484,9 +488,7 @@ RECURSIVE SUBROUTINE SGGQRCS( JOBU1, JOBU2, JOBX, M, N, P, L,
484488
RETURN
485489
ENDIF
486490
* Finish initialization
487-
IF( WANTX ) THEN
488-
VT = WORK( Z + 1 )
489-
ELSE
491+
IF( .NOT.WANTX ) THEN
490492
LDVT = 0
491493
END IF
492494
*
@@ -506,8 +508,8 @@ RECURSIVE SUBROUTINE SGGQRCS( JOBU1, JOBU2, JOBX, M, N, P, L,
506508
*
507509
* Copy matrices A, B into the (M+P) x N matrix G
508510
*
509-
CALL SLACPY( 'A', M, N, A, LDA, G( 1, 1 ), LDG )
510-
CALL SLACPY( 'A', P, N, B, LDB, G( M + 1, 1 ), LDG )
511+
CALL SLACPY( 'A', M, N, A, LDA, WORK( IG11 ), LDG )
512+
CALL SLACPY( 'A', P, N, B, LDB, WORK( IG21 ), LDG )
511513
*
512514
* DEBUG
513515
*
@@ -533,16 +535,16 @@ RECURSIVE SUBROUTINE SGGQRCS( JOBU1, JOBU2, JOBX, M, N, P, L,
533535
*
534536
* Compute the QR factorization with column pivoting GΠ = Q1 R1
535537
*
536-
CALL SGEQP3( M + P, N, G, LDG, IWORK, ALPHA,
537-
$ WORK( Z + 1 ), LWORK - Z, INFO )
538+
CALL SGEQP3( M + P, N, WORK( IG ), LDG, IWORK, ALPHA,
539+
$ WORK( IVT ), LWORK - IVT + 1, INFO )
538540
IF( INFO.NE.0 ) THEN
539541
RETURN
540542
END IF
541543
*
542544
* Determine the rank of G
543545
*
544546
DO I = 1, MIN( M + P, N )
545-
IF( ABS( G( I, I ) ).LE.TOL ) THEN
547+
IF( ABS( WORK( (I-1) * LDG + I ) ).LE.TOL ) THEN
546548
EXIT
547549
END IF
548550
L = L + 1
@@ -566,11 +568,11 @@ RECURSIVE SUBROUTINE SGGQRCS( JOBU1, JOBU2, JOBX, M, N, P, L,
566568
*
567569
IF( WANTX ) THEN
568570
IF( L.LE.M ) THEN
569-
CALL SLACPY( 'U', L, N, G, LDG, A, LDA )
571+
CALL SLACPY( 'U', L, N, WORK( IG ), LDG, A, LDA )
570572
CALL SLASET( 'L', L - 1, N, 0.0E0, 0.0E0, A( 2, 1 ), LDA )
571573
ELSE
572-
CALL SLACPY( 'U', M, N, G, LDG, A, LDA )
573-
CALL SLACPY( 'U', L - M, N - M, G( M+1,M+1 ), LDG, B, LDB )
574+
CALL SLACPY( 'U', M, N, WORK( IG ), LDG, A, LDA )
575+
CALL SLACPY( 'U', L - M, N - M, WORK( IG22 ), LDG, B, LDB )
574576
*
575577
CALL SLASET( 'L', M - 1, N, 0.0E0, 0.0E0, A( 2, 1 ), LDA )
576578
CALL SLASET( 'L', L-M-1, N, 0.0E0, 0.0E0, B( 2, 1 ), LDB )
@@ -579,8 +581,8 @@ RECURSIVE SUBROUTINE SGGQRCS( JOBU1, JOBU2, JOBX, M, N, P, L,
579581
*
580582
* Explicitly form Q1 so that we can compute the CS decomposition
581583
*
582-
CALL SORGQR( M + P, L, L, G, LDG, ALPHA,
583-
$ WORK( Z + 1 ), LWORK - Z, INFO )
584+
CALL SORGQR( M + P, L, L, WORK( IG ), LDG, ALPHA,
585+
$ WORK( IVT ), LWORK - IVT + 1, INFO )
584586
IF ( INFO.NE.0 ) THEN
585587
RETURN
586588
END IF
@@ -595,10 +597,10 @@ RECURSIVE SUBROUTINE SGGQRCS( JOBU1, JOBU2, JOBX, M, N, P, L,
595597
K = MIN( M, P, L, M + P - L )
596598
K1 = MAX( L - P, 0 )
597599
CALL SORCSD2BY1( JOBU1, JOBU2, JOBX, M + P, M, L,
598-
$ G( 1, 1 ), LDG, G( M + 1, 1 ), LDG,
600+
$ WORK( IG11 ), LDG, WORK( IG21 ), LDG,
599601
$ ALPHA,
600-
$ U1, LDU1, U2, LDU2, VT, LDVT,
601-
$ WORK( Z + LDVT*N + 1 ), LWORK - Z - LDVT*N,
602+
$ U1, LDU1, U2, LDU2, WORK( IVT ), LDVT,
603+
$ WORK( IVT + LDVT*N ), LWORK - IVT - LDVT*N + 1,
602604
$ IWORK( N + 1 ), INFO )
603605
IF( INFO.NE.0 ) THEN
604606
RETURN
@@ -614,14 +616,14 @@ RECURSIVE SUBROUTINE SGGQRCS( JOBU1, JOBU2, JOBX, M, N, P, L,
614616
LDX = L
615617
IF ( L.LE.M ) THEN
616618
CALL SGEMM( 'N', 'N', L, N, L,
617-
$ 1.0E0, VT, LDVT, A, LDA,
619+
$ 1.0E0, WORK( IVT ), LDVT, A, LDA,
618620
$ 0.0E0, WORK( 2 ), LDX )
619621
ELSE
620622
CALL SGEMM( 'N', 'N', L, N, M,
621-
$ 1.0E0, VT( 1, 1 ), LDVT, A, LDA,
623+
$ 1.0E0, WORK( IVT ), LDVT, A, LDA,
622624
$ 0.0E0, WORK( 2 ), LDX )
623625
CALL SGEMM( 'N', 'N', L, N - M, L - M,
624-
$ 1.0E0, VT( 1, M + 1 ), LDVT, B, LDB,
626+
$ 1.0E0, WORK( IVT12 ), LDVT, B, LDB,
625627
$ 1.0E0, WORK( L*M + 2 ), LDX )
626628
END IF
627629
* Revert column permutation Π by permuting the columns of X
@@ -640,7 +642,7 @@ RECURSIVE SUBROUTINE SGGQRCS( JOBU1, JOBU2, JOBX, M, N, P, L,
640642
ALPHA( I ) = 0.0E0
641643
BETA( I ) = 1.0E0
642644
IF( WANTX ) THEN
643-
WORK( Z + I + 1 ) = 1.0E0
645+
WORK( IVT + I ) = 1.0E0
644646
END IF
645647
ELSE
646648
* iota comes in the greek alphabet after theta
@@ -651,13 +653,13 @@ RECURSIVE SUBROUTINE SGGQRCS( JOBU1, JOBU2, JOBX, M, N, P, L,
651653
ALPHA( I ) = ( SIN( IOTA ) / TAN( THETA ) ) / W
652654
BETA( I ) = SIN( IOTA )
653655
IF( WANTX ) THEN
654-
WORK( Z + I + 1 ) = SIN( THETA ) / SIN( IOTA )
656+
WORK( IVT + I ) = SIN( THETA ) / SIN( IOTA )
655657
END IF
656658
ELSE
657659
ALPHA( I ) = COS( IOTA )
658660
BETA( I ) = SIN( IOTA )
659661
IF( WANTX ) THEN
660-
WORK( Z + I + 1 ) = COS( THETA ) / COS( IOTA ) / W
662+
WORK( IVT + I ) = COS( THETA ) / COS( IOTA ) / W
661663
END IF
662664
END IF
663665
END IF
@@ -670,7 +672,7 @@ RECURSIVE SUBROUTINE SGGQRCS( JOBU1, JOBU2, JOBX, M, N, P, L,
670672
END DO
671673
DO I = 1, K
672674
WORK( LDX*J + I + K1 + 1 ) =
673-
$ WORK( LDX*J + I + K1 + 1 ) * WORK( Z + I + 1 )
675+
$ WORK( LDX*J + I + K1 + 1 ) * WORK( IVT + I )
674676
END DO
675677
END DO
676678
END IF

0 commit comments

Comments
 (0)