[SciPy-Dev] calc_lwork.f different from Lapack?

Josh Lawrence josh.k.lawrence@gmail....
Wed Apr 13 14:10:59 CDT 2011


Hello,

After receiving an error in using linalg.lstsq, I started some 
investigating. I found that scipy/linalg/src/calc_lwork.f is essentially 
a stripped down version of the appropriate functions in Lapack (for 
linalg.lstsq using double precision complex numbers, it is ZGELSS). It 
seems that the portions for gelss in calc_lwork.f do not match the 
current version of lapack. I believe this is the source of my error (it 
doesn't occur if I use linalg.pinv2 for least squares instead). Is there 
a particular rationale for these being different or perhaps an incorrect 
copy and paste? The relevant portions are shown below (a line of % 
characters are prepended before the codes from the two locations).

For example, the first calculation of MAXWRK under Path 1 uses 3*N in 
calc_lwork and 2*N in lapack. If things were just being over-estimated, 
I wouldn't really care, but I was able to isolate down that the variable 
LWORK (calculated from MINWRK) is too small, which if I read the rest of 
the code right means that some work arrays are not long enough.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
calc_lwork.f (from scipy trunk on github, lines 191-261):

       IF( M.GE.N .AND. M.GE.MNTHR ) THEN
*
*     Path 1a - overdetermined, with many more rows than columns
*
          MM = N
          MAXWRK = MAX( MAXWRK, N+N*ILAENV( 1, prefix //  'GEQRF', ' ',
      $        M, N, -1, -1 ) )
          MAXWRK = MAX( MAXWRK, N+NRHS*
      $        ILAENV( 1, prefix //  'ORMQR', 'LT', M, NRHS, N, -1 ) )
       END IF
       IF( M.GE.N ) THEN
*
*     Path 1 - overdetermined or exactly determined
*
*     Compute workspace neede for BDSQR
*
          BDSPAC = MAX( 1, 5*N )
          MAXWRK = MAX( MAXWRK, 3*N+( MM+N )*
      $        ILAENV( 1, prefix // 'GEBRD', ' ', MM, N, -1, -1 ) )
          MAXWRK = MAX( MAXWRK, 3*N+NRHS*
      $        ILAENV( 1, prefix // 'ORMBR', 'QLT', MM, NRHS, N, -1 ) )
          MAXWRK = MAX( MAXWRK, 3*N+( N-1 )*
      $        ILAENV( 1, prefix // 'ORGBR', 'P', N, N, N, -1 ) )
          MAXWRK = MAX( MAXWRK, BDSPAC )
          MAXWRK = MAX( MAXWRK, N*NRHS )
          MINWRK = MAX( 3*N+MM, 3*N+NRHS, BDSPAC )
          MAXWRK = MAX( MINWRK, MAXWRK )
       END IF

       IF( N.GT.M ) THEN
*
*     Compute workspace neede for DBDSQR
*
          BDSPAC = MAX( 1, 5*M )
          MINWRK = MAX( 3*M+NRHS, 3*M+N, BDSPAC )
          IF( N.GE.MNTHR ) THEN
*
*     Path 2a - underdetermined, with many more columns
*     than rows
*
             MAXWRK = M + M*ILAENV( 1, prefix // 'GELQF', ' ',
      $           M, N, -1, -1 )
             MAXWRK = MAX( MAXWRK, M*M+4*M+2*M*
      $           ILAENV( 1, prefix //  'GEBRD', ' ', M, M, -1, -1 ) )
             MAXWRK = MAX( MAXWRK, M*M+4*M+NRHS*
      $           ILAENV( 1, prefix // 'ORMBR', 'QLT', M, NRHS, M, -1 ))
             MAXWRK = MAX( MAXWRK, M*M+4*M+( M-1 )*
      $           ILAENV( 1, prefix // 'ORGBR', 'P', M, M, M, -1 ) )
             MAXWRK = MAX( MAXWRK, M*M+M+BDSPAC )
             IF( NRHS.GT.1 ) THEN
                MAXWRK = MAX( MAXWRK, M*M+M+M*NRHS )
             ELSE
                MAXWRK = MAX( MAXWRK, M*M+2*M )
             END IF
             MAXWRK = MAX( MAXWRK, M+NRHS*
      $           ILAENV( 1, prefix // 'ORMLQ', 'LT', N, NRHS, M, -1 ) )

          ELSE
*
*     Path 2 - underdetermined
*
             MAXWRK = 3*M + ( N+M )*ILAENV( 1, prefix // 'GEBRD', ' ',
      $           M, N, -1, -1 )
             MAXWRK = MAX( MAXWRK, 3*M+NRHS*
      $           ILAENV( 1, prefix // 'ORMBR', 'QLT', M, NRHS, M, -1 ) )
             MAXWRK = MAX( MAXWRK, 3*M+M*
      $           ILAENV( 1, prefix // 'ORGBR', 'P', M, N, M, -1 ) )
             MAXWRK = MAX( MAXWRK, BDSPAC )
             MAXWRK = MAX( MAXWRK, N*NRHS )
          END IF
       END IF

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ZGELSS.F (from LAPACK on netlib, lines 165-236):

       IF( INFO.EQ.0 ) THEN
          MINWRK = 1
          MAXWRK = 1
          IF( MINMN.GT.0 ) THEN
             MM = M
             MNTHR = ILAENV( 6, 'ZGELSS', ' ', M, N, NRHS, -1 )
             IF( M.GE.N .AND. M.GE.MNTHR ) THEN
*
*              Path 1a - overdetermined, with many more rows than
*                        columns
*
                MM = N
                MAXWRK = MAX( MAXWRK, N + N*ILAENV( 1, 'ZGEQRF', ' ', M,
      $                       N, -1, -1 ) )
                MAXWRK = MAX( MAXWRK, N + NRHS*ILAENV( 1, 'ZUNMQR', 'LC',
      $                       M, NRHS, N, -1 ) )
             END IF
             IF( M.GE.N ) THEN
*
*              Path 1 - overdetermined or exactly determined
*
                MAXWRK = MAX( MAXWRK, 2*N + ( MM + N )*ILAENV( 1,
      $                       'ZGEBRD', ' ', MM, N, -1, -1 ) )
                MAXWRK = MAX( MAXWRK, 2*N + NRHS*ILAENV( 1, 'ZUNMBR',
      $                       'QLC', MM, NRHS, N, -1 ) )
                MAXWRK = MAX( MAXWRK, 2*N + ( N - 1 )*ILAENV( 1,
      $                       'ZUNGBR', 'P', N, N, N, -1 ) )
                MAXWRK = MAX( MAXWRK, N*NRHS )
                MINWRK = 2*N + MAX( NRHS, M )
             END IF
             IF( N.GT.M ) THEN
                MINWRK = 2*M + MAX( NRHS, N )
                IF( N.GE.MNTHR ) THEN
*
*                 Path 2a - underdetermined, with many more columns
*                 than rows
*
                   MAXWRK = M + M*ILAENV( 1, 'ZGELQF', ' ', M, N, -1,
      $                     -1 )
                   MAXWRK = MAX( MAXWRK, 3*M + M*M + 2*M*ILAENV( 1,
      $                          'ZGEBRD', ' ', M, M, -1, -1 ) )
                   MAXWRK = MAX( MAXWRK, 3*M + M*M + NRHS*ILAENV( 1,
      $                          'ZUNMBR', 'QLC', M, NRHS, M, -1 ) )
                   MAXWRK = MAX( MAXWRK, 3*M + M*M + ( M - 1 )*ILAENV( 1,
      $                          'ZUNGBR', 'P', M, M, M, -1 ) )
                   IF( NRHS.GT.1 ) THEN
                      MAXWRK = MAX( MAXWRK, M*M + M + M*NRHS )
                   ELSE
                      MAXWRK = MAX( MAXWRK, M*M + 2*M )
                   END IF
                   MAXWRK = MAX( MAXWRK, M + NRHS*ILAENV( 1, 'ZUNMLQ',
      $                          'LC', N, NRHS, M, -1 ) )
                ELSE
*
*                 Path 2 - underdetermined
*
                   MAXWRK = 2*M + ( N + M )*ILAENV( 1, 'ZGEBRD', ' ', M,
      $                     N, -1, -1 )
                   MAXWRK = MAX( MAXWRK, 2*M + NRHS*ILAENV( 1, 'ZUNMBR',
      $                          'QLC', M, NRHS, M, -1 ) )
                   MAXWRK = MAX( MAXWRK, 2*M + M*ILAENV( 1, 'ZUNGBR',
      $                          'P', M, N, M, -1 ) )
                   MAXWRK = MAX( MAXWRK, N*NRHS )
                END IF
             END IF
             MAXWRK = MAX( MINWRK, MAXWRK )
          END IF
          WORK( 1 ) = MAXWRK
*
          IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY )
      $      INFO = -12
       END IF




More information about the SciPy-Dev mailing list