Skip to content

Instantly share code, notes, and snippets.

@pv

pv/arpack-osx-failure.f

Last active December 19, 2015 00:59
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save pv/5873146 to your computer and use it in GitHub Desktop.
Save pv/5873146 to your computer and use it in GitHub Desktop.
Download ZIP
An ARPACK-using program that fails on OSX if linked with Accelerate.
Raw
arpack-osx-failure.f
c ==============================
c Arpack failure with Accelerate
c ==============================
c
c Example of an ARPACK using program that fails when linked with
c Accelerate on OSX.
c
c This computes two eigenvalues of a 6x6 matrix (a randomly chosen
c one), via shift-invert with sigma=0.5. On OSX 10.8.4 (and probably
c others), the eigenvector computed for the second eigenvalue is false.
c
c
c Build (Accelerate)
c ==================
c
c 1. Download Arpack-NG 3.1.3 source code:
c http://forge.scilab.org/index.php/p/arpack-ng/downloads/get/arpack-ng-3.1.3.tar.gz
c
c 2. Unpack, and save this file as failure.f under arpack-ng-3.1.3/
c
c 3. rm UTIL/second_NONE.f
c
c 4. gfortran -m64 -ff2c -o failure failure.f SRC/*.f UTIL/*.f -Wl,-framework -Wl,Accelerate
c
c 5. ./failure
c
c
c Result (Accelerate)
c ===================
c
c Col 1 Col 2
c Row 1: 4.53230D-01 3.95137D-15
c Row 2: 4.58707D-01 3.72325D+10
c
c Note that the second eigenvector is miscomputed (huge residual).
c
c
c Result (reference BLAS)
c =======================
c
c Col 1 Col 2
c Row 1: 4.53230D-01 3.49843D-15
c Row 2: 4.58707D-01 2.42642D-15
c
c
c The main routine is copypasted from EXAMPLES/SYM/dsdrv2.f
c
program failure
integer maxn, maxnev, maxncv, ldv
parameter (maxn=6, maxnev=5, maxncv=6,
& ldv=maxn )
c
c %--------------%
c | Local Arrays |
c %--------------%
c
Double precision
& v(ldv,maxncv), workl(maxncv*(maxncv+8)),
& workd(3*maxn), d(maxncv,2), resid(maxn),
& ad(maxn), adl(maxn), adu(maxn), adu2(maxn),
& ax(maxn)
logical select(maxncv)
integer iparam(11), ipntr(11), ipiv(maxn)
c
c %---------------%
c | Local Scalars |
c %---------------%
c
character bmat*1, which*2
integer ido, n, nev, ncv, lworkl, info, j, ierr,
& nconv, maxitr, ishfts, mode
logical rvec
Double precision
& sigma, tol, h2
c
c %------------%
c | Parameters |
c %------------%
c
Double precision
& zero, one, two
parameter (zero = 0.0D+0, one = 1.0D+0,
& two = 2.0D+0)
c
c %-----------------------------%
c | BLAS & LAPACK routines used |
c %-----------------------------%
c
Double precision
& dnrm2
external daxpy, dnrm2, dgttrf, dgttrs
c
c %--------------------%
c | Intrinsic function |
c %--------------------%
c
intrinsic abs
c
c %-----------------------%
c | Executable Statements |
c %-----------------------%
c
c %----------------------------------------------------%
c | The number N is the dimension of the matrix. A |
c | standard eigenvalue problem is solved (BMAT = 'I'. |
c | NEV is the number of eigenvalues (closest to |
c | SIGMA) to be approximated. Since the shift-invert |
c | mode is used, WHICH is set to 'LM'. The user can |
c | modify NEV, NCV, SIGMA to solve problems of |
c | different sizes, and to get different parts of the |
c | spectrum. However, The following conditions must |
c | be satisfied: |
c | N <= MAXN, |
c | NEV <= MAXNEV, |
c | NEV + 1 <= NCV <= MAXNCV |
c %----------------------------------------------------%
c
n = 6
nev = 2
ncv = 5
if ( n .gt. maxn ) then
print *, ' ERROR with _SDRV2: N is greater than MAXN '
go to 9000
else if ( nev .gt. maxnev ) then
print *, ' ERROR with _SDRV2: NEV is greater than MAXNEV '
go to 9000
else if ( ncv .gt. maxncv ) then
print *, ' ERROR with _SDRV2: NCV is greater than MAXNCV '
go to 9000
end if
c
bmat = 'I'
which = 'LM'
sigma = 0.5
c
c %--------------------------------------------------%
c | The work array WORKL is used in DSAUPD as |
c | workspace. Its dimension LWORKL is set as |
c | illustrated below. The parameter TOL determines |
c | the stopping criterion. If TOL<=0, machine |
c | precision is used. The variable IDO is used for |
c | reverse communication and is initially set to 0. |
c | Setting INFO=0 indicates that a random vector is |
c | generated in DSAUPD to start the Arnoldi |
c | iteration. |
c %--------------------------------------------------%
c
lworkl = ncv*(ncv+8)
tol = zero
ido = 0
info = 0
c
c %---------------------------------------------------%
c | This program uses exact shifts with respect to |
c | the current Hessenberg matrix (IPARAM(1) = 1). |
c | IPARAM(3) specifies the maximum number of Arnoldi |
c | iterations allowed. Mode 3 of DSAUPD is used |
c | (IPARAM(7) = 3). All these options may be |
c | changed by the user. For details, see the |
c | documentation in DSAUPD. |
c %---------------------------------------------------%
c
ishfts = 1
maxitr = 300
mode = 3
c
iparam(1) = ishfts
iparam(3) = maxitr
iparam(7) = mode
c
c %-------------------------------------------%
c | M A I N L O O P (Reverse communication) |
c %-------------------------------------------%
c
10 continue
c
c %---------------------------------------------%
c | Repeatedly call the routine DSAUPD and take |
c | actions indicated by parameter IDO until |
c | either convergence is indicated or maxitr |
c | has been exceeded. |
c %---------------------------------------------%
c
call dsaupd ( ido, bmat, n, which, nev, tol, resid,
& ncv, v, ldv, iparam, ipntr, workd, workl,
& lworkl, info )
c
if (ido .eq. -1 .or. ido .eq. 1) then
c
c %----------------------------------------%
c | Perform y <-- OP*x = inv[A-sigma*I]*x. |
c | The user only need the linear system |
c | solver here that takes workd(ipntr(1)) |
c | as input, and returns the result to |
c | workd(ipntr(2)). |
c %----------------------------------------%
c
call dcopy (n, workd(ipntr(1)), 1, workd(ipntr(2)), 1)
c
call tis(sigma, workd(ipntr(2)))
c
c %-----------------------------------------%
c | L O O P B A C K to call DSAUPD again. |
c %-----------------------------------------%
c
go to 10
c
end if
c
c %----------------------------------------%
c | Either we have convergence or there is |
c | an error. |
c %----------------------------------------%
c
if ( info .lt. 0 ) then
c
c %----------------------------%
c | Error message. Check the |
c | documentation in DSAUPD |
c %----------------------------%
c
print *, ' '
print *, ' Error with _saupd, info = ',info
print *, ' Check documentation of _saupd '
print *, ' '
c
else
c
c %-------------------------------------------%
c | No fatal errors occurred. |
c | Post-Process using DSEUPD. |
c | |
c | Computed eigenvalues may be extracted. |
c | |
c | Eigenvectors may also be computed now if |
c | desired. (indicated by rvec = .true.) |
c %-------------------------------------------%
c
rvec = .true.
c
call dseupd ( rvec, 'All', select, d, v, ldv, sigma,
& bmat, n, which, nev, tol, resid, ncv, v, ldv,
& iparam, ipntr, workd, workl, lworkl, ierr )
c
c %----------------------------------------------%
c | Eigenvalues are returned in the first column |
c | of the two dimensional array D and the |
c | corresponding eigenvectors are returned in |
c | the first NEV columns of the two dimensional |
c | array V if requested. Otherwise, an |
c | orthogonal basis for the invariant subspace |
c | corresponding to the eigenvalues in D is |
c | returned in V. |
c %----------------------------------------------%
if ( ierr .ne. 0 ) then
c
c %------------------------------------%
c | Error condition: |
c | Check the documentation of DSEUPD. |
c %------------------------------------%
c
print *, ' '
print *, ' Error with _seupd, info = ', ierr
print *, ' Check the documentation of _seupd '
print *, ' '
c
else
c
nconv = iparam(5)
do 30 j=1, nconv
c
c %---------------------------%
c | Compute the residual norm |
c | |
c | || A*x - lambda*x || |
c | |
c | for the NCONV accurately |
c | computed eigenvalues and |
c | eigenvectors. (iparam(5) |
c | indicates how many are |
c | accurate to the requested |
c | tolerance) |
c %---------------------------%
c
call tv(v(1,j), ax)
call daxpy(n, -d(j,1), v(1,j), 1, ax, 1)
d(j,2) = dnrm2(n, ax, 1)
d(j,2) = d(j,2) / abs(d(j,1))
c
30 continue
c
c %-------------------------------%
c | Display computed residuals |
c %-------------------------------%
c
call dmout(6, nconv, 2, d, maxncv, -6,
& 'Ritz values and relative residuals')
end if
c
c %------------------------------------------%
c | Print additional convergence information |
c %------------------------------------------%
c
if ( info .eq. 1) then
print *, ' '
print *, ' Maximum number of iterations reached.'
print *, ' '
else if ( info .eq. 3) then
print *, ' '
print *, ' No shifts could be applied during implicit',
& ' Arnoldi update, try increasing NCV.'
print *, ' '
end if
c
print *, ' '
print *, ' _SDRV2 '
print *, ' ====== '
print *, ' '
print *, ' Size of the matrix is ', n
print *, ' The number of Ritz values requested is ', nev
print *, ' The number of Arnoldi vectors generated',
& ' (NCV) is ', ncv
print *, ' What portion of the spectrum: ', which
print *, ' The number of converged Ritz values is ',
& nconv
print *, ' The number of Implicit Arnoldi update',
& ' iterations taken is ', iparam(3)
print *, ' The number of OP*x is ', iparam(9)
print *, ' The convergence criterion is ', tol
print *, ' '
c
end if
c
c %---------------------------%
c | Done with program dsdrv2. |
c %---------------------------%
c
9000 continue
c
end
c-------------------------------------------------------------------
c
c This is the problematic matrix:
c
subroutine gett(T)
Double precision T(6,6)
T(1,1) = 1.717005133628845215d+00
T(1,2) = 1.300954699516296387d+00
T(1,3) = 9.104259014129638672d-01
T(1,4) = 9.700312614440917969d-01
T(1,5) = 1.463213443756103516d+00
T(1,6) = 1.825662374496459961d+00
T(2,1) = 1.300954699516296387d+00
T(2,2) = 2.253338098526000977d+00
T(2,3) = 1.517550826072692871d+00
T(2,4) = 1.559799194335937500d+00
T(2,5) = 2.288973808288574219d+00
T(2,6) = 2.074316501617431641d+00
T(3,1) = 9.104259014129638672d-01
T(3,2) = 1.517550826072692871d+00
T(3,3) = 1.277888774871826172d+00
T(3,4) = 1.135284543037414551d+00
T(3,5) = 1.688696026802062988d+00
T(3,6) = 1.150150656700134277d+00
T(4,1) = 9.700312614440917969d-01
T(4,2) = 1.559799194335937500d+00
T(4,3) = 1.135284543037414551d+00
T(4,4) = 1.538701534271240234d+00
T(4,5) = 1.518137097358703613d+00
T(4,6) = 1.316925644874572754d+00
T(5,1) = 1.463213443756103516d+00
T(5,2) = 2.288973808288574219d+00
T(5,3) = 1.688696026802062988d+00
T(5,4) = 1.518137097358703613d+00
T(5,5) = 2.559108495712280273d+00
T(5,6) = 1.976615428924560547d+00
T(6,1) = 1.825662374496459961d+00
T(6,2) = 2.074316501617431641d+00
T(6,3) = 1.150150656700134277d+00
T(6,4) = 1.316925644874572754d+00
T(6,5) = 1.976615428924560547d+00
T(6,6) = 2.817672491073608398d+00
return
end
subroutine tis(sigma, v)
integer i, j, info, ipiv(6)
double precision sigma, v(6)
double precision T(6,6)
call gett(T)
do 10 i = 1, 6
T(i,i) = T(i,i) - sigma
10 continue
call dgesv(6, 1, T, 6, ipiv, v, 6, info)
if (info.ne.0) then
write(*,*) 'dgesv failed'
stop
end if
return
end
subroutine tv (x, y)
double precision x(6), y(6)
Double precision T(6,6)
call gett(T)
call dgemv('N', 6, 6, 1d0, T, 6, x, 1, 0d0, y, 1)
return
end
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment