Merge pull request #2684 from pv/arpack-fix

BUG: arpack-ng: fix bug occurring on OSX or LAPACK 3.2.1
scipy · Aug 6, 2013 · 862eb08 · 862eb08
2 parents aa8ec9b + 54073ef
commit 862eb08
Show file tree

Hide file tree

Showing 2 changed files with 24 additions and 4 deletions.
diff --git a/scipy/sparse/linalg/eigen/arpack/ARPACK/SRC/dseupd.f b/scipy/sparse/linalg/eigen/arpack/ARPACK/SRC/dseupd.f
@@ -760,6 +760,16 @@ subroutine dseupd (rvec  , howmny, select, d    ,
      &                ldq   , workl(iw+ncv), workl(ihb),
      &                ncv   , temp         , ierr)
 c
+c        %-----------------------------------------------------%
+c        | Make a copy of the last row into                    |
+c        | workl(iw+ncv:iw+2*ncv), as it is needed again in    |
+c        | the Ritz vector purification step below             |
+c        %-----------------------------------------------------%
+c
+         do 67 j = 1, nconv
+            workl(iw+ncv+j-1) = workl(ihb+j-1)
+ 67      continue
+
       else if (rvec .and. howmny .eq. 'S') then
 c
 c     Not yet implemented. See remark 2 above.
@@ -830,14 +840,14 @@ subroutine dseupd (rvec  , howmny, select, d    ,
       if (rvec .and. (type .eq. 'SHIFTI' .or. type .eq. 'CAYLEY')) then
 c
          do 110 k=0, nconv-1
-            workl(iw+k) = workl(iq+k*ldq+ncv-1)
+            workl(iw+k) = workl(iw+ncv+k)
      &                  / workl(iw+k)
  110     continue
 c
       else if (rvec .and. type .eq. 'BUCKLE') then
 c
          do 120 k=0, nconv-1
-            workl(iw+k) = workl(iq+k*ldq+ncv-1)
+            workl(iw+k) = workl(iw+ncv+k)
      &                  / (workl(iw+k)-one)
  120     continue
 c

diff --git a/scipy/sparse/linalg/eigen/arpack/ARPACK/SRC/sseupd.f b/scipy/sparse/linalg/eigen/arpack/ARPACK/SRC/sseupd.f
@@ -760,6 +760,16 @@ subroutine sseupd(rvec  , howmny, select, d    ,
      &                ldq   , workl(iw+ncv), workl(ihb),
      &                ncv   , temp         , ierr)
 c
+c        %-----------------------------------------------------%
+c        | Make a copy of the last row into                    |
+c        | workl(iw+ncv:iw+2*ncv), as it is needed again in    |
+c        | the Ritz vector purification step below             |
+c        %-----------------------------------------------------%
+c
+         do 67 j = 1, nconv
+            workl(iw+ncv+j-1) = workl(ihb+j-1)
+ 67      continue
+
       else if (rvec .and. howmny .eq. 'S') then
 c
 c     Not yet implemented. See remark 2 above.
@@ -830,14 +840,14 @@ subroutine sseupd(rvec  , howmny, select, d    ,
       if (rvec .and. (type .eq. 'SHIFTI' .or. type .eq. 'CAYLEY')) then
 c
          do 110 k=0, nconv-1
-            workl(iw+k) = workl(iq+k*ldq+ncv-1)
+            workl(iw+k) = workl(iw+ncv+k)
      &                  / workl(iw+k)
  110     continue
 c
       else if (rvec .and. type .eq. 'BUCKLE') then
 c
          do 120 k=0, nconv-1
-            workl(iw+k) = workl(iq+k*ldq+ncv-1)
+            workl(iw+k) = workl(iw+ncv+k)
      &                  / (workl(iw+k)-one)
  120     continue
 c