0

I'm trying to write a driver in C++ to calculate the eigenvalues for an asymmetric, real-valued sparse matrix using the fortran functions offered by ARPACK, but I am having a bit of trouble with the reverse communication approach.

Generally, I am trying to solve the normal eigenvalue equation:

A*v = lambda*v

and any interaction with the matrix A is done in ARPACK via a function 'av':

av(n, workd[ipntr[0]], workd[ipntr[1]])

which multiplies the vector held in the array 'workd' beginning at location 'ipntr[0]' and inserts the result into the array 'workd' beginning at location 'ipntr[1]'. Examples of this approach are given in the manual at http://www.caam.rice.edu/software/ARPACK/ and also in the ARPACK/EXAMPLES/SIMPLE/dnsimp.f code.

What I would like to know is how do I actually involve the matrix A? If it is not passed to the routine then how is it possible to find its action on the vector provided?

In the example code dnsimp.f their matrix A is calculated within the function 'av', and is 'derived from the standard central difference discretisation of the 2 dimensional convection-diffusion operator'. However, I believe this is problem specific? It also doesn't seem too useful to have to code the derivation of the matrix A into the function. I can't find much information on this from the manual either.

It doesn't seem to be too much of a problem, since as it is a user defined function I am able to just change the definition of 'av' to include the matrix A as a parameter. However I would like to know how it is done properly in case of any potential compatibility issues.

Thank you!

user3023621
  • 103
  • 1
  • 9

1 Answers1

2

You don't have to supply the matrix to ARPACK.

All you have to do, is to multiply the matrix with the returned vectors (thus, reverse communication) till the desired convergence is reached.

For information on the algorithms, you should take a look at the users guide and especially on the chapter about the underlying algorithms.

Response to comment: The underlying algorithm is a form of Arnoldi Iteration. The basic algorithm is shown in wikipedia and shows, that the matrix A won't be accessed. Neither directly, nor indirectly.

In particular, the algorithms starts with an arbitrary normalized vector q_1. This vector is returned to the user. The user multiplies this vector with the matrix A using their favourite routine (usually some efficient sparse matrix-vector-multiplication) and returns the result to the Arnoldi Iteration to calculate a part of the Hessenberg matrix H (whose eigenvalues typically converge to the extreme eigenvalues of A) and the next vector q_2. This has to be iterated, till your results are converged.

Stefan
  • 2,460
  • 1
  • 17
  • 33
  • I understand that ARPACK doesn't need the matrix directly, I was wondering as to how the function 'av' which handles the reverse communication section actually accesses the matrix A to perform the matrix-vector product. Without passing it to the function (and so changing their suggested declaration format), the only way I can see it working is to declare A as a global variable. Changing the declarations is fine for my personal use, although I was curious to see if there was a better way that might then be compatible with other code that uses their method. – user3023621 Sep 29 '14 at 20:42
  • My apologies, I'm probably being very unclear about this! My problem is not the underlying algorithm, but is in actually getting the matrix A to my 'favourite routine'. In the examples, their favourite routine is the function 'av' which is defined above, and A has been coded into the function (not ideal). As you say, the matrix A is needed to calculate the product, but the question is how does one use the matrix A if it is not passed as a parameter to 'av'. I guess a global variable would work, but I will just stick with changing the declaration. Thank you for your help though! – user3023621 Sep 30 '14 at 11:46
  • Ok, I have to apologize, too... Now, I understand, what you mean. They defined the matrix-vector-product for their example matrix in the shortest possible way. That's the reason, why it doesn't appear explicitely. Normally, you would replace this routine which is *NOT* part of the ARPACK library (but only an example) with a routine, that multiplies *YOUR* matrix with the returned vectors. – Stefan Sep 30 '14 at 13:36
  • Fantastic, thanks for clearing that up. I wasn't sure if there was a standard approach in order to ensure compatibility, but it makes sense now that I am free to change the format. Cheers! – user3023621 Sep 30 '14 at 15:24