Library

eigs(Symmetric(A), k; kwargs...)  # maps to symeigs 
eigs(Hermitian(A), k; kwargs...)  # maps to hermeigs 
symeigs(A, k; kwargs...)          # A must be symmetric, although we don't check. 
symeigs(A, B, k; kwargs...)       # A must be sym, B must be sym. pos. def, not checked 
symeigs(op, k; kwargs...)         # uses ArpackOp to handle multiple problem types 
hermeigs(A, k; kwargs...)         # forces complex valued vectors; A not checked.     
hermeigs(A, B k; kwargs...)       # forces complex valued vectors; A not checked.     
symeigs(TV, TF, op, k; kwargs...) # the most general implmementation

These functions compute Symmetric and Hermitian Eigenspaces hence, symeigs for symmetric eigenspace.

Calling sequence: the call is implemented with everything mapped to the single, most general call with symeigs.

Type parameters

• TV: The element type of the vector information; we need TV <: Complex for Hermitian problems
• TF: The element type of the factorization information; we need TV <: Real for all problems.

Note that it is okay to mix different precisions, although the lowest precision will determine the default tolerance; e.g. TV = Float32, TF = Float64 will use a default tolerance from Float32.

Arguments

• Symmetric(A) or Hermitian(A) a matrix whose type indicates it is symmetric or Hermitian
• k::Integer the size of the eigenspace to compute. (i.e. the number of eigenvectors)
• A a square matrix that should be symmetric (but we won't check)
• B or Symmetric(B) or Hermitian(B) like A, but for a generalized eigenvalue problem. The matrix B will be factorized via factorize(B), if you need more control, use ArpackSymmetricGeneralizedOp, or see one of the more advanced generalized interfaces below.
• op::ArpackOp An Arpack instance that holds the minimal information we need to run Arpack

Optional arguments

• which: which part of the spectrum to compute. There are a few choices:
  • :LM (the default) find the largest magnitude eigenspace
  • :SM find the smallest magnitude eigenspace
  • :LA find the largest algebraic value eigenspace
  • :SA find the smallest algebraic value eigenspace
  • :BE find both smallest and largest algebraic eigenspace
• ncv: the number of vectors used in the Arnoldi basis. The default value is min(2k, n) where n is the dimension of the matrix. Increasing this value can help if there are invariant subspaces with many nearby eigenvalues.
• tol: the tolerance on the error bound on the estimated eigenspace. The default value is max(eps(TF), eps(TV))/2; if TV is a complex valued type, this means the underlying element type.
• v0: a starting vector to build the initial basis. This should be of the same dimension as the matrix itself. The default choice is to use a random vector by setting v0=nothing. Note that GenericArpack.jl includes it's own random number generator to match the Fortran Arpack library. If you wish to match a sequence of calls, please save the ArpackState information between subsequent calls.
• maxiter: the maximum number of restarted Arnoldi bases to build. Default value is 300. this is highly conservative and may be changed in future
• failonmaxiter: (default true) if the maximum number of iterations is hit and this is true, then we force a hard error. Otherwise, we return information that can be used to run additional iterations.
• ritzvec: This determines if we compute eigenvector information or not. If not, then the matrix of eigenvectors will have zero columns.
• state: the ArpackState structure that stores information between calls. By default, we create a new state for each call. This stores information on the random number generator and information for the reverse communication interface.

Output and results optional parameters

• debug: This allow debugging output. This must be an instance of ArpackDebug(). See the field information as well as the Arpack debugging information to understand what is logged.
• stats: This tracks timing statistics and the number of various operations performed; it must be an instance of ArpackStats().

Advanced optional parameters

• mode: by default, the mode is set by the input or type of arpack op. While you can override this by setting it, this may result in incorrect behavior. In particular, we will throw a warning if the op seems inappropriate for the mode. The best way to use this is to use one of the more advanced ops, such as ArpackShiftInvertOp
• bmat: if bmat=Val(:I), then this is a standard eigenvalue problem. If bmat = Val(:G) then it is a generalized eigenvalue problem. By default, the choice of call and ArpackOp determines bmat and this should be changed rarely.

Sample calls

julia> using LinearAlgebra, SparseArrays, GenericArpack
julia> # standard real symmetric eigenvalue problem
julia> n=100; symeigs(Tridiagonal(-ones(n-1),2*ones(n),-ones(n-1)).*(n+1), 5)
julia> # real symmetric generalized eigenvalue problem
julia> n=100; symeigs(Tridiagonal(-ones(n-1),2*ones(n),-ones(n-1)).*(n+1),
                      SymTridiagonal(4*ones(n),ones(n-1))*(1/(6*(n+1))), 5)
julia> # use Float32 type for information on smallest eigenvalue problem
julia> n=100; symeigs(Float32, ArpackSimpleOp(Diagonal(1.0:n)), 5; which=:SA)
julia> # complex Hermitian eigenvalue problem 
julia> n=100; symeigs(Tridiagonal(im*ones(n-1),2*ones(n).+0im,im*ones(n-1)).*(n+1), 5)
julia> # mixed precision eigenvalue problem 
julia> n=100; symeigs(Float16, Float64, ArpackSimpleOp(Diagonal(1.0:n)), 5; which=:SA)

See also

If you wish more control over the interface and what is provided, your best solution is to work with the ArpackOp types, or write your own (it is very easy!)

• ArpackSimpleFunctionOp: This wraps a general input to output function.
• ArpackNormalOp: This is the op type used in the svds function
• ArpackShiftInvertOp: This is your best way to use the Shift Inver
• ArpackBucklingOp:
• svds

svds(A, k; kwargs...)
svds(T::Type, A, k; kwargs...)
svds(TV::Type, TF::Type, A, k; kwargs...)
svds(TV::Type, TF::Type, op, k; kwargs...)
complexsvds(op, k; kwargs...)

TODO: Write documentation See examples...

ArpackSimpleOp(A)

This corresponds to a simple eigenvalue problem Ax = lambda x, and builds a Julia object that represents the minimal information Arpack needs about the matrix to run an eigenvalue problem.

Arguments

• A: Anything that implements Base.size, LinearAlgebra.mul!

Examples

julia> op = ArpackSimpleOp(A)
juila> size(op) == size(A,1)

See also ArpackSymmetricGeneralizedOp, ArpackSimpleFunctionOp

ArpackSimpleFunctionOp(F::Function, n::Integer)

This corresponds to a simple eigenvalue problem Ax = lambda x, but takes a functional operator that we apply.

Arguments

• F::Function this is a function (y,x) -> mul!(y,A,x), i.e. a function that writes $A*x$ into $y$
• n::Integer the dimension of the problem

Examples

julia> using GenericArpack, SparseArrays
julia> function myf(y,x) 
       fill!(y, 0)
       y[1] += x[2] + x[100]
       for i in 2:99
         y[i] += x[i+1] + x[i-1]
       end
       y[100] += x[99] + x[1]
       end
julia> op = ArpackSimpleFunctionOp(myf, 100)

ArpackSymmetricGeneralizedOp(A,invB,B)

We need three operations: Ax, invBx, B*x B must also be symmetric, pos. def. Note that if B can be factorized efficiently via Cholesky, then there is a better way to proceed.

Examples

```julia-repl julia> using GenericArpack, LinearAlgebra julia> n = 100 julia> A = Tridiagonal(-ones(n-1),2ones(n),-ones(n-1)).(n+1) julia> B = Tridiagonal(ones(n-1),4ones(n),ones(n-1)).(1/(6*(n+1))) julia> op = ArpackSymmetricGeneralizedOp(A, lu!(copy(B)), B)

See also ArpackSimpleOp

ArpackNormalOp(A)

This provides an ArpackSVDOp that is used as a bridge in svds to call the symeigs routine. It maniuplates the so-called normal matrix $AA^H$ or $A^H A$ depending on which is smaller. Note that this function allocates memory to compute the operation, so it will not be safe to use with multiple threads.

Arguments

• A: Anything that implements Base.size, LinearAlgebra.mul!(y,A,x), LinearAlgebra.adjoint and LinearAlgebra.mul!(y,adjoint(A),x)`

Examples

julia> op = ArpackNormalOp(A)

See also ArpackNormalFunctionOp

ArpackNormalFunctionOp(ax::Function, atx::Function, m::Integer, n::Integer)

Arguments

• ax A function to compute (y,x) and write y = A*x into the memory for y
• atx A function to compute (y,x) and write y = A^H*x into the memory for y
• m the number of rows of A
• n the number of rows of A

See also ArpackNormalOp, ArpackSimpleFunctionOp

ArpackShiftInvertOp

ArpackBucklingOp

c  Mode 4:  K*x = lambda*KG*x, K symmetric positive semi-definite,
c           KG symmetric indefinite
c           ===> OP = (inv[K - sigma*KG])*K  and  B = K.
c           ===> Buckling mode

ArpackOp

The general abstract ArpackOp interface.

See also ArpackSimpleOp, ArpackSymmetricGeneralizedOp, ArpackSimpleFunctionOp

ArpackSVDSOp

The general abstract ArpackSVDOp interface that bridges between eigenvaules and singular values.

svd_residuals(A, U, s, V)
svd_residuals(A, SVD)
svd_residuals(A, U, s, V, k) # compute for only the top k 
svd_residuals!(r, A, U, s, V) # write result in place

Compute the residuals of an SVD computation $||A*V[:,i] - U[:,i]*sigma[i]||$, and return the result in a vector.

Using a matvec function

Note that A can also be a function A(y,x), where y = A*x is updated in place. e.g. svd(A,U,s,V) == svd((y,x)->mul!(y,A,x), U,s,V) are equivalent.

call dsaupd ( IDO, BMAT, N, WHICH, NEV, TOL, RESID, NCV, V, LDV, IPARAM, IPNTR, WORKD, WORKL, LWORKL, INFO )

Arguments

IDO Integer. (INPUT/OUTPUT) Reverse communication flag. IDO must be zero on the first call to dsaupd . IDO will be set internally to indicate the type of operation to be performed. Control is then given back to the calling routine which has the responsibility to carry out the requested operation and call dsaupd with the result. The operand is given in WORKD(IPNTR(1)), the result must be put in WORKD(IPNTR(2)). (If Mode = 2 see remark 5 below) ––––––––––––––––––––––––––––––- IDO = 0: first call to the reverse communication interface IDO = -1: compute Y = OP * X where IPNTR(1) is the pointer into WORKD for X, IPNTR(2) is the pointer into WORKD for Y. This is for the initialization phase to force the starting vector into the range of OP. IDO = 1: compute Y = OP * X where IPNTR(1) is the pointer into WORKD for X, IPNTR(2) is the pointer into WORKD for Y. In mode 3,4 and 5, the vector B * X is already available in WORKD(ipntr(3)). It does not need to be recomputed in forming OP * X. IDO = 2: compute Y = B * X where IPNTR(1) is the pointer into WORKD for X, IPNTR(2) is the pointer into WORKD for Y. IDO = 3: compute the IPARAM(8) shifts where IPNTR(11) is the pointer into WORKL for placing the shifts. See remark 6 below. IDO = 99: done ––––––––––––––––––––––––––––––-

BMAT Character1. (INPUT) BMAT specifies the type of the matrix B that defines the semi-inner product for the operator OP. B = 'I' -> standard eigenvalue problem Ax = lambdax B = 'G' -> generalized eigenvalue problem Ax = lambdaBx

N Integer. (INPUT) Dimension of the eigenproblem.

WHICH Character*2. (INPUT) Specify which of the Ritz values of OP to compute.

     'LA' - compute the NEV largest (algebraic) eigenvalues.
     'SA' - compute the NEV smallest (algebraic) eigenvalues.
     'LM' - compute the NEV largest (in magnitude) eigenvalues.
     'SM' - compute the NEV smallest (in magnitude) eigenvalues.
     'BE' - compute NEV eigenvalues, half from each end of the
            spectrum.  When NEV is odd, compute one more from the
            high end than from the low end.
      (see remark 1 below)

NEV Integer. (INPUT) Number of eigenvalues of OP to be computed. 0 < NEV < N.

TOL Double precision scalar. (INPUT) Stopping criterion: the relative accuracy of the Ritz value is considered acceptable if BOUNDS(I) .LE. TOL*ABS(RITZ(I)). If TOL .LE. 0. is passed a default is set: DEFAULT = DLAMCH ('EPS') (machine precision as computed by the LAPACK auxiliary subroutine DLAMCH ).

RESID Double precision array of length N. (INPUT/OUTPUT) On INPUT: If INFO .EQ. 0, a random initial residual vector is used. If INFO .NE. 0, RESID contains the initial residual vector, possibly from a previous run. On OUTPUT: RESID contains the final residual vector.

NCV Integer. (INPUT) Number of columns of the matrix V (less than or equal to N). This will indicate how many Lanczos vectors are generated at each iteration. After the startup phase in which NEV Lanczos vectors are generated, the algorithm generates NCV-NEV Lanczos vectors at each subsequent update iteration. Most of the cost in generating each Lanczos vector is in the matrix-vector product OP*x. (See remark 4 below).

V Double precision N by NCV array. (OUTPUT) The NCV columns of V contain the Lanczos basis vectors.

LDV Integer. (INPUT) Leading dimension of V exactly as declared in the calling program.

IPARAM Integer array of length 11. (INPUT/OUTPUT) IPARAM(1) = ISHIFT: method for selecting the implicit shifts. The shifts selected at each iteration are used to restart the Arnoldi iteration in an implicit fashion. ––––––––––––––––––––––––––––––- ISHIFT = 0: the shifts are provided by the user via reverse communication. The NCV eigenvalues of the current tridiagonal matrix T are returned in the part of WORKL array corresponding to RITZ. See remark 6 below. ISHIFT = 1: exact shifts with respect to the reduced tridiagonal matrix T. This is equivalent to restarting the iteration with a starting vector that is a linear combination of Ritz vectors associated with the "wanted" Ritz values. ––––––––––––––––––––––––––––––-

     IPARAM(2) = LEVEC
     No longer referenced. See remark 2 below.

     IPARAM(3) = MXITER
     On INPUT:  maximum number of Arnoldi update iterations allowed.
     On OUTPUT: actual number of Arnoldi update iterations taken.

     IPARAM(4) = NB: blocksize to be used in the recurrence.
     The code currently works only for NB = 1.

     IPARAM(5) = NCONV: number of "converged" Ritz values.
     This represents the number of Ritz values that satisfy
     the convergence criterion.

     IPARAM(6) = IUPD
     No longer referenced. Implicit restarting is ALWAYS used.

     IPARAM(7) = MODE
     On INPUT determines what type of eigenproblem is being solved.
     Must be 1,2,3,4,5; See under Description of dsaupd  for the
     five modes available.

     IPARAM(8) = NP
     When ido = 3 and the user provides shifts through reverse
     communication (IPARAM(1)=0), dsaupd  returns NP, the number
     of shifts the user is to provide. 0 < NP <=NCV-NEV. See Remark
     6 below.

     IPARAM(9) = NUMOP, IPARAM(10) = NUMOPB, IPARAM(11) = NUMREO,
     OUTPUT: NUMOP  = total number of OP*x operations,
             NUMOPB = total number of B*x operations if BMAT='G",
             NUMREO = total number of steps of re-orthogonalization.

IPNTR Integer array of length 11. (OUTPUT) Pointer to mark the starting locations in the WORKD and WORKL arrays for matrices/vectors used by the Lanczos iteration. ––––––––––––––––––––––––––––––- IPNTR(1): pointer to the current operand vector X in WORKD. IPNTR(2): pointer to the current result vector Y in WORKD. IPNTR(3): pointer to the vector B * X in WORKD when used in the shift-and-invert mode. IPNTR(4): pointer to the next available location in WORKL that is untouched by the program. IPNTR(5): pointer to the NCV by 2 tridiagonal matrix T in WORKL. IPNTR(6): pointer to the NCV RITZ values array in WORKL. IPNTR(7): pointer to the Ritz estimates in array WORKL associated with the Ritz values located in RITZ in WORKL. IPNTR(11): pointer to the NP shifts in WORKL. See Remark 6 below.

     Note: IPNTR(8:10) is only referenced by dseupd . See Remark 2.
     IPNTR(8): pointer to the NCV RITZ values of the original system.
     IPNTR(9): pointer to the NCV corresponding error bounds.
     IPNTR(10): pointer to the NCV by NCV matrix of eigenvectors
                of the tridiagonal matrix T. Only referenced by
                dseupd  if RVE= .TRUE. See Remarks.
     -------------------------------------------------------------

WORKD Double precision work array of length 3N. (REVERSE COMMUNICATION) Distributed array to be used in the basiArnoldi iteration for reverse communication. The user should not use WORKD as temporary workspace during the iteration. Upon termination WORKD(1:N) contains BRESID(1:N). If the Ritz vectors are desired subroutine dseupd uses this output. See Data Distribution Note below.

WORKL Double precision work array of length LWORKL. (OUTPUT/WORKSPACE) Private (replicated) array on each PE or array allocated on the front end. See Data Distribution Note below.

LWORKL Integer. (INPUT) LWORKL must be at least NCV**2 + 8*NCV .

Changes:

Note that info is only used for Input. INFO Integer. (INPUT/OUTPUT) If INFO .EQ. 0, a randomly initial residual vector is used. If INFO .NE. 0, RESID contains the initial residual vector, possibly from a previous run.

Return value (this is info in Fortran)

The return value is a pair Error flag on output. = 0: Normal exit. = 1: Maximum number of iterations taken. All possible eigenvalues of OP has been found. IPARAM(5) returns the number of wanted converged Ritz values. = 2: No longer an informational error. Deprecated starting with release 2 of ARPACK. = 3: No shifts could be applied during a cycle of the Implicitly restarted Arnoldi iteration. One possibility is to increase the size of NCV relative to NEV. See remark 4 below. = -1: N must be positive. = -2: NEV must be positive. = -3: NCV must be greater than NEV and less than or equal to N. = -4: The maximum number of Arnoldi update iterations allowed must be greater than zero. = -5: WHICH must be one of 'LM", 'SM", 'LA", 'SA' or 'BE'. = -6: BMAT must be one of 'I' or 'G'. = -7: Length of private work array WORKL is not sufficient. = -8: Error return from trid. eigenvalue calculation; Informatinal error from LAPACK routine dsteqr . = -9: Starting vector is zero. = -10: IPARAM(7) must be 1,2,3,4,5. = -11: IPARAM(7) = 1 and BMAT = 'G' are incompatible. = -12: IPARAM(1) must be equal to 0 or 1. = -13: NEV and WHICH = 'BE' are incompatible. = -9999: Could not build an Arnoldi factorization. IPARAM(5) returns the size of the current Arnoldi factorization. The user is advised to check that enough workspace and array storage has been allocated.

Usage: call dseupd ( RVEC, HOWMNY, SELECT, D, Z, LDZ, SIGMA, BMAT, N, WHICH, NEV, TOL, RESID, NCV, V, LDV, IPARAM, IPNTR, WORKD, WORKL, LWORKL, INFO )

RVEC LOGICAL (INPUT) Specifies whether Ritz vectors corresponding to the Ritz value approximations to the eigenproblem Az = lambdaB*z are computed.

         RVEC = .FALSE.     Compute Ritz values only.

         RVEC = .TRUE.      Compute Ritz vectors.

HOWMNY Character*1 (INPUT) Specifies how many Ritz vectors are wanted and the form of Z the matrix of Ritz vectors. See remark 1 below. = 'A': compute NEV Ritz vectors; = 'S': compute some of the Ritz vectors, specified by the logical array SELECT.

SELECT Logical array of dimension NCV. (INPUT/WORKSPACE) If HOWMNY = 'S', SELECT specifies the Ritz vectors to be computed. To select the Ritz vector corresponding to a Ritz value D(j), SELECT(j) must be set to .TRUE.. If HOWMNY = 'A' , SELECT is used as a workspace for reordering the Ritz values.

D Double precision array of dimension NEV. (OUTPUT) On exit, D contains the Ritz value approximations to the eigenvalues of Az = lambdaBz. The values are returned in ascending order. If IPARAM(7) = 3,4,5 then D represents the Ritz values of OP computed by dsaupd transformed to those of the original eigensystem Az = lambdaBz. If IPARAM(7) = 1,2 then the Ritz values of OP are the same as the those of Az = lambdaB*z.

Z Double precision N by NEV array if HOWMNY = 'A'. (OUTPUT) On exit, Z contains the B-orthonormal Ritz vectors of the eigensystem Az = lambdaB*z corresponding to the Ritz value approximations. If RVEC = .FALSE. then Z is not referenced. NOTE: The array Z may be set equal to first NEV columns of the Arnoldi/Lanczos basis array V computed by DSAUPD .

LDZ Integer. (INPUT) The leading dimension of the array Z. If Ritz vectors are desired, then LDZ .ge. max( 1, N ). In any case, LDZ .ge. 1.

SIGMA Double precision (INPUT) If IPARAM(7) = 3,4,5 represents the shift. Not referenced if IPARAM(7) = 1 or 2.

**** The remaining arguments MUST be the same as for the **** **** call to DSAUPD that was just completed. ****

NOTE: The remaining arguments

       BMAT, N, WHICH, NEV, TOL, RESID, NCV, V, LDV, IPARAM, IPNTR,
       WORKD, WORKL, LWORKL, INFO

     must be passed directly to DSEUPD  following the last call
     to DSAUPD .  These arguments MUST NOT BE MODIFIED between
     the the last call to DSAUPD  and the call to DSEUPD .

Two of these parameters (WORKL, INFO) are also output parameters:

WORKL Double precision work array of length LWORKL. (OUTPUT/WORKSPACE) WORKL(1:4ncv) contains information obtained in dsaupd . They are not changed by dseupd . WORKL(4ncv+1:ncvncv+8ncv) holds the untransformed Ritz values, the computed error estimates, and the associated eigenvector matrix of H.

      Note: IPNTR(8:10) contains the pointer into WORKL for addresses
      of the above information computed by dseupd .
      -------------------------------------------------------------
      IPNTR(8): pointer to the NCV RITZ values of the original system.
      IPNTR(9): pointer to the NCV corresponding error bounds.
      IPNTR(10): pointer to the NCV by NCV matrix of eigenvectors
                 of the tridiagonal matrix T. Only referenced by
                 dseupd  if RVEC = .TRUE. See Remarks.
      -------------------------------------------------------------

INFO Integer. (OUTPUT) Error flag on output. = 0: Normal exit. = -1: N must be positive. = -2: NEV must be positive. = -3: NCV must be greater than NEV and less than or equal to N. = -5: WHICH must be one of 'LM', 'SM', 'LA', 'SA' or 'BE'. = -6: BMAT must be one of 'I' or 'G'. = -7: Length of private work WORKL array is not sufficient. = -8: Error return from trid. eigenvalue calculation; Information error from LAPACK routine dsteqr . = -9: Starting vector is zero. = -10: IPARAM(7) must be 1,2,3,4,5. = -11: IPARAM(7) = 1 and BMAT = 'G' are incompatible. = -12: NEV and WHICH = 'BE' are incompatible. = -14: DSAUPD did not find any eigenvalues to sufficient accuracy. = -15: HOWMNY must be one of 'A' or 'S' if RVEC = .true. = -16: HOWMNY = 'S' not yet implemented = -17: DSEUPD got a different count of the number of converged Ritz values than DSAUPD got. This indicates the user probably made an error in passing data from DSAUPD to DSEUPD or that the data was modified before entering DSEUPD .

Notes on Julia interface

• we eliminate 'howmny' and 'select' as they aren't used in the Arpack code.
• Julia doesn't need ldz, so we remove that too...