Compiling and using ARPACK on Mac OS X
Alec Jacobson
November 27, 2010
Computer programming paleontology
Recently we have been prototyping using MATLAB's eigs. It is an extremely easy to use eigen-decomposition tool that works on sparse matrices. It allows you to choose how many eigenvalues you want from which end of the spectrum (smallest or largest magnitude). Here's an example of how easy it is:
% Build a sparse 100 by 100 second-order finite-difference Laplacian matrix
A = delsq(numgrid('C',10));
% Get 5 smallest magnitude eigenvectors (columns of V) and eigenvalues
% (diagonal of D).
[V,D] = eigs(A,5,'sm')
After some digging around for a C/C++ equivalent I found that MATLAB is interfacing ARPACK. For some reason ARPACK feels ancient, but it's actually been around since only 1996 and seemed to be maintained at least as recently as 2008. Trying to compile it though felt like what it must have felt like for Otto Lidenbrock to find dinosaurs living at present day (only near the center of the earth).
Compiling ARPACK
I loosely followed these helpful instructions. This also assumes that you have a universal (32-bit and 64-bit) build of f2c, if not I have previously posted instructions.
Download the ARPACK source and patch as instructed on the ARPACK download site.
Use zcat to patch the directories. If Safari auto uncompresses the .zips you may have to start over with pure zips to be sure that it's being done as instructed here.
Edit the ARmake.inc file, to look like this:
home = .
BLASdir = $(home)/BLAS
LAPACKdir = $(home)/LAPACK
UTILdir = $(home)/UTIL
SRCdir = $(home)/SRC
DIRS = $(UTILdir) $(SRCdir)
ARPACKLIB = $(home)/libarpack.a
.SUFFIXES: .f .o
.DEFAULT:
@$(ECHO) "Unknown target $@, try: make help"
F2C = f2c
F2CFLAGS = -ARw8 -Nn802 -Nq300 -Nx400
CC = cc
CFLAGS = -arch i386 -arch x86_64 -O
CD = cd
ECHO = echo
MAKE = /usr/bin/make
RM = rm
RMFLAGS = -f
help:
@$(ECHO) "usage: make ?"
Edit Makefile, to look like this:
include ARmake.inc
PRECISIONS = single double complex complex16
all: lib
lib: arpacklib
clean: cleanlib
arpacklib:
@( \
for f in $(DIRS); \
do \
$(CD) $$f; \
$(ECHO) Making lib in $$f; \
$(MAKE) $(PRECISIONS); \
$(CD) ..; \
done );
cleanlib:
( cd $(BLASdir); $(MAKE) clean )
( cd $(LAPACKdir); $(MAKE) clean )
( cd $(UTILdir); $(MAKE) clean )
( cd $(SRCdir); $(MAKE) clean )
help:
@$(ECHO) "usage: make ?"
Edit SRC/Makefile to look like this:
include ../ARmake.inc
SOBJ = sgetv0.o slaqrb.o sstqrb.o ssortc.o ssortr.o sstatn.o sstats.o \
snaitr.o snapps.o snaup2.o snaupd.o snconv.o sneigh.o sngets.o \
ssaitr.o ssapps.o ssaup2.o ssaupd.o ssconv.o sseigt.o ssgets.o \
sneupd.o sseupd.o ssesrt.o
DOBJ = dgetv0.o dlaqrb.o dstqrb.o dsortc.o dsortr.o dstatn.o dstats.o \
dnaitr.o dnapps.o dnaup2.o dnaupd.o dnconv.o dneigh.o dngets.o \
dsaitr.o dsapps.o dsaup2.o dsaupd.o dsconv.o dseigt.o dsgets.o \
dneupd.o dseupd.o dsesrt.o
COBJ = cnaitr.o cnapps.o cnaup2.o cnaupd.o cneigh.o cneupd.o cngets.o \
cgetv0.o csortc.o cstatn.o
ZOBJ = znaitr.o znapps.o znaup2.o znaupd.o zneigh.o zneupd.o zngets.o \
zgetv0.o zsortc.o zstatn.o
.f.o:
$(F2C) $(F2CFLAGS) $*.f
$(CC) $(CFLAGS) -c $*.c
$(RM) $*.c
all: single double complex complex16
single: $(SOBJ)
double: $(DOBJ)
complex: $(COBJ)
complex16: $(ZOBJ)
#
# clean - remove all object files
#
clean:
rm -f *.o a.out core *.c
Edit UTIL/Makefile to look like this:
include ../ARmake.inc
OBJS = icnteq.o icopy.o iset.o iswap.o ivout.o second.o
SOBJ = svout.o smout.o
DOBJ = dvout.o dmout.o
COBJ = cvout.o cmout.o
ZOBJ = zvout.o zmout.o
.SUFFIXES: .o .F .f
.f.o:
$(F2C) $(F2CFLAGS) $*.f
$(CC) $(CFLAGS) -c $*.c
$(RM) $*.c
#
# make the library containing both single and double precision
#
all: single double complex complex16
single: $(SOBJ) $(OBJS)
double: $(DOBJ) $(OBJS) $(ZOBJ)
complex: $(SOBJ) $(OBJS) $(COBJ)
complex16: $(DOBJ) $(OBJS) $(ZOBJ)
#
# clean - remove all object files
#
clean:
rm -f *.o a.out core *.c
Now you can build all of the .o object files with:
make lib
So far we have not assembled the static library, to do this issue:
libtool -o libarpack.a SRC/*.o UTIL/*.o
Using ARPACK
Travel to EXAMPLES/SIMPLE
Then you can convert the fortran file to a c file with:
f2c sssimp.f
Then compile with:
gcc -arch i386 -arch x86_64 -framework Accelerate -lf2c ../../libarpack.a sssimp.c -o sssimp
And finally run with:
./sssimp
You should see something like:
_saupd: number of update iterations taken
-----------------------------------------
1 - 1: 1
_saupd: number of "converged" Ritz values
-----------------------------------------
1 - 1: 4
_saupd: final Ritz values
-------------------------
1 - 4: 5.040E+02 5.050E+02 5.177E+02 5.475E+02
_saupd: corresponding error bounds
----------------------------------
1 - 4: 0.000E+00 0.000E+00 0.000E+00 0.000E+00
==========================================
= Symmetric implicit Arnoldi update code =
= Version Number: 2.4 =
= Version Date: 07/31/96 =
==========================================
= Summary of timing statistics =
==========================================
Total number update iterations = 1
Total number of OP*x operations = 20
Total number of B*x operations = 0
Total number of reorthogonalization steps = 20
Total number of iterative refinement steps = 0
Total number of restart steps = 0
Total time in user OP*x operation = .000000
Total time in user B*x operation = .000000
Total time in Arnoldi update routine = .000000
Total time in saup2 routine = .000000
Total time in basic Arnoldi iteration loop = .000000
Total time in reorthogonalization phase = .000000
Total time in (re)start vector generation = .000000
Total time in trid eigenvalue subproblem = .000000
Total time in getting the shifts = .000000
Total time in applying the shifts = .000000
Total time in convergence testing = .000000
Ritz values and relative residuals
----------------------------------
Col 1 Col 2
Row 1: 8.63063E+02 0.00000E+00
Row 2: 8.86061E+02 0.00000E+00
Row 3: 9.19768E+02 0.00000E+00
Row 4: 9.48391E+02 0.00000E+00
_SSIMP
======
Size of the matrix is 100
The number of Ritz values requested is 4
The number of Arnoldi vectors generated (NCV) is 20
What portion of the spectrum: LM
The number of converged Ritz values is 4
The number of Implicit Arnoldi update iterations taken is 1
The number of OP*x is 20
The convergence criterion is 0.