Patch provided by Lawrence Mulholland from Nag on Nov 1st 2013
authorjulie <julielangou@users.noreply.github.com>
Sun, 17 Nov 2013 00:01:25 +0000 (00:01 +0000)
committerjulie <julielangou@users.noreply.github.com>
Sun, 17 Nov 2013 00:01:25 +0000 (00:01 +0000)
Email below:
============
I have been incorporating some routines into the NAG Library,
which means some automatic code translation and writing some
example and test programs.

The routines I have been adding are:
?geqrt, ?gemqrt, ?tpqrt, ?tpmqrt, ?orcsd, ?uncsd

At the end of this message I will give you my current svn status
and svn diff for consideration and approval before I commit.

In each case, when testing immediate exits, my tests failed because
constraints were mutually exclusive for the immediate return case.
I have already committed changes to the constraints for some of
the above to allow immediate exit.
I have completed this for the remainder of this set.

Less importantly, there are things in the code that trip up a checking
compiler:
  a) an
           IF ( clause1(i) .AND. clause2(array(i)) ) THEN

      where array(i) is either not initialized or is out of bounds if
      clause1(i) is .FALSE.

      This is wrong since a Fortran compiler is at liberty to test clause2 first.
      In my changes this has been split into two as best suits the case.

 b) an
          CALL SUB (i, array(N-i+2))
     with i = 1 and array(N+1) either not initialized or out of bounds, but
     internally array(N+1) is not referenced.

     In this case I don't think the Fortran standard is clear, but it trips up the
     nagfor compiler with checking on. So in the NAG incorporated versions
     of Lapack routines such calls are protected and/or
     a special i=1 call is made.
     The changes I want to commit also do this.

 c) workspace queries passing zero instead of array references
         e.g.
                lwork = -1
                call barf(n,m,0,0,0,0,0,-1,info)

     a checking compiler won't like this.
     I have changed cases like this to pass available arrays of sufficient size
     and the right shape in place of the zeros.

25 files changed:
SRC/cbbcsd.f
SRC/cgemqrt.f
SRC/clahef_rook.f
SRC/ctpmqrt.f
SRC/ctpqrt.f
SRC/cunbdb.f
SRC/cuncsd.f
SRC/dbbcsd.f
SRC/dgemqrt.f
SRC/dorbdb.f
SRC/dorcsd.f
SRC/dtpmqrt.f
SRC/dtpqrt.f
SRC/sbbcsd.f
SRC/sgemqrt.f
SRC/sorbdb.f
SRC/sorcsd.f
SRC/stpmqrt.f
SRC/stpqrt.f
SRC/zbbcsd.f
SRC/zgemqrt.f
SRC/ztpmqrt.f
SRC/ztpqrt.f
SRC/zunbdb.f
SRC/zuncsd.f

index fdc2bab..9085821 100644 (file)
 *     Initial deflation
 *
       IMAX = Q
-      DO WHILE((IMAX.GT.1).AND.(PHI(MAX(IMAX-1,1)).EQ.ZERO))
+      DO WHILE( IMAX .GT. 1 )
+         IF( PHI(IMAX-1) .NE. ZERO ) THEN
+            EXIT
+         END IF
          IMAX = IMAX - 1
       END DO
       IMIN = IMAX - 1
index 76f345c..d24e88c 100644 (file)
          INFO = -4
       ELSE IF( K.LT.0 .OR. K.GT.Q ) THEN
          INFO = -5
-      ELSE IF( NB.LT.1 .OR. NB.GT.K ) THEN
+      ELSE IF( NB.LT.1 .OR. (NB.GT.K .AND. K.GT.0)) THEN
          INFO = -6
       ELSE IF( LDV.LT.MAX( 1, Q ) ) THEN
          INFO = -8
index d189376..e5b5a03 100644 (file)
       EXTERNAL           CCOPY, CSSCAL, CGEMM, CGEMV, CLACGV, CSWAP
 *     ..
 *     .. Intrinsic Functions ..
-      INTRINSIC          ABS, CONJG, IMAG, MAX, MIN, REAL, SQRT
+      INTRINSIC          ABS, CONJG, AIMAG, MAX, MIN, REAL, SQRT
 *     ..
 *     .. Statement Functions ..
       REAL               CABS1
 *     ..
 *     .. Statement Function definitions ..
-      CABS1( Z ) = ABS( REAL( Z ) ) + ABS( IMAG( Z ) )
+      CABS1( Z ) = ABS( REAL( Z ) ) + ABS( AIMAG( Z ) )
 *     ..
 *     .. Executable Statements ..
 *
index abfdae8..7c767f8 100644 (file)
 *     ..
 *     .. Local Scalars ..
       LOGICAL            LEFT, RIGHT, TRAN, NOTRAN
-      INTEGER            I, IB, MB, LB, KF, Q
+      INTEGER            I, IB, MB, LB, KF, LDAQ, LDVQ
 *     ..
 *     .. External Functions ..
       LOGICAL            LSAME
       TRAN   = LSAME( TRANS, 'C' )
       NOTRAN = LSAME( TRANS, 'N' )
 *      
-      IF( LEFT ) THEN
-         Q = M
+      IF ( LEFT ) THEN
+         LDVQ = MAX( 1, M )
+         LDAQ = MAX( 1, K )
       ELSE IF ( RIGHT ) THEN
-         Q = N
+         LDVQ = MAX( 1, N )
+         LDAQ = MAX( 1, M )
       END IF
       IF( .NOT.LEFT .AND. .NOT.RIGHT ) THEN
          INFO = -1
          INFO = -5
       ELSE IF( L.LT.0 .OR. L.GT.K ) THEN
          INFO = -6         
-      ELSE IF( NB.LT.1 .OR. NB.GT.K ) THEN
+      ELSE IF( NB.LT.1 .OR. (NB.GT.K .AND. K.GT.0) ) THEN
          INFO = -7
-      ELSE IF( LDV.LT.MAX( 1, Q ) ) THEN
+      ELSE IF( LDV.LT.LDVQ ) THEN
          INFO = -9
       ELSE IF( LDT.LT.NB ) THEN
          INFO = -11
-      ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
+      ELSE IF( LDA.LT.LDAQ ) THEN
          INFO = -13
       ELSE IF( LDB.LT.MAX( 1, M ) ) THEN
          INFO = -15
index f29ddf9..c5e828b 100644 (file)
          INFO = -1
       ELSE IF( N.LT.0 ) THEN
          INFO = -2
-      ELSE IF( L.LT.0 .OR. L.GT.MIN(M,N) ) THEN
+      ELSE IF( L.LT.0 .OR. (L.GT.MIN(M,N) .AND. MIN(M,N).GT.0)) THEN
          INFO = -3
-      ELSE IF( NB.LT.1 .OR. NB.GT.N ) THEN
+      ELSE IF( NB.LT.1 .OR. (NB.GT.N .AND. N.GT.0)) THEN
          INFO = -4
       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
          INFO = -6
index 248070a..4af6817 100644 (file)
             THETA(I) = ATAN2( SCNRM2( M-P-I+1, X21(I,I), 1 ),
      $                 SCNRM2( P-I+1, X11(I,I), 1 ) )
 *
-            CALL CLARFGP( P-I+1, X11(I,I), X11(I+1,I), 1, TAUP1(I) )
+            IF( P .GT. I ) THEN
+               CALL CLARFGP( P-I+1, X11(I,I), X11(I+1,I), 1, TAUP1(I) )
+            ELSE IF ( P .EQ. I ) THEN
+               CALL CLARFGP( P-I+1, X11(I,I), X11(I,I), 1, TAUP1(I) )
+            END IF
             X11(I,I) = ONE
-            CALL CLARFGP( M-P-I+1, X21(I,I), X21(I+1,I), 1, TAUP2(I) )
+            IF ( M-P .GT. I ) THEN
+               CALL CLARFGP( M-P-I+1, X21(I,I), X21(I+1,I), 1, 
+     $                       TAUP2(I) )
+            ELSE IF ( M-P .EQ. I ) THEN
+               CALL CLARFGP( M-P-I+1, X21(I,I), X21(I,I), 1,
+     $                       TAUP2(I) )
+            END IF
             X21(I,I) = ONE
 *
