ARPACK driver translated from fortran to c++
Alec Jacobson
December 01, 2010
In an attempt to understand how to make a C++ interface to ARPACK to solve eigenvalue problems. I've translated the "real nonsymmetric generalized eigenvalue problem" example (EXAMPLES/NONSYM/dndrv3.f) from fortran into a cpp program.
Here's my translation, save it in a file called dndrv3.cpp:
/* dndrv3.f translated by Alec Jacobson
*/
#include <assert.h>
#include <cstdio>
#include "f2c.h"
/* %-----------------------------% */
/* | Define leading dimensions | */
/* | for all arrays. | */
/* | MAXN: Maximum dimension | */
/* | of the A allowed. | */
/* | MAXNEV: Maximum NEV allowed | */
/* | MAXNCV: Maximum NCV allowed | */
/* %-----------------------------% */
#define maxn 256
#define maxnev 10
#define maxncv 25
#define ldv maxn
// need pointer to ldv constant value
integer LDV = ldv;
/* %--------------% */
/* | Local Arrays | */
/* %--------------% */
integer iparam[11], ipntr[14];
logical select[maxncv];
// Fortran is column-major order!
doublereal ax[maxn],
mx[maxn],
d[3*maxncv],
resid[maxn],
v[maxncv*ldv],
workd[3*maxn],
workev[3*maxncv],
workl[3*maxncv*maxncv+6*maxncv],
md[maxn],
me[maxn-1];
/* %---------------% */
/* | Local Scalars | */
/* %---------------% */
char * bmat;
char * which;
integer ido,
n,
nev,
ncv,
lworkl,
info,
ierr,
j,
nconv,
maxitr,
ishfts,
mode;
doublereal tol,
sigmar,
sigmai,
h;
logical first,
rvec;
/* %------------% */
/* | Parameters | */
/* %------------% */
doublereal zero = 0.0,
one = 1.0;
/* %-----------------------------% */
/* | BLAS & LAPACK routines used | */
/* %-----------------------------% */
extern "C"{
extern doublereal dnrm2_(integer *, doublereal *, integer *);
extern doublereal dlapy2_(doublereal *, doublereal *);
extern /* Subroutine */ int dscal_(integer *, doublereal *, doublereal *,
integer *);
extern /* Subroutine */ int daxpy_(integer *, doublereal *, doublereal *,
integer *, doublereal *, integer *);
extern /* Subroutine */ int dpttrf_(integer *, doublereal *, doublereal *,
integer *);
extern /* Subroutine */ int dpttrs_(integer *, integer *, doublereal *,
doublereal *, doublereal *, integer *, integer *);
/* %-----------------------------% */
/* | other routines used | */
/* %-----------------------------% */
/* Subroutine */ int s_copy(char *, char *, ftnlen, ftnlen);
extern /* Subroutine */ int dmout_(integer *, integer *, integer *,
doublereal *, integer *, integer *, char *, ftnlen);
extern /* Subroutine */ int av_(integer *, doublereal *, doublereal *);
extern /* Subroutine */ int mv_(integer *, doublereal *, doublereal *);
extern /* Subroutine */ int dnaupd_(integer *, char *, integer *, char *,
integer *, doublereal *, doublereal *, integer *, doublereal *,
integer *, integer *, integer *, doublereal *, doublereal *,
integer *, integer *, ftnlen, ftnlen);
extern /* Subroutine */ int dneupd_(logical *, char *, logical *,
doublereal *, doublereal *, doublereal *, integer *, doublereal *,
doublereal *, doublereal *, char *, integer *, char *, integer *,
doublereal *, doublereal *, integer *, doublereal *, integer *,
integer *, integer *, doublereal *, doublereal *, integer *,
integer *, ftnlen, ftnlen, ftnlen);
}
/* matrix vector multiplication subroutine */
int av_(integer *n, doublereal *v, doublereal *w)
{
integer j;
doublereal one, two, dd, dl, du, s, h, rho;
rho = 10.;
one = 1.0;
two = 2.0;
/* Compute the matrix vector multiplication y<---A*x */
/* where A is stiffness matrix obtained from the finite element */
/* discretization of the 1-dimensional convection diffusion operator */
/* d^2u/dx^2 + rho*(du/dx) */
/* on the interval [0,1] with zero Dirichlet boundary condition using */
/* linear elements. */
h = one / (doublereal)(*n+1);
s = rho/two;
dd = two/h;
dl = -one/h -s;
du = -one/h + s;
w[0] = dd*v[0] + du * v[1];
for(int j = 1;j<*n-1;j++)
{
w[j] = dl*v[j-1] + dd*v[j] + du*v[j+1];
}
w[*n-1] = dl*v[*n-2] + dd*v[*n-1];
}
int mv_(integer *n, doublereal *v, doublereal *w)
{
integer j;
doublereal one, four, h;
one = 1.0;
four = 4.0;
/* Compute the matrix vector multiplication y<---M*x */
/* where M is the mass matrix formed by using piecewise linear */
/* elements on [0,1]. */
w[0] = four*v[0] + one*v[1];
for(int j = 1;j<*n-1;j++)
{
w[j] = one*v[j-1] + four*v[j] + one*v[j+1];
}
w[*n-1] = one*v[*n-2] + four*v[*n-1];
h = one / (doublereal)(*n + 1);
// increment of array, need as pointer by dscal_
integer incx = 1;
dscal_(n,&h,w,&incx);
}
int main(void){
/* %----------------------------------------------------% */
/* | The number N is the dimension of the matrix. A | */
/* | generalized eigenvalue problem is solved (BMAT = | */
/* | 'G'). NEV is the number of eigenvalues to be | */
/* | approximated. The user can modify NEV, NCV, WHICH | */
/* | to solve problems of different sizes, and to get | */
/* | different parts of the spectrum. However, The | */
/* | following conditions must be satisfied: | */
/* | N <= MAXN, | */
/* | NEV <= MAXNEV, | */
/* | NEV + 2 <= NCV <= MAXNCV | */
/* %----------------------------------------------------% */
// Set problem size
n = 100;
// Set number of eigen values desired
nev = 4;
// NCV is the largest number of basis vectors that will be used in
// the Implicitly Restarted Arnoldi Process. Work per major iteration
// is proportional to N*NCV*NCV.
ncv = 20;
assert(n<=maxn);
assert(nev<=maxnev);
assert(ncv<=maxncv);
bmat = "G";
which = "LM";
/* %------------------------------------------------% */
/* | M is the mass matrix formed by using piecewise | */
/* | linear elements on [0,1]. | */
/* %------------------------------------------------% */
h = one / (doublereal)(n+1);
for(int j = 0;j<n-1;j++)
{
md[j] = 4.0 * h;
me[j] = one*h;
}
md[n-1] = 4*h;
// Factor md using LAPACK's symmetric positive definite factor routine
dpttrf_(&n,md,me,&ierr);
if( ierr != 0 )
{
printf("ERROR with _pttrf\n");
exit(1);
}
/* %-----------------------------------------------------% */
/* | The work array WORKL is used in DNAUPD as | */
/* | workspace. Its dimension LWORKL is set as | */
/* | illustrated below. The parameter TOL determines | */
/* | the stopping criterion. If TOL<=0, machine | */
/* | precision is used. The variable IDO is used for | */
/* | reverse communication, and is initially set to 0. | */
/* | Setting INFO=0 indicates that a random vector is | */
/* | generated in DNAUPD to start the Arnoldi iteration. | */
/* %-----------------------------------------------------% */
lworkl = 3*ncv*ncv+6*ncv;
tol = 0.0;
ido = 0;
info = 0;
/* %---------------------------------------------------% */
/* | This program uses exact shifts with respect to | */
/* | the current Hessenberg matrix (IPARAM[0] = 1). | */
/* | IPARAM[2] specifies the maximum number of Arnoldi | */
/* | iterations allowed. Mode 2 of DNAUPD is used | */
/* | (IPARAM[6] = 2). All these options can be | */
/* | changed by the user. For details, see the | */
/* | documentation in DNAUPD. | */
/* %---------------------------------------------------% */
// use the exact shift strategy
// ISHIFT = 1: exact shift with respect to the current
// Hessenberg matrix H. This is equivalent to
// restarting the iteration from the beginning
// after updating the starting vector with a linear
// combination of Ritz vectors associated with the
// "wanted" eigenvalues.
ishfts = 1;
maxitr = 300;
mode = 2;
iparam[0] = ishfts;
iparam[2] = maxitr;
iparam[6] = mode;
/* %-------------------------------------------% */
/* | M A I N L O O P (Reverse communication) | */
/* %-------------------------------------------% */
while(true)
{
/* %---------------------------------------------% */
/* | Repeatedly call the routine DNAUPD and take | */
/* | actions indicated by parameter IDO until | */
/* | either convergence is indicated or maxitr | */
/* | has been exceeded. | */
/* %---------------------------------------------% */
dnaupd_(
&ido,
bmat,
&n,
which,
&nev,
&tol,
resid,
&ncv,
v,
&LDV,
iparam,
ipntr,
workd,
workl,
&lworkl,
&info,
(ftnlen)1,
(ftnlen)2);
if (ido == -1 || ido == 1)
{
/* %----------------------------------------% */
/* | Perform y <--- OP*x = inv[M]*A*x | */
/* | The user should supply his/her own | */
/* | matrix vector routine and a linear | */
/* | system solver. The matrix-vector | */
/* | subroutine should take workd(ipntr(1)) | */
/* | as input, and the final result should | */
/* | be returned to workd(ipntr(2)). | */
/* %----------------------------------------% */
// Perform matrix vector multiply: A*word[ipntr[0]-1] into
// workd[ipntr[1]-1]
av_(&n, &workd[ipntr[0] - 1], &workd[ipntr[1] - 1]);
// dpttrs_ needs pointer to NRHS, number of right hand sides
integer num_rhs= 1;
// DPTTRS solves a tridiagonal system of the form
// M * Y = A*X
// you can read this as computing y:
// Y = M^-1 * A * X
// workd[ipntr[1]-1] is both the input (A*X) and the output (Y)
// using the L*D*L' factorization of A computed by DPTTRF. D is a
// diagonal matrix specified in the vector D, L is a unit bidiagonal
// matrix whose subdiagonal is specified in the vector E, and X and B
// are N by NRHS matrices.
dpttrs_(&n, &num_rhs, md, me, &workd[ipntr[1] - 1], &n, &ierr);
if( ierr != 0 )
{
printf("ERROR with _pttrs\n");
exit(1);
}
}else if ( ido == 2 )
{
/* %-------------------------------------% */
/* | Perform y <--- M*x | */
/* | The matrix vector multiplication | */
/* | routine should take workd(ipntr(1)) | */
/* | as input and return the result to | */
/* | workd(ipntr(2)). | */
/* %-------------------------------------% */
mv_(&n, &workd[ipntr[0] - 1], &workd[ipntr[1] - 1]);
}
/* %-----------------------------------------% */
/* | Either we have convergence, or there is | */
/* | an error. | */
/* %-----------------------------------------% */
else if (info < 0)
{
/* %--------------------------% */
/* | Error message. Check the | */
/* | documentation in DNAUPD. | */
/* %--------------------------% */
printf("Error with _naupd, info = %d\n"
"Check the documentation of _naupd.\n",
info);
exit(1);
}else{
/* %-------------------------------------------% */
/* | No fatal errors occurred. | */
/* %-------------------------------------------% */
// Break out of reverse communication loop
break;
}
}
/* %-------------------------------------------% */
/* | Post-Process using DNEUPD. | */
/* | | */
/* | Computed eigenvalues may be extracted. | */
/* | | */
/* | Eigenvectors may also be computed now if | */
/* | desired. (indicated by rvec = .true.) | */
/* %-------------------------------------------% */
rvec = true;
dneupd_(
&rvec,
"A",
select,
d,
&d[0+1*maxncv],
v,
&LDV,
&sigmar,
&sigmai,
workev,
bmat,
&n,
which,
&nev,
&tol,
resid,
&ncv,
v,
&LDV,
iparam,
ipntr,
workd,
workl,
&lworkl,
&ierr,
(ftnlen) 1,
(ftnlen)1,
(ftnlen)2);
/* %-----------------------------------------------% */
/* | The real part of the eigenvalue is returned | */
/* | in the first column of the two dimensional | */
/* | array D, and the IMAGINARY part is returned | */
/* | in the second column of D. The corresponding | */
/* | eigenvectors are returned in the first NEV | */
/* | columns of the two dimensional array V if | */
/* | requested. Otherwise, an orthogonal basis | */
/* | for the invariant subspace corresponding to | */
/* | the eigenvalues in D is returned in V. | */
/* %-----------------------------------------------% */
if( ierr != 0 )
{
/* %------------------------------------% */
/* | Error condition: | */
/* | Check the documentation of DNEUPD. | */
/* %------------------------------------% */
printf("Error with _neupd, info = %d\n"
"Check the doumentation of _neupd\n",
ierr);
}else
{
first = true;
nconv = iparam[4];
for(int j = 0; j<iparam[4];j++)
{
/* %---------------------------% */
/* | Compute the residual norm | */
/* | | */
/* | || A*x - lambda*M*x || | */
/* | | */
/* | for the NCONV accurately | */
/* | computed eigenvalues and | */
/* | eigenvectors. (iparam(5) | */
/* | indicates how many are | */
/* | accurate to the requested | */
/* | tolerance) | */
/* %---------------------------% */
// need point to x,y increment values for daxpy
integer incx = 1;
integer incy = 1;
if( d[1*maxncv+j] == zero)
{
/* %--------------------% */
/* | Ritz value is real | */
/* %--------------------% */
// Multiply A * v(:,j) -> ax
av_(&n, &v[j*ldv + 0], ax);
// Multiply M * v(:,j) -> mx
mv_(&n, &v[j*ldv + 0], mx);
// need pointer to minus d[0,j]
doublereal minus_d_0_j = -d[0*maxncv+j];
// ax <-- -d[0,j] * mx + ax
daxpy_(&n, &minus_d_0_j, mx, &incx, ax, &incy);
// d[2,j] <-- ||ax||
d[2*maxncv+j] = dnrm2_(&n,ax,&incx);
d[2*maxncv+j] /= dabs(d[0*maxncv + j]);
}else if (first)
{
/* %------------------------% */
/* | Ritz value is complex | */
/* | Residual of one Ritz | */
/* | value of the conjugate | */
/* | pair is computed. | */
/* %------------------------% */
// Multiply A * v(:,j) -> ax
av_(&n, &v[j*ldv + 0], ax);
// Multiply M * v(:,j) -> mx
mv_(&n, &v[j*ldv + 0], mx);
// need pointer to minus d[0,j]
doublereal minus_d_0_j = -d[0*maxncv+j];
// ax <-- -d[0,j] * mx + ax
daxpy_(&n, &minus_d_0_j, mx, &incx, ax, &incy);
// Multiply M * v(:,j+1) -> mx
mv_(&n, &v[(j+1)*ldv + 0], mx);
// ax <-- d[1,j] * mx + ax
daxpy_(&n, &d[1*maxncv +j], mx, &incx, ax, &incy);
// d[2,j] <-- ||ax||
d[2*maxncv+j] = dnrm2_(&n,ax,&incx);
// d[2,j] <-- d[2,j]^2
d[2*maxncv+j] *= d[2*maxncv+j];
// Multiply A * v(:,j+1) -> ax
av_(&n, &v[(j+1)*ldv + 0], ax);
// Multiply M * v(:,j+1) -> mx
mv_(&n, &v[(j+1)*ldv + 0], mx);
// need pointer to minus d[0,j]
minus_d_0_j = -d[0*maxncv+j];
// ax <-- -d[0,j] * mx + ax
daxpy_(&n, &minus_d_0_j, mx, &incx, ax, &incy);
// Multiply M * v(:,j) -> mx
mv_(&n, &v[(j)*ldv + 0], mx);
// need pointer to minus d[1,j]
doublereal minus_d_1_j = -d[1*maxncv+j];
// ax <-- d[1,j] * mx + ax
daxpy_(&n, &minus_d_1_j, mx, &incx, ax, &incy);
// d[2,j] <-- sqrt(d[2,j]^2 + ||ax||^2)
doublereal ax_sqr = dnrm2_(&n,ax,&incx);
d[2*maxncv+j] = dlapy2_(&d[2*maxncv+j],&ax_sqr);
// d[2,j] <-- d[2,j] / sqrt(d[0,j]^2 + d[1,j]^2)
d[2*maxncv+j] = d[2*maxncv+j] / dlapy2_(&d[0*maxncv+j],&d[1*maxncv+j]);
d[2*maxncv+j+1] = d[2*maxncv+j];
first = false;
}else
{
first = true;
}
}
/* %-----------------------------% */
/* | Display computed residuals. | */
/* %-----------------------------% */
// need pointer to maxncv
integer MAXNCV = maxncv;
integer m = 6, lda = 3, idigit = -6;
dmout_(&m, &nconv, &lda, d, &MAXNCV, &idigit,
"Ritz values (Real,Imag) and relative residuals",
(ftnlen)46);
}
/* %------------------------------------------% */
/* | Print additional convergence information | */
/* %------------------------------------------% */
if( info == 1 )
{
printf("Maximum number of iterations reached.\n");
}else if(info == 3)
{
printf("No shifts could be applied during implicit"
" Arnoldi update, try increasing NCV.\n");
}
printf("_NRDV3\n"
" ====== \n"
" \n"
" Size of the matrix is %d\n"
" The number of Ritz values requested is %d\n"
" The number of Arnoldi vectors generated"
" (NCV) is %d\n"
" What portion of the spectrum: %s\n"
" The number of converged Ritz values is %d\n"
" The number of Implicit Arnoldi update,"
" iterations taken is %d\n"
" The number of OP*x is %d\n"
" The convergence criterion is %g\n"
" \n",
n,
nev,
ncv,
which,
nconv,
iparam[2],
iparam[8],
tol);
/* %---------------------------% */
/* | Done with program dndrv3. | */
/* %---------------------------% */
return 0;
}
This assumes that you used f2c to compile ARPACK, so you should be able to compile this on a mac with something like:
g++ dndrv3.cpp -o dndrv3_cpp -framework Accelerate -L[path to libarpack.a] -larpack -lf2c
Then if you issue:
./dndrv3_cpp
You should see:
Ritz values (Real,Imag) and relative residuals
----------------------------------------------
Col 1 Col 2 Col 3
Row 1: 4.96055E+01 0.00000E+00 9.79801E-01
Row 2: 4.10548E+02 0.00000E+00 2.73244E-01
Row 3: 2.00256E+03 0.00000E+00 1.38874E-01
Row 4: 3.01419E+03 0.00000E+00 7.42380E-02
_NRDV3
======
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 14
The number of OP*x is 218
The convergence criterion is 1.11022e-16
Which is the same as what you should see if you ran the original fortran version (up to double precision print accuracy).