-            CALL CLARF( 'L', P-I+1, Q-I, X11(I,I), 1, CONJG(TAUP1(I)),
-     $                  X11(I,I+1), LDX11, WORK )
-            CALL CLARF( 'L', P-I+1, M-Q-I+1, X11(I,I), 1,
-     $                  CONJG(TAUP1(I)), X12(I,I), LDX12, WORK )
-            CALL CLARF( 'L', M-P-I+1, Q-I, X21(I,I), 1, CONJG(TAUP2(I)),
-     $                  X21(I,I+1), LDX21, WORK )
-            CALL CLARF( 'L', M-P-I+1, M-Q-I+1, X21(I,I), 1,
-     $                  CONJG(TAUP2(I)), X22(I,I), LDX22, WORK )
+            IF ( Q .GT. I ) THEN
+               CALL CLARF( 'L', P-I+1, Q-I, X11(I,I), 1, 
+     $                     CONJG(TAUP1(I)), X11(I,I+1), LDX11, WORK )
+               CALL CLARF( 'L', M-P-I+1, Q-I, X21(I,I), 1,
+     $                     CONJG(TAUP2(I)), X21(I,I+1), LDX21, WORK )
+            END IF
+            IF ( M-Q+1 .GT. I ) THEN
+               CALL CLARF( 'L', P-I+1, M-Q-I+1, X11(I,I), 1,
+     $                     CONJG(TAUP1(I)), X12(I,I), LDX12, WORK )
+               CALL CLARF( 'L', M-P-I+1, M-Q-I+1, X21(I,I), 1,
+     $                     CONJG(TAUP2(I)), X22(I,I), LDX22, WORK )
+            END IF
 *
             IF( I .LT. Q ) THEN
                CALL CSCAL( Q-I, CMPLX( -Z1*Z3*SIN(THETA(I)), 0.0E0 ),
 *
             IF( I .LT. Q ) THEN
                CALL CLACGV( Q-I, X11(I,I+1), LDX11 )
-               CALL CLARFGP( Q-I, X11(I,I+1), X11(I,I+2), LDX11,
-     $                       TAUQ1(I) )
+               IF ( I .EQ. Q-1 ) THEN
+                  CALL CLARFGP( Q-I, X11(I,I+1), X11(I,I+1), LDX11,
+     $                          TAUQ1(I) )
+               ELSE
+                  CALL CLARFGP( Q-I, X11(I,I+1), X11(I,I+2), LDX11,
+     $                          TAUQ1(I) )
+               END IF
                X11(I,I+1) = ONE
             END IF
-            CALL CLACGV( M-Q-I+1, X12(I,I), LDX12 )
-            CALL CLARFGP( M-Q-I+1, X12(I,I), X12(I,I+1), LDX12,
-     $                    TAUQ2(I) )
+            IF ( M-Q+1 .GT. I ) THEN
+               CALL CLACGV( M-Q-I+1, X12(I,I), LDX12 )
+               IF ( M-Q .EQ. I ) THEN
+                  CALL CLARFGP( M-Q-I+1, X12(I,I), X12(I,I), LDX12,
+     $                          TAUQ2(I) )
+               ELSE
+                  CALL CLARFGP( M-Q-I+1, X12(I,I), X12(I,I+1), LDX12,
+     $                          TAUQ2(I) )
+               END IF
+            END IF
             X12(I,I) = ONE
 *
             IF( I .LT. Q ) THEN
                CALL CLARF( 'R', M-P-I, Q-I, X11(I,I+1), LDX11, TAUQ1(I),
      $                     X21(I+1,I+1), LDX21, WORK )
             END IF
-            CALL CLARF( 'R', P-I, M-Q-I+1, X12(I,I), LDX12, TAUQ2(I),
-     $                  X12(I+1,I), LDX12, WORK )
-            CALL CLARF( 'R', M-P-I, M-Q-I+1, X12(I,I), LDX12, TAUQ2(I),
-     $                  X22(I+1,I), LDX22, WORK )
+            IF ( P .GT. I ) THEN
+               CALL CLARF( 'R', P-I, M-Q-I+1, X12(I,I), LDX12, TAUQ2(I),
+     $                     X12(I+1,I), LDX12, WORK )
+            END IF
+            IF ( M-P .GT. I ) THEN
+               CALL CLARF( 'R', M-P-I, M-Q-I+1, X12(I,I), LDX12,
+     $                     TAUQ2(I), X22(I+1,I), LDX22, WORK )
+            END IF
 *
             IF( I .LT. Q )
      $         CALL CLACGV( Q-I, X11(I,I+1), LDX11 )
             CALL CSCAL( M-Q-I+1, CMPLX( -Z1*Z4, 0.0E0 ), X12(I,I),
      $                  LDX12 )
             CALL CLACGV( M-Q-I+1, X12(I,I), LDX12 )
-            CALL CLARFGP( M-Q-I+1, X12(I,I), X12(I,I+1), LDX12,
-     $                    TAUQ2(I) )
+            IF ( I .GE. M-Q ) THEN
+               CALL CLARFGP( M-Q-I+1, X12(I,I), X12(I,I), LDX12,
+     $                       TAUQ2(I) )
+            ELSE
+               CALL CLARFGP( M-Q-I+1, X12(I,I), X12(I,I+1), LDX12,
+     $                       TAUQ2(I) )
+            END IF
             X12(I,I) = ONE
 *
-            CALL CLARF( 'R', P-I, M-Q-I+1, X12(I,I), LDX12, TAUQ2(I),
-     $                  X12(I+1,I), LDX12, WORK )
+            IF ( P .GT. I ) THEN
+               CALL CLARF( 'R', P-I, M-Q-I+1, X12(I,I), LDX12, TAUQ2(I),
+     $                     X12(I+1,I), LDX12, WORK )
+            END IF
             IF( M-P-Q .GE. 1 )
      $         CALL CLARF( 'R', M-P-Q, M-Q-I+1, X12(I,I), LDX12,
      $                     TAUQ2(I), X22(Q+1,I), LDX22, WORK )
 *
             CALL CLARFGP( P-I+1, X11(I,I), X11(I,I+1), LDX11, TAUP1(I) )
             X11(I,I) = ONE
-            CALL CLARFGP( M-P-I+1, X21(I,I), X21(I,I+1), LDX21,
-     $                    TAUP2(I) )
+            IF ( I .EQ. M-P ) THEN
+               CALL CLARFGP( M-P-I+1, X21(I,I), X21(I,I), LDX21,
+     $                       TAUP2(I) )
+            ELSE
+               CALL CLARFGP( M-P-I+1, X21(I,I), X21(I,I+1), LDX21,
+     $                       TAUP2(I) )
+            END IF
             X21(I,I) = ONE
 *
             CALL CLARF( 'R', Q-I, P-I+1, X11(I,I), LDX11, TAUP1(I),
             END IF
             CALL CLARF( 'L', M-Q-I+1, P-I, X12(I,I), 1, CONJG(TAUQ2(I)),
      $                  X12(I,I+1), LDX12, WORK )
-            CALL CLARF( 'L', M-Q-I+1, M-P-I, X12(I,I), 1,
-     $                  CONJG(TAUQ2(I)), X22(I,I+1), LDX22, WORK )
-*
+
+            IF ( M-P .GT. I ) THEN
+               CALL CLARF( 'L', M-Q-I+1, M-P-I, X12(I,I), 1,
+     $                     CONJG(TAUQ2(I)), X22(I,I+1), LDX22, WORK )
+            END IF
          END DO
 *
 *        Reduce columns Q + 1, ..., P of X12, X22
             CALL CLARFGP( M-Q-I+1, X12(I,I), X12(I+1,I), 1, TAUQ2(I) )
             X12(I,I) = ONE
 *
-            CALL CLARF( 'L', M-Q-I+1, P-I, X12(I,I), 1, CONJG(TAUQ2(I)),
-     $                  X12(I,I+1), LDX12, WORK )
+            IF ( P .GT. I ) THEN
+               CALL CLARF( 'L', M-Q-I+1, P-I, X12(I,I), 1,
+     $                     CONJG(TAUQ2(I)), X12(I,I+1), LDX12, WORK )
+            END IF
             IF( M-P-Q .GE. 1 )
      $         CALL CLARF( 'L', M-Q-I+1, M-P-Q, X12(I,I), 1,
      $                     CONJG(TAUQ2(I)), X22(I,Q+1), LDX22, WORK )
             CALL CLARFGP( M-P-Q-I+1, X22(P+I,Q+I), X22(P+I+1,Q+I), 1,
      $                    TAUQ2(P+I) )
             X22(P+I,Q+I) = ONE
-*
-            CALL CLARF( 'L', M-P-Q-I+1, M-P-Q-I, X22(P+I,Q+I), 1,
-     $                  CONJG(TAUQ2(P+I)), X22(P+I,Q+I+1), LDX22, WORK )
-*
+            IF ( M-P-Q .NE. I ) THEN
+               CALL CLARF( 'L', M-P-Q-I+1, M-P-Q-I, X22(P+I,Q+I), 1,
+     $                     CONJG(TAUQ2(P+I)), X22(P+I,Q+I+1), LDX22,
+     $                     WORK )
+            END IF
          END DO
 *
       END IF
index fef03dd..c9b36fb 100644 (file)
      $                   LBBCSDWORKOPT, LORBDBWORK, LORBDBWORKMIN,
      $                   LORBDBWORKOPT, LORGLQWORK, LORGLQWORKMIN,
      $                   LORGLQWORKOPT, LORGQRWORK, LORGQRWORKMIN,
-     $                   LORGQRWORKOPT, LWORKMIN, LWORKOPT
+     $                   LORGQRWORKOPT, LWORKMIN, LWORKOPT, P1, Q1
       LOGICAL            COLMAJOR, DEFAULTSIGNS, LQUERY, WANTU1, WANTU2,
      $                   WANTV1T, WANTV2T
       INTEGER            LRWORKMIN, LRWORKOPT
          INFO = -8
       ELSE IF( Q .LT. 0 .OR. Q .GT. M ) THEN
          INFO = -9
-      ELSE IF( ( COLMAJOR .AND. LDX11 .LT. MAX(1,P) ) .OR.
-     $         ( .NOT.COLMAJOR .AND. LDX11 .LT. MAX(1,Q) ) ) THEN
-         INFO = -11
+      ELSE IF ( COLMAJOR .AND.  LDX11 .LT. MAX( 1, P ) ) THEN
+        INFO = -11
+      ELSE IF (.NOT. COLMAJOR .AND. LDX11 .LT. MAX( 1, Q ) ) THEN
+        INFO = -11
+      ELSE IF (COLMAJOR .AND. LDX12 .LT. MAX( 1, P ) ) THEN
+        INFO = -13
+      ELSE IF (.NOT. COLMAJOR .AND. LDX12 .LT. MAX( 1, M-Q ) ) THEN
+        INFO = -13
+      ELSE IF (COLMAJOR .AND. LDX21 .LT. MAX( 1, M-P ) ) THEN
+        INFO = -15
+      ELSE IF (.NOT. COLMAJOR .AND. LDX21 .LT. MAX( 1, Q ) ) THEN
+        INFO = -15
+      ELSE IF (COLMAJOR .AND. LDX22 .LT. MAX( 1, M-P ) ) THEN
+        INFO = -17
+      ELSE IF (.NOT. COLMAJOR .AND. LDX22 .LT. MAX( 1, M-Q ) ) THEN
+        INFO = -17
       ELSE IF( WANTU1 .AND. LDU1 .LT. P ) THEN
          INFO = -20
       ELSE IF( WANTU2 .AND. LDU2 .LT. M-P ) THEN
          IB22D = IB21E + MAX( 1, Q - 1 )
          IB22E = IB22D + MAX( 1, Q )
          IBBCSD = IB22E + MAX( 1, Q - 1 )
-         CALL CBBCSD( JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS, M, P, Q, 0,
-     $                0, U1, LDU1, U2, LDU2, V1T, LDV1T, V2T, LDV2T, 0,
-     $                0, 0, 0, 0, 0, 0, 0, RWORK, -1, CHILDINFO )
+         CALL CBBCSD( JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS, M, P, Q, 
+     $                THETA, THETA, U1, LDU1, U2, LDU2, V1T, LDV1T,
+     $                V2T, LDV2T, THETA, THETA, THETA, THETA, THETA,
+     $                THETA, THETA, THETA, RWORK, -1, CHILDINFO )
          LBBCSDWORKOPT = INT( RWORK(1) )
          LBBCSDWORKMIN = LBBCSDWORKOPT
          LRWORKOPT = IBBCSD + LBBCSDWORKOPT - 1
          ITAUQ1 = ITAUP2 + MAX( 1, M - P )
          ITAUQ2 = ITAUQ1 + MAX( 1, Q )
          IORGQR = ITAUQ2 + MAX( 1, M - Q )
-         CALL CUNGQR( M-Q, M-Q, M-Q, 0, MAX(1,M-Q), 0, WORK, -1,
+         CALL CUNGQR( M-Q, M-Q, M-Q, 0, MAX(1,M-Q), U1, WORK, -1,
      $                CHILDINFO )
          LORGQRWORKOPT = INT( WORK(1) )
          LORGQRWORKMIN = MAX( 1, M - Q )
          IORGLQ = ITAUQ2 + MAX( 1, M - Q )
-         CALL CUNGLQ( M-Q, M-Q, M-Q, 0, MAX(1,M-Q), 0, WORK, -1,
+         CALL CUNGLQ( M-Q, M-Q, M-Q, 0, MAX(1,M-Q), U1, WORK, -1,
      $                CHILDINFO )
          LORGLQWORKOPT = INT( WORK(1) )
          LORGLQWORKMIN = MAX( 1, M - Q )
          IORBDB = ITAUQ2 + MAX( 1, M - Q )
          CALL CUNBDB( TRANS, SIGNS, M, P, Q, X11, LDX11, X12, LDX12,
-     $                X21, LDX21, X22, LDX22, 0, 0, 0, 0, 0, 0, WORK,
-     $                -1, CHILDINFO )
+     $                X21, LDX21, X22, LDX22, THETA, THETA, U1, U2,
+     $                V1T, V2T, WORK, -1, CHILDINFO )
          LORBDBWORKOPT = INT( WORK(1) )
          LORBDBWORKMIN = LORBDBWORKOPT
          LWORKOPT = MAX( IORGQR + LORGQRWORKOPT, IORGLQ + LORGLQWORKOPT,
          END IF
          IF( WANTV2T .AND. M-Q .GT. 0 ) THEN
             CALL CLACPY( 'U', P, M-Q, X12, LDX12, V2T, LDV2T )
-            CALL CLACPY( 'U', M-P-Q, M-P-Q, X22(Q+1,P+1), LDX22,
-     $                   V2T(P+1,P+1), LDV2T )
-            CALL CUNGLQ( M-Q, M-Q, M-Q, V2T, LDV2T, WORK(ITAUQ2),
-     $                   WORK(IORGLQ), LORGLQWORK, INFO )
+            IF( M-P .GT. Q ) THEN
+               CALL CLACPY( 'U', M-P-Q, M-P-Q, X22(Q+1,P+1), LDX22,
+     $                      V2T(P+1,P+1), LDV2T )
+            END IF
+            IF( M .GT. Q ) THEN
+               CALL CUNGLQ( M-Q, M-Q, M-Q, V2T, LDV2T, WORK(ITAUQ2),
+     $                      WORK(IORGLQ), LORGLQWORK, INFO )
+            END IF
          END IF
       ELSE
          IF( WANTU1 .AND. P .GT. 0 ) THEN
      $                   WORK(IORGQR), LORGQRWORK, INFO )
          END IF
          IF( WANTV2T .AND. M-Q .GT. 0 ) THEN
+            P1 = MIN( P+1, M )
+            Q1 = MIN( Q+1, M )
             CALL CLACPY( 'L', M-Q, P, X12, LDX12, V2T, LDV2T )
-            CALL CLACPY( 'L', M-P-Q, M-P-Q, X22(P+1,Q+1), LDX22,
-     $                   V2T(P+1,P+1), LDV2T )
+            IF ( M .GT. P+Q ) THEN
+               CALL CLACPY( 'L', M-P-Q, M-P-Q, X22(P1,Q1), LDX22,
+     $                      V2T(P+1,P+1), LDV2T )
+            END IF
             CALL CUNGQR( M-Q, M-Q, M-Q, V2T, LDV2T, WORK(ITAUQ2),
      $                   WORK(IORGQR), LORGQRWORK, INFO )
          END IF
index 775643b..bb60bcd 100644 (file)
       PARAMETER          ( HUNDRED = 100.0D0, MEIGHTH = -0.125D0,
      $                     ONE = 1.0D0, PIOVER2 = 1.57079632679489662D0,
      $                     TEN = 10.0D0, ZERO = 0.0D0 )
-      DOUBLE PRECISION   NEGONECOMPLEX
-      PARAMETER          ( NEGONECOMPLEX = -1.0D0 )
+      DOUBLE PRECISION   NEGONE
+      PARAMETER          ( NEGONE = -1.0D0 )
 *     ..
 *     .. Local Scalars ..
       LOGICAL            COLMAJOR, LQUERY, RESTART11, RESTART12,
 *     Initial deflation
 *
       IMAX = Q
-      DO WHILE((IMAX.GT.1).AND.(PHI(MAX(IMAX-1,1)).EQ.ZERO))
+      DO WHILE( IMAX .GT. 1 )
+         IF( PHI(IMAX-1) .NE. ZERO ) THEN
+            EXIT
+         END IF
          IMAX = IMAX - 1
       END DO
       IMIN = IMAX - 1
             B21D(IMAX) = -B21D(IMAX)
             IF( WANTV1T ) THEN
                IF( COLMAJOR ) THEN
-                  CALL DSCAL( Q, NEGONECOMPLEX, V1T(IMAX,1), LDV1T )
+                  CALL DSCAL( Q, NEGONE, V1T(IMAX,1), LDV1T )
                ELSE
-                  CALL DSCAL( Q, NEGONECOMPLEX, V1T(1,IMAX), 1 )
+                  CALL DSCAL( Q, NEGONE, V1T(1,IMAX), 1 )
                END IF
             END IF
          END IF
             B12D(IMAX) = -B12D(IMAX)
             IF( WANTU1 ) THEN
                IF( COLMAJOR ) THEN
-                  CALL DSCAL( P, NEGONECOMPLEX, U1(1,IMAX), 1 )
+                  CALL DSCAL( P, NEGONE, U1(1,IMAX), 1 )
                ELSE
-                  CALL DSCAL( P, NEGONECOMPLEX, U1(IMAX,1), LDU1 )
+                  CALL DSCAL( P, NEGONE, U1(IMAX,1), LDU1 )
                END IF
             END IF
          END IF
             B22D(IMAX) = -B22D(IMAX)
             IF( WANTU2 ) THEN
                IF( COLMAJOR ) THEN
-                  CALL DSCAL( M-P, NEGONECOMPLEX, U2(1,IMAX), 1 )
+                  CALL DSCAL( M-P, NEGONE, U2(1,IMAX), 1 )
                ELSE
-                  CALL DSCAL( M-P, NEGONECOMPLEX, U2(IMAX,1), LDU2 )
+                  CALL DSCAL( M-P, NEGONE, U2(IMAX,1), LDU2 )
                END IF
             END IF
          END IF
          IF( B12D(IMAX)+B22D(IMAX) .LT. 0 ) THEN
             IF( WANTV2T ) THEN
                IF( COLMAJOR ) THEN
-                  CALL DSCAL( M-Q, NEGONECOMPLEX, V2T(IMAX,1), LDV2T )
+                  CALL DSCAL( M-Q, NEGONE, V2T(IMAX,1), LDV2T )
                ELSE
-                  CALL DSCAL( M-Q, NEGONECOMPLEX, V2T(1,IMAX), 1 )
+                  CALL DSCAL( M-Q, NEGONE, V2T(1,IMAX), 1 )
                END IF
             END IF
          END IF
index 99b3968..753c843 100644 (file)
          INFO = -4
       ELSE IF( K.LT.0 .OR. K.GT.Q ) THEN
          INFO = -5
-      ELSE IF( NB.LT.1 .OR. NB.GT.K ) THEN
+      ELSE IF( NB.LT.1 .OR. (NB.GT.K .AND. K.GT.0)) THEN
          INFO = -6
       ELSE IF( LDV.LT.MAX( 1, Q ) ) THEN
          INFO = -8
index 5580f4e..6a5072c 100644 (file)
             THETA(I) = ATAN2( DNRM2( M-P-I+1, X21(I,I), 1 ),
      $                 DNRM2( P-I+1, X11(I,I), 1 ) )
 *
-            CALL DLARFGP( P-I+1, X11(I,I), X11(I+1,I), 1, TAUP1(I) )
+            IF( P .GT. I ) THEN
+               CALL DLARFGP( P-I+1, X11(I,I), X11(I+1,I), 1, TAUP1(I) )
+            ELSE IF( P .EQ. I ) THEN
+               CALL DLARFGP( P-I+1, X11(I,I), X11(I,I), 1, TAUP1(I) )
+            END IF
             X11(I,I) = ONE
-            CALL DLARFGP( M-P-I+1, X21(I,I), X21(I+1,I), 1, TAUP2(I) )
+            IF ( M-P .GT. I ) THEN
+               CALL DLARFGP( M-P-I+1, X21(I,I), X21(I+1,I), 1,
+     $                       TAUP2(I) )
+            ELSE IF ( M-P .EQ. I ) THEN
+               CALL DLARFGP( M-P-I+1, X21(I,I), X21(I,I), 1, TAUP2(I) )
+            END IF
             X21(I,I) = ONE
 *
-            CALL DLARF( 'L', P-I+1, Q-I, X11(I,I), 1, TAUP1(I),
-     $                  X11(I,I+1), LDX11, WORK )
-            CALL DLARF( 'L', P-I+1, M-Q-I+1, X11(I,I), 1, TAUP1(I),
-     $                  X12(I,I), LDX12, WORK )
-            CALL DLARF( 'L', M-P-I+1, Q-I, X21(I,I), 1, TAUP2(I),
-     $                  X21(I,I+1), LDX21, WORK )
-            CALL DLARF( 'L', M-P-I+1, M-Q-I+1, X21(I,I), 1, TAUP2(I),
-     $                  X22(I,I), LDX22, WORK )
+            IF ( Q .GT. I ) THEN
+               CALL DLARF( 'L', P-I+1, Q-I, X11(I,I), 1, TAUP1(I),
+     $                     X11(I,I+1), LDX11, WORK )
+            END IF
+            IF ( M-Q+1 .GT. I ) THEN
+               CALL DLARF( 'L', P-I+1, M-Q-I+1, X11(I,I), 1, TAUP1(I),
+     $                     X12(I,I), LDX12, WORK )
+            END IF
+            IF ( Q .GT. I ) THEN
+               CALL DLARF( 'L', M-P-I+1, Q-I, X21(I,I), 1, TAUP2(I),
+     $                     X21(I,I+1), LDX21, WORK )
+            END IF
+            IF ( M-Q+1 .GT. I ) THEN
+               CALL DLARF( 'L', M-P-I+1, M-Q-I+1, X21(I,I), 1, TAUP2(I),
+     $                     X22(I,I), LDX22, WORK )
+            END IF
 *
             IF( I .LT. Q ) THEN
                CALL DSCAL( Q-I, -Z1*Z3*SIN(THETA(I)), X11(I,I+1),
      $                  DNRM2( M-Q-I+1, X12(I,I), LDX12 ) )
 *
             IF( I .LT. Q ) THEN
-               CALL DLARFGP( Q-I, X11(I,I+1), X11(I,I+2), LDX11,
-     $                       TAUQ1(I) )
+               IF ( Q-I .EQ. 1 ) THEN
+                  CALL DLARFGP( Q-I, X11(I,I+1), X11(I,I+1), LDX11,
+     $                          TAUQ1(I) )
+               ELSE
+                  CALL DLARFGP( Q-I, X11(I,I+1), X11(I,I+2), LDX11,
+     $                          TAUQ1(I) )
+               END IF
                X11(I,I+1) = ONE
             END IF
-            CALL DLARFGP( M-Q-I+1, X12(I,I), X12(I,I+1), LDX12,
-     $                    TAUQ2(I) )
+            IF ( Q+I-1 .LT. M ) THEN
+               IF ( M-Q .EQ. I ) THEN
+                  CALL DLARFGP( M-Q-I+1, X12(I,I), X12(I,I), LDX12,
+     $                          TAUQ2(I) )
+               ELSE
+                  CALL DLARFGP( M-Q-I+1, X12(I,I), X12(I,I+1), LDX12,
+     $                          TAUQ2(I) )
+               END IF
+            END IF
             X12(I,I) = ONE
 *
             IF( I .LT. Q ) THEN
                CALL DLARF( 'R', M-P-I, Q-I, X11(I,I+1), LDX11, TAUQ1(I),
      $                     X21(I+1,I+1), LDX21, WORK )
             END IF
-            CALL DLARF( 'R', P-I, M-Q-I+1, X12(I,I), LDX12, TAUQ2(I),
-     $                  X12(I+1,I), LDX12, WORK )
-            CALL DLARF( 'R', M-P-I, M-Q-I+1, X12(I,I), LDX12, TAUQ2(I),
-     $                  X22(I+1,I), LDX22, WORK )
+            IF ( P .GT. I ) THEN
+               CALL DLARF( 'R', P-I, M-Q-I+1, X12(I,I), LDX12, TAUQ2(I),
+     $                     X12(I+1,I), LDX12, WORK )
+            END IF
+            IF ( M-P .GT. I ) THEN
+               CALL DLARF( 'R', M-P-I, M-Q-I+1, X12(I,I), LDX12,
+     $                     TAUQ2(I), X22(I+1,I), LDX22, WORK )
+            END IF
 *
          END DO
 *
          DO I = Q + 1, P
 *
             CALL DSCAL( M-Q-I+1, -Z1*Z4, X12(I,I), LDX12 )
-            CALL DLARFGP( M-Q-I+1, X12(I,I), X12(I,I+1), LDX12,
-     $                    TAUQ2(I) )
+            IF ( I .GE. M-Q ) THEN
+               CALL DLARFGP( M-Q-I+1, X12(I,I), X12(I,I), LDX12,
+     $                       TAUQ2(I) )
+            ELSE
+               CALL DLARFGP( M-Q-I+1, X12(I,I), X12(I,I+1), LDX12,
+     $                       TAUQ2(I) )
+            END IF
             X12(I,I) = ONE
 *
-            CALL DLARF( 'R', P-I, M-Q-I+1, X12(I,I), LDX12, TAUQ2(I),
-     $                  X12(I+1,I), LDX12, WORK )
+            IF ( P. GT. I ) THEN
+               CALL DLARF( 'R', P-I, M-Q-I+1, X12(I,I), LDX12, TAUQ2(I),
+     $                     X12(I+1,I), LDX12, WORK )
+            END IF
             IF( M-P-Q .GE. 1 )
      $         CALL DLARF( 'R', M-P-Q, M-Q-I+1, X12(I,I), LDX12,
      $                     TAUQ2(I), X22(Q+1,I), LDX22, WORK )
          DO I = 1, M - P - Q
 *
             CALL DSCAL( M-P-Q-I+1, Z2*Z4, X22(Q+I,P+I), LDX22 )
-            CALL DLARFGP( M-P-Q-I+1, X22(Q+I,P+I), X22(Q+I,P+I+1),
-     $                    LDX22, TAUQ2(P+I) )
+            IF ( I .EQ. M-P-Q ) THEN
+               CALL DLARFGP( M-P-Q-I+1, X22(Q+I,P+I), X22(Q+I,P+I),
+     $                       LDX22, TAUQ2(P+I) )
+            ELSE
+               CALL DLARFGP( M-P-Q-I+1, X22(Q+I,P+I), X22(Q+I,P+I+1),
+     $                       LDX22, TAUQ2(P+I) )
+            END IF
             X22(Q+I,P+I) = ONE
-            CALL DLARF( 'R', M-P-Q-I, M-P-Q-I+1, X22(Q+I,P+I), LDX22,
-     $                  TAUQ2(P+I), X22(Q+I+1,P+I), LDX22, WORK )
+            IF ( I .LT. M-P-Q ) THEN
+               CALL DLARF( 'R', M-P-Q-I, M-P-Q-I+1, X22(Q+I,P+I), LDX22,
+     $                     TAUQ2(P+I), X22(Q+I+1,P+I), LDX22, WORK )
+            END IF
 *
          END DO
 *
 *
             CALL DLARFGP( P-I+1, X11(I,I), X11(I,I+1), LDX11, TAUP1(I) )
             X11(I,I) = ONE
-            CALL DLARFGP( M-P-I+1, X21(I,I), X21(I,I+1), LDX21,
+            IF ( I .EQ. M-P ) THEN
+               CALL DLARFGP( M-P-I+1, X21(I,I), X21(I,I), LDX21,
+     $                    TAUP2(I) )
+            ELSE
+               CALL DLARFGP( M-P-I+1, X21(I,I), X21(I,I+1), LDX21,
      $                    TAUP2(I) )
+            END IF
             X21(I,I) = ONE
 *
-            CALL DLARF( 'R', Q-I, P-I+1, X11(I,I), LDX11, TAUP1(I),
-     $                  X11(I+1,I), LDX11, WORK )
-            CALL DLARF( 'R', M-Q-I+1, P-I+1, X11(I,I), LDX11, TAUP1(I),
-     $                  X12(I,I), LDX12, WORK )
-            CALL DLARF( 'R', Q-I, M-P-I+1, X21(I,I), LDX21, TAUP2(I),
-     $                  X21(I+1,I), LDX21, WORK )
-            CALL DLARF( 'R', M-Q-I+1, M-P-I+1, X21(I,I), LDX21,
-     $                  TAUP2(I), X22(I,I), LDX22, WORK )
+            IF ( Q .GT. I ) THEN
+               CALL DLARF( 'R', Q-I, P-I+1, X11(I,I), LDX11, TAUP1(I),
+     $                     X11(I+1,I), LDX11, WORK )
+            END IF
+            IF ( M-Q+1 .GT. I ) THEN
+               CALL DLARF( 'R', M-Q-I+1, P-I+1, X11(I,I), LDX11,
+     $                     TAUP1(I), X12(I,I), LDX12, WORK )
+            END IF
+            IF ( Q .GT. I ) THEN
+               CALL DLARF( 'R', Q-I, M-P-I+1, X21(I,I), LDX21, TAUP2(I),
+     $                     X21(I+1,I), LDX21, WORK )
+            END IF
+            IF ( M-Q+1 .GT. I ) THEN
+               CALL DLARF( 'R', M-Q-I+1, M-P-I+1, X21(I,I), LDX21,
+     $                     TAUP2(I), X22(I,I), LDX22, WORK )
+            END IF
 *
             IF( I .LT. Q ) THEN
                CALL DSCAL( Q-I, -Z1*Z3*SIN(THETA(I)), X11(I+1,I), 1 )
      $                  DNRM2( M-Q-I+1, X12(I,I), 1 ) )
 *
             IF( I .LT. Q ) THEN
-               CALL DLARFGP( Q-I, X11(I+1,I), X11(I+2,I), 1, TAUQ1(I) )
+               IF ( Q-I .EQ. 1) THEN
+                  CALL DLARFGP( Q-I, X11(I+1,I), X11(I+1,I), 1,
+     $                          TAUQ1(I) )
+               ELSE
+                  CALL DLARFGP( Q-I, X11(I+1,I), X11(I+2,I), 1,
+     $                          TAUQ1(I) )
+               END IF
                X11(I+1,I) = ONE
             END IF
-            CALL DLARFGP( M-Q-I+1, X12(I,I), X12(I+1,I), 1, TAUQ2(I) )
+            IF ( M-Q .GT. I ) THEN
+               CALL DLARFGP( M-Q-I+1, X12(I,I), X12(I+1,I), 1, 
+     $                       TAUQ2(I) )
+            ELSE
+               CALL DLARFGP( M-Q-I+1, X12(I,I), X12(I,I), 1, 
+     $                       TAUQ2(I) )
+            END IF               
             X12(I,I) = ONE
 *
             IF( I .LT. Q ) THEN
             END IF
             CALL DLARF( 'L', M-Q-I+1, P-I, X12(I,I), 1, TAUQ2(I),
      $                  X12(I,I+1), LDX12, WORK )
-            CALL DLARF( 'L', M-Q-I+1, M-P-I, X12(I,I), 1, TAUQ2(I),
-     $                  X22(I,I+1), LDX22, WORK )
+            IF ( M-P-I .GT. 0 ) THEN
+               CALL DLARF( 'L', M-Q-I+1, M-P-I, X12(I,I), 1, TAUQ2(I),
+     $                     X22(I,I+1), LDX22, WORK )
+            END IF
 *
          END DO
 *
             CALL DLARFGP( M-Q-I+1, X12(I,I), X12(I+1,I), 1, TAUQ2(I) )
             X12(I,I) = ONE
 *
-            CALL DLARF( 'L', M-Q-I+1, P-I, X12(I,I), 1, TAUQ2(I),
+            IF ( P .GT. I ) THEN 
+               CALL DLARF( 'L', M-Q-I+1, P-I, X12(I,I), 1, TAUQ2(I),
      $                  X12(I,I+1), LDX12, WORK )
+            END IF
             IF( M-P-Q .GE. 1 )
      $         CALL DLARF( 'L', M-Q-I+1, M-P-Q, X12(I,I), 1, TAUQ2(I),
      $                     X22(I,Q+1), LDX22, WORK )
          DO I = 1, M - P - Q
 *
             CALL DSCAL( M-P-Q-I+1, Z2*Z4, X22(P+I,Q+I), 1 )
-            CALL DLARFGP( M-P-Q-I+1, X22(P+I,Q+I), X22(P+I+1,Q+I), 1,
-     $                    TAUQ2(P+I) )
-            X22(P+I,Q+I) = ONE
-*
-            CALL DLARF( 'L', M-P-Q-I+1, M-P-Q-I, X22(P+I,Q+I), 1,
+            IF ( M-P-Q .EQ. I ) THEN
+               CALL DLARFGP( M-P-Q-I+1, X22(P+I,Q+I), X22(P+I,Q+I), 1,
+     $                       TAUQ2(P+I) )
+            ELSE
+               CALL DLARFGP( M-P-Q-I+1, X22(P+I,Q+I), X22(P+I+1,Q+I), 1,
+     $                       TAUQ2(P+I) )
+               CALL DLARF( 'L', M-P-Q-I+1, M-P-Q-I, X22(P+I,Q+I), 1,
      $                  TAUQ2(P+I), X22(P+I,Q+I+1), LDX22, WORK )
+            END IF
+            X22(P+I,Q+I) = ONE
 *
          END DO
 *
index 927bf66..2e753a8 100644 (file)
          INFO = -8
       ELSE IF( Q .LT. 0 .OR. Q .GT. M ) THEN
          INFO = -9
-      ELSE IF( ( COLMAJOR .AND. LDX11 .LT. MAX(1,P) ) .OR.
-     $         ( .NOT.COLMAJOR .AND. LDX11 .LT. MAX(1,Q) ) ) THEN
-         INFO = -11
+      ELSE IF ( COLMAJOR .AND.  LDX11 .LT. MAX( 1, P ) ) THEN
+        INFO = -11
+      ELSE IF (.NOT. COLMAJOR .AND. LDX11 .LT. MAX( 1, Q ) ) THEN
+        INFO = -11
+      ELSE IF (COLMAJOR .AND. LDX12 .LT. MAX( 1, P ) ) THEN
+        INFO = -13
+      ELSE IF (.NOT. COLMAJOR .AND. LDX12 .LT. MAX( 1, M-Q ) ) THEN
+        INFO = -13
+      ELSE IF (COLMAJOR .AND. LDX21 .LT. MAX( 1, M-P ) ) THEN
+        INFO = -15
+      ELSE IF (.NOT. COLMAJOR .AND. LDX21 .LT. MAX( 1, Q ) ) THEN
+        INFO = -15
+      ELSE IF (COLMAJOR .AND. LDX22 .LT. MAX( 1, M-P ) ) THEN
+        INFO = -17
+      ELSE IF (.NOT. COLMAJOR .AND. LDX22 .LT. MAX( 1, M-Q ) ) THEN
+        INFO = -17
       ELSE IF( WANTU1 .AND. LDU1 .LT. P ) THEN
          INFO = -20
       ELSE IF( WANTU2 .AND. LDU2 .LT. M-P ) THEN
          ITAUQ1 = ITAUP2 + MAX( 1, M - P )
          ITAUQ2 = ITAUQ1 + MAX( 1, Q )
          IORGQR = ITAUQ2 + MAX( 1, M - Q )
-         CALL DORGQR( M-Q, M-Q, M-Q, 0, MAX(1,M-Q), 0, WORK, -1,
+         CALL DORGQR( M-Q, M-Q, M-Q, U1, MAX(1,M-Q), U1, WORK, -1,
      $                CHILDINFO )
          LORGQRWORKOPT = INT( WORK(1) )
          LORGQRWORKMIN = MAX( 1, M - Q )
          IORGLQ = ITAUQ2 + MAX( 1, M - Q )
-         CALL DORGLQ( M-Q, M-Q, M-Q, 0, MAX(1,M-Q), 0, WORK, -1,
+         CALL DORGLQ( M-Q, M-Q, M-Q, U1, MAX(1,M-Q), U1, WORK, -1,
      $                CHILDINFO )
          LORGLQWORKOPT = INT( WORK(1) )
          LORGLQWORKMIN = MAX( 1, M - Q )
          IORBDB = ITAUQ2 + MAX( 1, M - Q )
          CALL DORBDB( TRANS, SIGNS, M, P, Q, X11, LDX11, X12, LDX12,
-     $                X21, LDX21, X22, LDX22, 0, 0, 0, 0, 0, 0, WORK,
-     $                -1, CHILDINFO )
+     $                X21, LDX21, X22, LDX22, THETA, V1T, U1, U2, V1T,
+     $                V2T, WORK, -1, CHILDINFO )
          LORBDBWORKOPT = INT( WORK(1) )
          LORBDBWORKMIN = LORBDBWORKOPT
          IB11D = ITAUQ2 + MAX( 1, M - Q )
          IB22D = IB21E + MAX( 1, Q - 1 )
          IB22E = IB22D + MAX( 1, Q )
          IBBCSD = IB22E + MAX( 1, Q - 1 )
-         CALL DBBCSD( JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS, M, P, Q, 0,
-     $                0, U1, LDU1, U2, LDU2, V1T, LDV1T, V2T, LDV2T, 0,
-     $                0, 0, 0, 0, 0, 0, 0, WORK, -1, CHILDINFO )
+         CALL DBBCSD( JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS, M, P, Q, 
+     $                THETA, THETA, U1, LDU1, U2, LDU2, V1T, LDV1T, V2T,
+     $                LDV2T, U1, U1, U1, U1, U1, U1, U1, U1, WORK, -1,
+     $                CHILDINFO )
          LBBCSDWORKOPT = INT( WORK(1) )
          LBBCSDWORKMIN = LBBCSDWORKOPT
          LWORKOPT = MAX( IORGQR + LORGQRWORKOPT, IORGLQ + LORGLQWORKOPT,
          END IF
          IF( WANTV2T .AND. M-Q .GT. 0 ) THEN
             CALL DLACPY( 'U', P, M-Q, X12, LDX12, V2T, LDV2T )
-            CALL DLACPY( 'U', M-P-Q, M-P-Q, X22(Q+1,P+1), LDX22,
-     $                   V2T(P+1,P+1), LDV2T )
-            CALL DORGLQ( M-Q, M-Q, M-Q, V2T, LDV2T, WORK(ITAUQ2),
-     $                   WORK(IORGLQ), LORGLQWORK, INFO )
+            IF (M-P .GT. Q) Then
+               CALL DLACPY( 'U', M-P-Q, M-P-Q, X22(Q+1,P+1), LDX22,
+     $                      V2T(P+1,P+1), LDV2T )
+            END IF
+            IF (M .GT. Q) THEN
+               CALL DORGLQ( M-Q, M-Q, M-Q, V2T, LDV2T, WORK(ITAUQ2),
+     $                      WORK(IORGLQ), LORGLQWORK, INFO )
+            END IF
          END IF
       ELSE
          IF( WANTU1 .AND. P .GT. 0 ) THEN
index baeffa0..709da33 100644 (file)
 *     ..
 *     .. Local Scalars ..
       LOGICAL            LEFT, RIGHT, TRAN, NOTRAN
-      INTEGER            I, IB, MB, LB, KF, Q
+      INTEGER            I, IB, MB, LB, KF, LDAQ, LDVQ
 *     ..
 *     .. External Functions ..
       LOGICAL            LSAME
       TRAN   = LSAME( TRANS, 'T' )
       NOTRAN = LSAME( TRANS, 'N' )
 *      
-      IF( LEFT ) THEN
-         Q = M
+      IF ( LEFT ) THEN
+         LDVQ = MAX( 1, M )
+         LDAQ = MAX( 1, K )
       ELSE IF ( RIGHT ) THEN
-         Q = N
+         LDVQ = MAX( 1, N )
+         LDAQ = MAX( 1, M )
       END IF
       IF( .NOT.LEFT .AND. .NOT.RIGHT ) THEN
          INFO = -1
          INFO = -5
       ELSE IF( L.LT.0 .OR. L.GT.K ) THEN
          INFO = -6         
-      ELSE IF( NB.LT.1 .OR. NB.GT.K ) THEN
+      ELSE IF( NB.LT.1 .OR. (NB.GT.K .AND. K.GT.0) ) THEN
          INFO = -7
-      ELSE IF( LDV.LT.MAX( 1, Q ) ) THEN
+      ELSE IF( LDV.LT.LDVQ ) THEN
          INFO = -9
       ELSE IF( LDT.LT.NB ) THEN
          INFO = -11
-      ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
+      ELSE IF( LDA.LT.LDAQ ) THEN
          INFO = -13
       ELSE IF( LDB.LT.MAX( 1, M ) ) THEN
          INFO = -15
index 95f6fde..25acef5 100644 (file)
          INFO = -1
       ELSE IF( N.LT.0 ) THEN
          INFO = -2
-      ELSE IF( L.LT.0 .OR. L.GT.MIN(M,N) ) THEN
+      ELSE IF( L.LT.0 .OR. (L.GT.MIN(M,N) .AND. MIN(M,N).GT.0)) THEN
          INFO = -3
-      ELSE IF( NB.LT.1 .OR. NB.GT.N ) THEN
+      ELSE IF( NB.LT.1 .OR. (NB.GT.N .AND. N.GT.0)) THEN
          INFO = -4
       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
          INFO = -6
index c304d07..b572681 100644 (file)
       PARAMETER          ( HUNDRED = 100.0E0, MEIGHTH = -0.125E0,
      $                     ONE = 1.0E0, PIOVER2 = 1.57079632679489662E0,
      $                     TEN = 10.0E0, ZERO = 0.0E0 )
-      REAL               NEGONECOMPLEX
-      PARAMETER          ( NEGONECOMPLEX = -1.0E0 )
+      REAL               NEGONE
+      PARAMETER          ( NEGONE = -1.0E0 )
 *     ..
 *     .. Local Scalars ..
       LOGICAL            COLMAJOR, LQUERY, RESTART11, RESTART12,
 *     Initial deflation
 *
       IMAX = Q
-      DO WHILE((IMAX.GT.1).AND.(PHI(MAX(IMAX-1,1)).EQ.ZERO))
+      DO WHILE( IMAX .GT. 1 )
+         IF( PHI(IMAX-1) .NE. ZERO ) THEN
+            EXIT
+         END IF
          IMAX = IMAX - 1
       END DO
       IMIN = IMAX - 1
             B21D(IMAX) = -B21D(IMAX)
             IF( WANTV1T ) THEN
                IF( COLMAJOR ) THEN
-                  CALL SSCAL( Q, NEGONECOMPLEX, V1T(IMAX,1), LDV1T )
+                  CALL SSCAL( Q, NEGONE, V1T(IMAX,1), LDV1T )
                ELSE
-                  CALL SSCAL( Q, NEGONECOMPLEX, V1T(1,IMAX), 1 )
+                  CALL SSCAL( Q, NEGONE, V1T(1,IMAX), 1 )
                END IF
             END IF
          END IF
             B12D(IMAX) = -B12D(IMAX)
             IF( WANTU1 ) THEN
                IF( COLMAJOR ) THEN
-                  CALL SSCAL( P, NEGONECOMPLEX, U1(1,IMAX), 1 )
+                  CALL SSCAL( P, NEGONE, U1(1,IMAX), 1 )
                ELSE
-                  CALL SSCAL( P, NEGONECOMPLEX, U1(IMAX,1), LDU1 )
+                  CALL SSCAL( P, NEGONE, U1(IMAX,1), LDU1 )
                END IF
             END IF
          END IF
             B22D(IMAX) = -B22D(IMAX)
             IF( WANTU2 ) THEN
                IF( COLMAJOR ) THEN
-                  CALL SSCAL( M-P, NEGONECOMPLEX, U2(1,IMAX), 1 )
+                  CALL SSCAL( M-P, NEGONE, U2(1,IMAX), 1 )
                ELSE
-                  CALL SSCAL( M-P, NEGONECOMPLEX, U2(IMAX,1), LDU2 )
+                  CALL SSCAL( M-P, NEGONE, U2(IMAX,1), LDU2 )
                END IF
             END IF
          END IF
          IF( B12D(IMAX)+B22D(IMAX) .LT. 0 ) THEN
             IF( WANTV2T ) THEN
                IF( COLMAJOR ) THEN
-                  CALL SSCAL( M-Q, NEGONECOMPLEX, V2T(IMAX,1), LDV2T )
+                  CALL SSCAL( M-Q, NEGONE, V2T(IMAX,1), LDV2T )
                ELSE
-                  CALL SSCAL( M-Q, NEGONECOMPLEX, V2T(1,IMAX), 1 )
+                  CALL SSCAL( M-Q, NEGONE, V2T(1,IMAX), 1 )
                END IF
             END IF
          END IF
index 443c9d1..4826f3b 100644 (file)
          INFO = -4
       ELSE IF( K.LT.0 .OR. K.GT.Q ) THEN
          INFO = -5
-      ELSE IF( NB.LT.1 .OR. NB.GT.K ) THEN
+      ELSE IF( NB.LT.1 .OR. (NB.GT.K .AND. K.GT.0)) THEN
          INFO = -6
       ELSE IF( LDV.LT.MAX( 1, Q ) ) THEN
          INFO = -8
index 433b3c0..9f406bf 100644 (file)
             THETA(I) = ATAN2( SNRM2( M-P-I+1, X21(I,I), 1 ),
      $                 SNRM2( P-I+1, X11(I,I), 1 ) )
 *
-            CALL SLARFGP( P-I+1, X11(I,I), X11(I+1,I), 1, TAUP1(I) )
+            IF( P .GT. I ) THEN
+               CALL SLARFGP( P-I+1, X11(I,I), X11(I+1,I), 1, TAUP1(I) )
+            ELSE IF( P .EQ. I ) THEN
+               CALL SLARFGP( P-I+1, X11(I,I), X11(I,I), 1, TAUP1(I) )
+            END IF
             X11(I,I) = ONE
-            CALL SLARFGP( M-P-I+1, X21(I,I), X21(I+1,I), 1, TAUP2(I) )
+            IF ( M-P .GT. I ) THEN
+               CALL SLARFGP( M-P-I+1, X21(I,I), X21(I+1,I), 1,
+     $                       TAUP2(I) )
+            ELSE IF ( M-P .EQ. I ) THEN
+               CALL SLARFGP( M-P-I+1, X21(I,I), X21(I,I), 1, TAUP2(I) )
+            END IF
             X21(I,I) = ONE
 *
-            CALL SLARF( 'L', P-I+1, Q-I, X11(I,I), 1, TAUP1(I),
-     $                  X11(I,I+1), LDX11, WORK )
-            CALL SLARF( 'L', P-I+1, M-Q-I+1, X11(I,I), 1, TAUP1(I),
-     $                  X12(I,I), LDX12, WORK )
-            CALL SLARF( 'L', M-P-I+1, Q-I, X21(I,I), 1, TAUP2(I),
-     $                  X21(I,I+1), LDX21, WORK )
-            CALL SLARF( 'L', M-P-I+1, M-Q-I+1, X21(I,I), 1, TAUP2(I),
-     $                  X22(I,I), LDX22, WORK )
+            IF ( Q .GT. I ) THEN
+               CALL SLARF( 'L', P-I+1, Q-I, X11(I,I), 1, TAUP1(I),
+     $                     X11(I,I+1), LDX11, WORK )
+            END IF
+            IF ( M-Q+1 .GT. I ) THEN
+               CALL SLARF( 'L', P-I+1, M-Q-I+1, X11(I,I), 1, TAUP1(I),
+     $                     X12(I,I), LDX12, WORK )
+            END IF
+            IF ( Q .GT. I ) THEN
+               CALL SLARF( 'L', M-P-I+1, Q-I, X21(I,I), 1, TAUP2(I),
+     $                     X21(I,I+1), LDX21, WORK )
+            END IF
+            IF ( M-Q+1 .GT. I ) THEN
+               CALL SLARF( 'L', M-P-I+1, M-Q-I+1, X21(I,I), 1, TAUP2(I),
+     $                     X22(I,I), LDX22, WORK )
+            END IF
 *
             IF( I .LT. Q ) THEN
                CALL SSCAL( Q-I, -Z1*Z3*SIN(THETA(I)), X11(I,I+1),
      $                  SNRM2( M-Q-I+1, X12(I,I), LDX12 ) )
 *
             IF( I .LT. Q ) THEN
-               CALL SLARFGP( Q-I, X11(I,I+1), X11(I,I+2), LDX11,
-     $                       TAUQ1(I) )
+               IF ( Q-I .EQ. 1 ) THEN
+                  CALL SLARFGP( Q-I, X11(I,I+1), X11(I,I+1), LDX11,
+     $                          TAUQ1(I) )
+               ELSE
+                  CALL SLARFGP( Q-I, X11(I,I+1), X11(I,I+2), LDX11,
+     $                          TAUQ1(I) )
+               END IF
                X11(I,I+1) = ONE
             END IF
-            CALL SLARFGP( M-Q-I+1, X12(I,I), X12(I,I+1), LDX12,
-     $                    TAUQ2(I) )
+            IF ( Q+I-1 .LT. M ) THEN
+               IF ( M-Q .EQ. I ) THEN
+                  CALL SLARFGP( M-Q-I+1, X12(I,I), X12(I,I), LDX12,
+     $                          TAUQ2(I) )
+               ELSE
+                  CALL SLARFGP( M-Q-I+1, X12(I,I), X12(I,I+1), LDX12,
+     $                          TAUQ2(I) )
+               END IF
+            END IF
             X12(I,I) = ONE
 *
             IF( I .LT. Q ) THEN
                CALL SLARF( 'R', M-P-I, Q-I, X11(I,I+1), LDX11, TAUQ1(I),
      $                     X21(I+1,I+1), LDX21, WORK )
             END IF
-            CALL SLARF( 'R', P-I, M-Q-I+1, X12(I,I), LDX12, TAUQ2(I),
-     $                  X12(I+1,I), LDX12, WORK )
-            CALL SLARF( 'R', M-P-I, M-Q-I+1, X12(I,I), LDX12, TAUQ2(I),
-     $                  X22(I+1,I), LDX22, WORK )
+            IF ( P .GT. I ) THEN
+               CALL SLARF( 'R', P-I, M-Q-I+1, X12(I,I), LDX12, TAUQ2(I),
+     $                     X12(I+1,I), LDX12, WORK )
+            END IF
+            IF ( M-P .GT. I ) THEN
+               CALL SLARF( 'R', M-P-I, M-Q-I+1, X12(I,I), LDX12,
+     $                     TAUQ2(I), X22(I+1,I), LDX22, WORK )
+            END IF
 *
          END DO
 *
          DO I = Q + 1, P
 *
             CALL SSCAL( M-Q-I+1, -Z1*Z4, X12(I,I), LDX12 )
-            CALL SLARFGP( M-Q-I+1, X12(I,I), X12(I,I+1), LDX12,
-     $                    TAUQ2(I) )
+            IF ( I .GE. M-Q ) THEN
+               CALL SLARFGP( M-Q-I+1, X12(I,I), X12(I,I), LDX12,
+     $                       TAUQ2(I) )
+            ELSE
+               CALL SLARFGP( M-Q-I+1, X12(I,I), X12(I,I+1), LDX12,
+     $                       TAUQ2(I) )
+            END IF
             X12(I,I) = ONE
 *
-            CALL SLARF( 'R', P-I, M-Q-I+1, X12(I,I), LDX12, TAUQ2(I),
-     $                  X12(I+1,I), LDX12, WORK )
+            IF ( P. GT. I ) THEN
+               CALL SLARF( 'R', P-I, M-Q-I+1, X12(I,I), LDX12, TAUQ2(I),
+     $                     X12(I+1,I), LDX12, WORK )
+            END IF
             IF( M-P-Q .GE. 1 )
      $         CALL SLARF( 'R', M-P-Q, M-Q-I+1, X12(I,I), LDX12,
      $                     TAUQ2(I), X22(Q+1,I), LDX22, WORK )
          DO I = 1, M - P - Q
 *
             CALL SSCAL( M-P-Q-I+1, Z2*Z4, X22(Q+I,P+I), LDX22 )
-            CALL SLARFGP( M-P-Q-I+1, X22(Q+I,P+I), X22(Q+I,P+I+1),
-     $                    LDX22, TAUQ2(P+I) )
+            IF ( I .EQ. M-P-Q ) THEN
+               CALL SLARFGP( M-P-Q-I+1, X22(Q+I,P+I), X22(Q+I,P+I),
+     $                       LDX22, TAUQ2(P+I) )
+            ELSE
+               CALL SLARFGP( M-P-Q-I+1, X22(Q+I,P+I), X22(Q+I,P+I+1),
+     $                       LDX22, TAUQ2(P+I) )
+            END IF
             X22(Q+I,P+I) = ONE
-            CALL SLARF( 'R', M-P-Q-I, M-P-Q-I+1, X22(Q+I,P+I), LDX22,
-     $                  TAUQ2(P+I), X22(Q+I+1,P+I), LDX22, WORK )
+            IF ( I .LT. M-P-Q ) THEN
+               CALL SLARF( 'R', M-P-Q-I, M-P-Q-I+1, X22(Q+I,P+I), LDX22,
+     $                     TAUQ2(P+I), X22(Q+I+1,P+I), LDX22, WORK )
+            END IF
 *
          END DO
 *
 *
             CALL SLARFGP( P-I+1, X11(I,I), X11(I,I+1), LDX11, TAUP1(I) )
             X11(I,I) = ONE
-            CALL SLARFGP( M-P-I+1, X21(I,I), X21(I,I+1), LDX21,
-     $                    TAUP2(I) )
+            IF ( I .EQ. M-P ) THEN
+               CALL SLARFGP( M-P-I+1, X21(I,I), X21(I,I), LDX21,
+     $                       TAUP2(I) )
+            ELSE
+               CALL SLARFGP( M-P-I+1, X21(I,I), X21(I,I+1), LDX21,
+     $                       TAUP2(I) )
+            END IF
             X21(I,I) = ONE
 *
-            CALL SLARF( 'R', Q-I, P-I+1, X11(I,I), LDX11, TAUP1(I),
-     $                  X11(I+1,I), LDX11, WORK )
-            CALL SLARF( 'R', M-Q-I+1, P-I+1, X11(I,I), LDX11, TAUP1(I),
-     $                  X12(I,I), LDX12, WORK )
-            CALL SLARF( 'R', Q-I, M-P-I+1, X21(I,I), LDX21, TAUP2(I),
-     $                  X21(I+1,I), LDX21, WORK )
-            CALL SLARF( 'R', M-Q-I+1, M-P-I+1, X21(I,I), LDX21,
-     $                  TAUP2(I), X22(I,I), LDX22, WORK )
+            IF ( Q .GT. I ) THEN
+               CALL SLARF( 'R', Q-I, P-I+1, X11(I,I), LDX11, TAUP1(I),
+     $                     X11(I+1,I), LDX11, WORK )
+            END IF
+            IF ( M-Q+1 .GT. I ) THEN
+               CALL SLARF( 'R', M-Q-I+1, P-I+1, X11(I,I), LDX11,
+     $                     TAUP1(I), X12(I,I), LDX12, WORK )
+            END IF
+            IF ( Q .GT. I ) THEN
+               CALL SLARF( 'R', Q-I, M-P-I+1, X21(I,I), LDX21, TAUP2(I),
+     $                     X21(I+1,I), LDX21, WORK )
+            END IF
+            IF ( M-Q+1 .GT. I ) THEN
+               CALL SLARF( 'R', M-Q-I+1, M-P-I+1, X21(I,I), LDX21,
+     $                     TAUP2(I), X22(I,I), LDX22, WORK )
+            END IF
 *
             IF( I .LT. Q ) THEN
                CALL SSCAL( Q-I, -Z1*Z3*SIN(THETA(I)), X11(I+1,I), 1 )
      $                  SNRM2( M-Q-I+1, X12(I,I), 1 ) )
 *
             IF( I .LT. Q ) THEN
-               CALL SLARFGP( Q-I, X11(I+1,I), X11(I+2,I), 1, TAUQ1(I) )
+               IF ( Q-I .EQ. 1) THEN
+                  CALL SLARFGP( Q-I, X11(I+1,I), X11(I+1,I), 1,
+     $                          TAUQ1(I) )
+               ELSE
+                  CALL SLARFGP( Q-I, X11(I+1,I), X11(I+2,I), 1,
+     $                          TAUQ1(I) )
+               END IF
                X11(I+1,I) = ONE
             END IF
-            CALL SLARFGP( M-Q-I+1, X12(I,I), X12(I+1,I), 1, TAUQ2(I) )
+            IF ( M-Q .GT. I ) THEN
+               CALL SLARFGP( M-Q-I+1, X12(I,I), X12(I+1,I), 1,
+     $                       TAUQ2(I) )
+            ELSE
+               CALL SLARFGP( M-Q-I+1, X12(I,I), X12(I,I), 1,
+     $                       TAUQ2(I) )
+            END IF
             X12(I,I) = ONE
 *
             IF( I .LT. Q ) THEN
             END IF
             CALL SLARF( 'L', M-Q-I+1, P-I, X12(I,I), 1, TAUQ2(I),
      $                  X12(I,I+1), LDX12, WORK )
-            CALL SLARF( 'L', M-Q-I+1, M-P-I, X12(I,I), 1, TAUQ2(I),
-     $                  X22(I,I+1), LDX22, WORK )
+            IF ( M-P-I .GT. 0 ) THEN
+               CALL SLARF( 'L', M-Q-I+1, M-P-I, X12(I,I), 1, TAUQ2(I),
+     $                     X22(I,I+1), LDX22, WORK )
+            END IF
 *
          END DO
 *
             CALL SLARFGP( M-Q-I+1, X12(I,I), X12(I+1,I), 1, TAUQ2(I) )
             X12(I,I) = ONE
 *
-            CALL SLARF( 'L', M-Q-I+1, P-I, X12(I,I), 1, TAUQ2(I),
-     $                  X12(I,I+1), LDX12, WORK )
+            IF ( P .GT. I ) THEN
+               CALL SLARF( 'L', M-Q-I+1, P-I, X12(I,I), 1, TAUQ2(I),
+     $                     X12(I,I+1), LDX12, WORK )
+            END IF
             IF( M-P-Q .GE. 1 )
      $         CALL SLARF( 'L', M-Q-I+1, M-P-Q, X12(I,I), 1, TAUQ2(I),
      $                     X22(I,Q+1), LDX22, WORK )
          DO I = 1, M - P - Q
 *
             CALL SSCAL( M-P-Q-I+1, Z2*Z4, X22(P+I,Q+I), 1 )
-            CALL SLARFGP( M-P-Q-I+1, X22(P+I,Q+I), X22(P+I+1,Q+I), 1,
-     $                    TAUQ2(P+I) )
-            X22(P+I,Q+I) = ONE
+            IF ( M-P-Q .EQ. I ) THEN
+               CALL SLARFGP( M-P-Q-I+1, X22(P+I,Q+I), X22(P+I,Q+I), 1,
+     $                       TAUQ2(P+I) )
+               X22(P+I,Q+I) = ONE
+            ELSE
+               CALL SLARFGP( M-P-Q-I+1, X22(P+I,Q+I), X22(P+I+1,Q+I), 1,
+     $                       TAUQ2(P+I) )
+               X22(P+I,Q+I) = ONE
+               CALL SLARF( 'L', M-P-Q-I+1, M-P-Q-I, X22(P+I,Q+I), 1,
+     $                     TAUQ2(P+I), X22(P+I,Q+I+1), LDX22, WORK )
+            END IF
 *
-            CALL SLARF( 'L', M-P-Q-I+1, M-P-Q-I, X22(P+I,Q+I), 1,
-     $                  TAUQ2(P+I), X22(P+I,Q+I+1), LDX22, WORK )
 *
          END DO
 *
index f9a47d6..ac0a7df 100644 (file)
          INFO = -8
       ELSE IF( Q .LT. 0 .OR. Q .GT. M ) THEN
          INFO = -9
-      ELSE IF( ( COLMAJOR .AND. LDX11 .LT. MAX(1,P) ) .OR.
-     $         ( .NOT.COLMAJOR .AND. LDX11 .LT. MAX(1,Q) ) ) THEN
-         INFO = -11
+      ELSE IF ( COLMAJOR .AND.  LDX11 .LT. MAX( 1, P ) ) THEN
+        INFO = -11
+      ELSE IF (.NOT. COLMAJOR .AND. LDX11 .LT. MAX( 1, Q ) ) THEN
+        INFO = -11
+      ELSE IF (COLMAJOR .AND. LDX12 .LT. MAX( 1, P ) ) THEN
+        INFO = -13
+      ELSE IF (.NOT. COLMAJOR .AND. LDX12 .LT. MAX( 1, M-Q ) ) THEN
+        INFO = -13
+      ELSE IF (COLMAJOR .AND. LDX21 .LT. MAX( 1, M-P ) ) THEN
+        INFO = -15
+      ELSE IF (.NOT. COLMAJOR .AND. LDX21 .LT. MAX( 1, Q ) ) THEN
+        INFO = -15
+      ELSE IF (COLMAJOR .AND. LDX22 .LT. MAX( 1, M-P ) ) THEN
+        INFO = -17
+      ELSE IF (.NOT. COLMAJOR .AND. LDX22 .LT. MAX( 1, M-Q ) ) THEN
+        INFO = -17
       ELSE IF( WANTU1 .AND. LDU1 .LT. P ) THEN
          INFO = -20
       ELSE IF( WANTU2 .AND. LDU2 .LT. M-P ) THEN
index 848507f..38cb9e3 100644 (file)
 *     ..
 *     .. Local Scalars ..
       LOGICAL            LEFT, RIGHT, TRAN, NOTRAN
-      INTEGER            I, IB, MB, LB, KF, Q
+      INTEGER            I, IB, MB, LB, KF, LDAQ, LDVQ
 *     ..
 *     .. External Functions ..
       LOGICAL            LSAME
       TRAN   = LSAME( TRANS, 'T' )
       NOTRAN = LSAME( TRANS, 'N' )
 *      
-      IF( LEFT ) THEN
-         Q = M
+      IF ( LEFT ) THEN
+         LDVQ = MAX( 1, M )
+         LDAQ = MAX( 1, K )
       ELSE IF ( RIGHT ) THEN
-         Q = N
+         LDVQ = MAX( 1, N )
+         LDAQ = MAX( 1, M )
       END IF
       IF( .NOT.LEFT .AND. .NOT.RIGHT ) THEN
          INFO = -1
          INFO = -5
       ELSE IF( L.LT.0 .OR. L.GT.K ) THEN
          INFO = -6         
-      ELSE IF( NB.LT.1 .OR. NB.GT.K ) THEN
+      ELSE IF( NB.LT.1 .OR. (NB.GT.K .AND. K.GT.0) ) THEN
          INFO = -7
-      ELSE IF( LDV.LT.MAX( 1, Q ) ) THEN
+      ELSE IF( LDV.LT.LDVQ ) THEN
          INFO = -9
       ELSE IF( LDT.LT.NB ) THEN
          INFO = -11
-      ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
+      ELSE IF( LDA.LT.LDAQ ) THEN
          INFO = -13
       ELSE IF( LDB.LT.MAX( 1, M ) ) THEN
          INFO = -15
index d9efb26..4cf2a02 100644 (file)
          INFO = -1
       ELSE IF( N.LT.0 ) THEN
          INFO = -2
-      ELSE IF( L.LT.0 .OR. L.GT.MIN(M,N) ) THEN
+      ELSE IF( L.LT.0 .OR. (L.GT.MIN(M,N) .AND. MIN(M,N).GT.0)) THEN
          INFO = -3
-      ELSE IF( NB.LT.1 .OR. NB.GT.N ) THEN
+      ELSE IF( NB.LT.1 .OR. (NB.GT.N .AND. N.GT.0)) THEN
          INFO = -4
       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
          INFO = -6
index 4a6b5c5..84090b1 100644 (file)
 *     Initial deflation
 *
       IMAX = Q
-      DO WHILE((IMAX.GT.1).AND.(PHI(MAX(IMAX-1,1)).EQ.ZERO))
+      DO WHILE( IMAX .GT. 1 )
+         IF( PHI(IMAX-1) .NE. ZERO ) THEN
+            EXIT
+         END IF
          IMAX = IMAX - 1
       END DO
       IMIN = IMAX - 1
index 71db24f..8f28ffc 100644 (file)
          INFO = -4
       ELSE IF( K.LT.0 .OR. K.GT.Q ) THEN
          INFO = -5
-      ELSE IF( NB.LT.1 .OR. NB.GT.K ) THEN
+      ELSE IF( NB.LT.1 .OR. (NB.GT.K .AND. K.GT.0)) THEN
          INFO = -6
       ELSE IF( LDV.LT.MAX( 1, Q ) ) THEN
          INFO = -8
index 684f6d9..b378ac8 100644 (file)
 *     ..
 *     .. Local Scalars ..
       LOGICAL            LEFT, RIGHT, TRAN, NOTRAN
-      INTEGER            I, IB, MB, LB, KF, Q
+      INTEGER            I, IB, MB, LB, KF, LDAQ, LDVQ
 *     ..
 *     .. External Functions ..
       LOGICAL            LSAME
       TRAN   = LSAME( TRANS, 'C' )
       NOTRAN = LSAME( TRANS, 'N' )
 *      
-      IF( LEFT ) THEN
-         Q = M
+      IF ( LEFT ) THEN
+         LDVQ = MAX( 1, M )
+         LDAQ = MAX( 1, K )
       ELSE IF ( RIGHT ) THEN
-         Q = N
+         LDVQ = MAX( 1, N )
+         LDAQ = MAX( 1, M )
       END IF
       IF( .NOT.LEFT .AND. .NOT.RIGHT ) THEN
          INFO = -1
          INFO = -5
       ELSE IF( L.LT.0 .OR. L.GT.K ) THEN
          INFO = -6         
-      ELSE IF( NB.LT.1 .OR. NB.GT.K ) THEN
+      ELSE IF( NB.LT.1 .OR. (NB.GT.K .AND. K.GT.0) ) THEN
          INFO = -7
-      ELSE IF( LDV.LT.MAX( 1, Q ) ) THEN
+      ELSE IF( LDV.LT.LDVQ ) THEN
          INFO = -9
       ELSE IF( LDT.LT.NB ) THEN
          INFO = -11
-      ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
+      ELSE IF( LDA.LT.LDAQ ) THEN
          INFO = -13
       ELSE IF( LDB.LT.MAX( 1, M ) ) THEN
          INFO = -15
index ac15f83..46067f5 100644 (file)
          INFO = -1
       ELSE IF( N.LT.0 ) THEN
          INFO = -2
-      ELSE IF( L.LT.0 .OR. L.GT.MIN(M,N) ) THEN
+      ELSE IF( L.LT.0 .OR. (L.GT.MIN(M,N) .AND. MIN(M,N).GT.0)) THEN
          INFO = -3
-      ELSE IF( NB.LT.1 .OR. NB.GT.N ) THEN
+      ELSE IF( NB.LT.1 .OR. (NB.GT.N .AND. N.GT.0)) THEN
          INFO = -4
       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
          INFO = -6
index 2455540..0bf2071 100644 (file)
 *     ..
 *     .. Local Scalars ..
       LOGICAL            COLMAJOR, LQUERY
-      INTEGER            I, LWORKMIN, LWORKOPT
+      INTEGER            I, LWORKMIN, LWORKOPT, PI1, QI1
       DOUBLE PRECISION   Z1, Z2, Z3, Z4
 *     ..
 *     .. External Subroutines ..
             THETA(I) = ATAN2( DZNRM2( M-P-I+1, X21(I,I), 1 ),
      $                 DZNRM2( P-I+1, X11(I,I), 1 ) )
 *
-            CALL ZLARFGP( P-I+1, X11(I,I), X11(I+1,I), 1, TAUP1(I) )
+            IF( P .GT. I ) THEN
+               CALL ZLARFGP( P-I+1, X11(I,I), X11(I+1,I), 1, TAUP1(I) )
+            ELSE IF ( P .EQ. I ) THEN
+               CALL ZLARFGP( P-I+1, X11(I,I), X11(I,I), 1, TAUP1(I) )
+            END IF
             X11(I,I) = ONE
-            CALL ZLARFGP( M-P-I+1, X21(I,I), X21(I+1,I), 1, TAUP2(I) )
+            IF ( M-P .GT. I ) THEN
+               CALL ZLARFGP( M-P-I+1, X21(I,I), X21(I+1,I), 1, 
+     $                       TAUP2(I) )
+            ELSE IF ( M-P .EQ. I ) THEN
+               CALL ZLARFGP( M-P-I+1, X21(I,I), X21(I,I), 1,
+     $                       TAUP2(I) )
+            END IF
             X21(I,I) = ONE
 *
-            CALL ZLARF( 'L', P-I+1, Q-I, X11(I,I), 1, DCONJG(TAUP1(I)),
-     $                  X11(I,I+1), LDX11, WORK )
-            CALL ZLARF( 'L', P-I+1, M-Q-I+1, X11(I,I), 1,
-     $                  DCONJG(TAUP1(I)), X12(I,I), LDX12, WORK )
-            CALL ZLARF( 'L', M-P-I+1, Q-I, X21(I,I), 1,
-     $                  DCONJG(TAUP2(I)), X21(I,I+1), LDX21, WORK )
-            CALL ZLARF( 'L', M-P-I+1, M-Q-I+1, X21(I,I), 1,
-     $                  DCONJG(TAUP2(I)), X22(I,I), LDX22, WORK )
+            IF ( Q .GT. I ) THEN
+               CALL ZLARF( 'L', P-I+1, Q-I, X11(I,I), 1, 
+     $                     DCONJG(TAUP1(I)), X11(I,I+1), LDX11, WORK )
+               CALL ZLARF( 'L', M-P-I+1, Q-I, X21(I,I), 1,
+     $                     DCONJG(TAUP2(I)), X21(I,I+1), LDX21, WORK )
+            END IF
+            IF ( M-Q+1 .GT. I ) THEN
+               CALL ZLARF( 'L', P-I+1, M-Q-I+1, X11(I,I), 1,
+     $                     DCONJG(TAUP1(I)), X12(I,I), LDX12, WORK )
+               CALL ZLARF( 'L', M-P-I+1, M-Q-I+1, X21(I,I), 1,
+     $                     DCONJG(TAUP2(I)), X22(I,I), LDX22, WORK )
+            END IF
 *
             IF( I .LT. Q ) THEN
                CALL ZSCAL( Q-I, DCMPLX( -Z1*Z3*SIN(THETA(I)), 0.0D0 ),
 *
             IF( I .LT. Q ) THEN
                CALL ZLACGV( Q-I, X11(I,I+1), LDX11 )
-               CALL ZLARFGP( Q-I, X11(I,I+1), X11(I,I+2), LDX11,
-     $                       TAUQ1(I) )
+               IF ( I .EQ. Q-1 ) THEN
+                  CALL ZLARFGP( Q-I, X11(I,I+1), X11(I,I+1), LDX11,
+     $                          TAUQ1(I) )
+               ELSE
+                  CALL ZLARFGP( Q-I, X11(I,I+1), X11(I,I+2), LDX11,
+     $                          TAUQ1(I) )
+               END IF
                X11(I,I+1) = ONE
             END IF
-            CALL ZLACGV( M-Q-I+1, X12(I,I), LDX12 )
-            CALL ZLARFGP( M-Q-I+1, X12(I,I), X12(I,I+1), LDX12,
-     $                    TAUQ2(I) )
+            IF ( M-Q+1 .GT. I ) THEN
+               CALL ZLACGV( M-Q-I+1, X12(I,I), LDX12 )
+               IF ( M-Q .EQ. I ) THEN
+                  CALL ZLARFGP( M-Q-I+1, X12(I,I), X12(I,I), LDX12,
+     $                          TAUQ2(I) )
+               ELSE
+                  CALL ZLARFGP( M-Q-I+1, X12(I,I), X12(I,I+1), LDX12,
+     $                          TAUQ2(I) )
+               END IF
+            END IF
             X12(I,I) = ONE
 *
             IF( I .LT. Q ) THEN
                CALL ZLARF( 'R', M-P-I, Q-I, X11(I,I+1), LDX11, TAUQ1(I),
      $                     X21(I+1,I+1), LDX21, WORK )
             END IF
-            CALL ZLARF( 'R', P-I, M-Q-I+1, X12(I,I), LDX12, TAUQ2(I),
-     $                  X12(I+1,I), LDX12, WORK )
-            CALL ZLARF( 'R', M-P-I, M-Q-I+1, X12(I,I), LDX12, TAUQ2(I),
-     $                  X22(I+1,I), LDX22, WORK )
+            IF ( P .GT. I ) THEN
+               CALL ZLARF( 'R', P-I, M-Q-I+1, X12(I,I), LDX12, TAUQ2(I),
+     $                     X12(I+1,I), LDX12, WORK )
+            END IF
+            IF ( M-P .GT. I ) THEN
+               CALL ZLARF( 'R', M-P-I, M-Q-I+1, X12(I,I), LDX12,
+     $                     TAUQ2(I), X22(I+1,I), LDX22, WORK )
+            END IF
 *
             IF( I .LT. Q )
      $         CALL ZLACGV( Q-I, X11(I,I+1), LDX11 )
             CALL ZSCAL( M-Q-I+1, DCMPLX( -Z1*Z4, 0.0D0 ), X12(I,I),
      $                  LDX12 )
             CALL ZLACGV( M-Q-I+1, X12(I,I), LDX12 )
-            CALL ZLARFGP( M-Q-I+1, X12(I,I), X12(I,I+1), LDX12,
-     $                    TAUQ2(I) )
+            IF ( I .GE. M-Q ) THEN
+               CALL ZLARFGP( M-Q-I+1, X12(I,I), X12(I,I), LDX12,
+     $                       TAUQ2(I) )
+            ELSE
+               CALL ZLARFGP( M-Q-I+1, X12(I,I), X12(I,I+1), LDX12,
+     $                       TAUQ2(I) )
+            END IF
             X12(I,I) = ONE
 *
-            CALL ZLARF( 'R', P-I, M-Q-I+1, X12(I,I), LDX12, TAUQ2(I),
-     $                  X12(I+1,I), LDX12, WORK )
+            IF ( P .GT. I ) THEN
+               CALL ZLARF( 'R', P-I, M-Q-I+1, X12(I,I), LDX12, TAUQ2(I),
+     $                     X12(I+1,I), LDX12, WORK )
+            END IF
             IF( M-P-Q .GE. 1 )
      $         CALL ZLARF( 'R', M-P-Q, M-Q-I+1, X12(I,I), LDX12,
      $                     TAUQ2(I), X22(Q+1,I), LDX22, WORK )
 *
             CALL ZLARFGP( P-I+1, X11(I,I), X11(I,I+1), LDX11, TAUP1(I) )
             X11(I,I) = ONE
-            CALL ZLARFGP( M-P-I+1, X21(I,I), X21(I,I+1), LDX21,
-     $                    TAUP2(I) )
+            IF ( I .EQ. M-P ) THEN
+               CALL ZLARFGP( M-P-I+1, X21(I,I), X21(I,I), LDX21,
+     $                       TAUP2(I) )
+            ELSE
+               CALL ZLARFGP( M-P-I+1, X21(I,I), X21(I,I+1), LDX21,
+     $                       TAUP2(I) )
+            END IF
             X21(I,I) = ONE
 *
             CALL ZLARF( 'R', Q-I, P-I+1, X11(I,I), LDX11, TAUP1(I),
             END IF
             CALL ZLARF( 'L', M-Q-I+1, P-I, X12(I,I), 1,
      $                  DCONJG(TAUQ2(I)), X12(I,I+1), LDX12, WORK )
-            CALL ZLARF( 'L', M-Q-I+1, M-P-I, X12(I,I), 1,
-     $                  DCONJG(TAUQ2(I)), X22(I,I+1), LDX22, WORK )
+            IF ( M-P .GT. I ) THEN
+               CALL ZLARF( 'L', M-Q-I+1, M-P-I, X12(I,I), 1,
+     $                     DCONJG(TAUQ2(I)), X22(I,I+1), LDX22, WORK )
+            END IF
 *
          END DO
 *
             CALL ZLARFGP( M-Q-I+1, X12(I,I), X12(I+1,I), 1, TAUQ2(I) )
             X12(I,I) = ONE
 *
-            CALL ZLARF( 'L', M-Q-I+1, P-I, X12(I,I), 1,
-     $                  DCONJG(TAUQ2(I)), X12(I,I+1), LDX12, WORK )
+            IF ( P .GT. I ) THEN
+               CALL ZLARF( 'L', M-Q-I+1, P-I, X12(I,I), 1,
+     $                     DCONJG(TAUQ2(I)), X12(I,I+1), LDX12, WORK )
+            END IF
             IF( M-P-Q .GE. 1 )
      $         CALL ZLARF( 'L', M-Q-I+1, M-P-Q, X12(I,I), 1,
      $                     DCONJG(TAUQ2(I)), X22(I,Q+1), LDX22, WORK )
      $                    TAUQ2(P+I) )
             X22(P+I,Q+I) = ONE
 *
-            CALL ZLARF( 'L', M-P-Q-I+1, M-P-Q-I, X22(P+I,Q+I), 1,
-     $                  DCONJG(TAUQ2(P+I)), X22(P+I,Q+I+1), LDX22,
-     $                  WORK )
+            IF ( M-P-Q .NE. I ) THEN
+               CALL ZLARF( 'L', M-P-Q-I+1, M-P-Q-I, X22(P+I,Q+I), 1,
+     $                     DCONJG(TAUQ2(P+I)), X22(P+I,Q+I+1), LDX22,
+     $                     WORK )
+            END IF
 *
          END DO
 *
index 3f7861d..356d9e9 100644 (file)
      $                   LBBCSDWORKOPT, LORBDBWORK, LORBDBWORKMIN,
      $                   LORBDBWORKOPT, LORGLQWORK, LORGLQWORKMIN,
      $                   LORGLQWORKOPT, LORGQRWORK, LORGQRWORKMIN,
-     $                   LORGQRWORKOPT, LWORKMIN, LWORKOPT
+     $                   LORGQRWORKOPT, LWORKMIN, LWORKOPT, P1, Q1
       LOGICAL            COLMAJOR, DEFAULTSIGNS, LQUERY, WANTU1, WANTU2,
      $                   WANTV1T, WANTV2T
       INTEGER            LRWORKMIN, LRWORKOPT
       EXTERNAL           LSAME
 *     ..
 *     .. Intrinsic Functions
-      INTRINSIC          COS, INT, MAX, MIN, SIN
+      INTRINSIC          INT, MAX, MIN
 *     ..
 *     .. Executable Statements ..
 *
          INFO = -8
       ELSE IF( Q .LT. 0 .OR. Q .GT. M ) THEN
          INFO = -9
-      ELSE IF( ( COLMAJOR .AND. LDX11 .LT. MAX(1,P) ) .OR.
-     $         ( .NOT.COLMAJOR .AND. LDX11 .LT. MAX(1,Q) ) ) THEN
-         INFO = -11
+      ELSE IF ( COLMAJOR .AND.  LDX11 .LT. MAX( 1, P ) ) THEN
+        INFO = -11
+      ELSE IF (.NOT. COLMAJOR .AND. LDX11 .LT. MAX( 1, Q ) ) THEN
+        INFO = -11
+      ELSE IF (COLMAJOR .AND. LDX12 .LT. MAX( 1, P ) ) THEN
+        INFO = -13
+      ELSE IF (.NOT. COLMAJOR .AND. LDX12 .LT. MAX( 1, M-Q ) ) THEN
+        INFO = -13
+      ELSE IF (COLMAJOR .AND. LDX21 .LT. MAX( 1, M-P ) ) THEN
+        INFO = -15
+      ELSE IF (.NOT. COLMAJOR .AND. LDX21 .LT. MAX( 1, Q ) ) THEN
+        INFO = -15
+      ELSE IF (COLMAJOR .AND. LDX22 .LT. MAX( 1, M-P ) ) THEN
+        INFO = -17
+      ELSE IF (.NOT. COLMAJOR .AND. LDX22 .LT. MAX( 1, M-Q ) ) THEN
+        INFO = -17
       ELSE IF( WANTU1 .AND. LDU1 .LT. P ) THEN
          INFO = -20
       ELSE IF( WANTU2 .AND. LDU2 .LT. M-P ) THEN
          IB22D = IB21E + MAX( 1, Q - 1 )
          IB22E = IB22D + MAX( 1, Q )
          IBBCSD = IB22E + MAX( 1, Q - 1 )
-         CALL ZBBCSD( JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS, M, P, Q, 0,
-     $                0, U1, LDU1, U2, LDU2, V1T, LDV1T, V2T, LDV2T, 0,
-     $                0, 0, 0, 0, 0, 0, 0, RWORK, -1, CHILDINFO )
+         CALL ZBBCSD( JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS, M, P, Q,
+     $                THETA, THETA, U1, LDU1, U2, LDU2, V1T, LDV1T,
+     $                V2T, LDV2T, THETA, THETA, THETA, THETA, THETA,
+     $                THETA, THETA, THETA, RWORK, -1, CHILDINFO )
          LBBCSDWORKOPT = INT( RWORK(1) )
          LBBCSDWORKMIN = LBBCSDWORKOPT
          LRWORKOPT = IBBCSD + LBBCSDWORKOPT - 1
          ITAUQ1 = ITAUP2 + MAX( 1, M - P )
          ITAUQ2 = ITAUQ1 + MAX( 1, Q )
          IORGQR = ITAUQ2 + MAX( 1, M - Q )
-         CALL ZUNGQR( M-Q, M-Q, M-Q, 0, MAX(1,M-Q), 0, WORK, -1,
+         CALL ZUNGQR( M-Q, M-Q, M-Q, U1, MAX(1,M-Q), U1, WORK, -1,
      $                CHILDINFO )
          LORGQRWORKOPT = INT( WORK(1) )
          LORGQRWORKMIN = MAX( 1, M - Q )
          IORGLQ = ITAUQ2 + MAX( 1, M - Q )
-         CALL ZUNGLQ( M-Q, M-Q, M-Q, 0, MAX(1,M-Q), 0, WORK, -1,
+         CALL ZUNGLQ( M-Q, M-Q, M-Q, U1, MAX(1,M-Q), U1, WORK, -1,
      $                CHILDINFO )
          LORGLQWORKOPT = INT( WORK(1) )
          LORGLQWORKMIN = MAX( 1, M - Q )
          IORBDB = ITAUQ2 + MAX( 1, M - Q )
          CALL ZUNBDB( TRANS, SIGNS, M, P, Q, X11, LDX11, X12, LDX12,
-     $                X21, LDX21, X22, LDX22, 0, 0, 0, 0, 0, 0, WORK,
-     $                -1, CHILDINFO )
+     $                X21, LDX21, X22, LDX22, THETA, THETA, U1, U2,
+     $                V1T, V2T, WORK, -1, CHILDINFO )
          LORBDBWORKOPT = INT( WORK(1) )
          LORBDBWORKMIN = LORBDBWORKOPT
          LWORKOPT = MAX( IORGQR + LORGQRWORKOPT, IORGLQ + LORGLQWORKOPT,
          END IF
          IF( WANTV2T .AND. M-Q .GT. 0 ) THEN
             CALL ZLACPY( 'U', P, M-Q, X12, LDX12, V2T, LDV2T )
-            CALL ZLACPY( 'U', M-P-Q, M-P-Q, X22(Q+1,P+1), LDX22,
-     $                   V2T(P+1,P+1), LDV2T )
-            CALL ZUNGLQ( M-Q, M-Q, M-Q, V2T, LDV2T, WORK(ITAUQ2),
-     $                   WORK(IORGLQ), LORGLQWORK, INFO )
+            IF( M-P .GT. Q) THEN
+               CALL ZLACPY( 'U', M-P-Q, M-P-Q, X22(Q+1,P+1), LDX22,
+     $                      V2T(P+1,P+1), LDV2T )
+            END IF
+            IF( M .GT. Q ) THEN
+               CALL ZUNGLQ( M-Q, M-Q, M-Q, V2T, LDV2T, WORK(ITAUQ2),
+     $                      WORK(IORGLQ), LORGLQWORK, INFO )
+            END IF
          END IF
       ELSE
          IF( WANTU1 .AND. P .GT. 0 ) THEN
      $                   WORK(IORGQR), LORGQRWORK, INFO )
          END IF
          IF( WANTV2T .AND. M-Q .GT. 0 ) THEN
+            P1 = MIN( P+1, M )
+            Q1 = MIN( Q+1, M )
             CALL ZLACPY( 'L', M-Q, P, X12, LDX12, V2T, LDV2T )
-            CALL ZLACPY( 'L', M-P-Q, M-P-Q, X22(P+1,Q+1), LDX22,
-     $                   V2T(P+1,P+1), LDV2T )
+            IF( M .GT. P+Q ) THEN
+               CALL ZLACPY( 'L', M-P-Q, M-P-Q, X22(P1,Q1), LDX22,
+     $                      V2T(P+1,P+1), LDV2T )
+            END IF
             CALL ZUNGQR( M-Q, M-Q, M-Q, V2T, LDV2T, WORK(ITAUQ2),
      $                   WORK(IORGQR), LORGQRWORK, INFO )
          END IF