SF_SP_LINALG

Description

SciFortran module for Arpack/Lanczos diagonalization

Quick access

Routines:

sp_dvdson_eigh(), sp_eigh(), sp_lanc_eigh(), sp_lanc_tridiag()

Used modules

Subroutines and functions

interface  sf_sp_linalg/sp_eigh(matvec, eval, evec, mpicomm[, nblock, nitermax, bmat, v0, tol, iverbose, vrandom, iexit, numop, niter])

This subroutine computes the lowest Neigen=size(eval) eigenpairs of a large sparse real symmetric (lanczos_arpack_d) or complex hermitian (lanczos_arpack_c) matrix \(H\), which is defined only through its action on a vector, vout = H*vin, by the user routine MatVec. It uses the implicitly restarted Lanczos method of the ARPACK library, dsaupd + dseupd (real) and znaupd + zneupd (complex), selecting the smallest algebraic eigenvalues (smallest real part in the complex case).

The eigenvalues are returned in eval in ascending order, and the eigenvectors in the columns of evec, which has shape [Ns,Neigen], evec(:,j) being associated to eval(j). The eigenvalues are real also for complex hermitian matrices. The starting vector is v0 if present, otherwise a random vector (a constant one if vrandom is false), and it is normalized. The optional arguments iexit, NumOp and Niter return the ARPACK exit status, the number of matrix-vector products and the number of iterations. The program stops if bmat is neither 'I' nor 'G', or if the post-processing ARPACK routine returns an error code.

The parallel versions lanczos_parpack_d and lanczos_parpack_c use PARPACK (pdsaupd, pdseupd, pznaupd, pzneupd) and are available only when compiling with _MPI. They take the MPI communicator MpiComm as first argument, and do not have the argument bmat. The vectors are distributed over the processes in chunks: MatVec acts on the local chunk of size nchunk, v0 and evec hold the local chunk, and the eigenvalues are the same on all the processes. Nothing is done if MpiComm is MPI_COMM_NULL.

Parameters:

  • matvec [external] – vout=H*vin on vectors of size Nloc (local chunk if MPI)

  • eval (•) [real] – lowest eigenvalues in ascending order: [Neigen]

  • evec (•, •) [real/complex, in,required] – eigenvectors, evec(:,j) for eval(j): [Ns,Neigen], Ns local chunk if MPI

  • mpicomm [integer] – MPI communicatorInterface to Matrix-Vector routine:

Options:

  • nblock [integer] – number of Lanczos vectors, NCV (default 10*Neigen, at least Neigen+2)

  • nitermax [integer] – maximum number of Arnoldi iterations (default 512)

  • bmat [character(len=1)] – ‘I’ standard (default) or ‘G’ generalized eigenproblem

  • v0 (size(evec, 1)) [real/complex] – starting vector, local chunk if MPI (default: random)

  • tol [real] – relative accuracy of the Ritz values (default 0: machine precision)

  • iverbose [logical] – if T print ARPACK debug and info messages (default F)

  • vrandom [logical] – if T (default) random starting vector, else constant, if v0 is absent

  • iexit [integer] – ARPACK exit status (0: success)

  • numop [integer] – number of matrix-vector products

  • niter [integer] – number of iterationsDimensions:

interface  sf_sp_linalg/sp_lanc_eigh(matvec, egs, vect, nitermax, mpicomm[, iverbose, threshold, ncheck, vrandom, numop])

This subroutine computes the lowest eigenvalue Egs and the corresponding eigenvector Vect of a large sparse real symmetric (lanczos_eigh_d) or complex hermitian (lanczos_eigh_c) matrix \(H\), defined only through its action on a vector, vout = H*vin, by the user routine MatVec. It uses the plain Lanczos method, without re-orthogonalization.

The tridiagonal matrix is diagonalized with eigh at each iteration, and the iterations stop when the Lanczos coefficient b falls below threshold, or when the change of the lowest eigenvalue between two consecutive iterations falls below threshold, a test that starts after ncheck iterations, or at Nitermax iterations. The eigenvector is then rebuilt with a second Lanczos pass, so the number of matrix-vector products, returned in NumOp, is about twice the number of iterations. On input Vect is the starting vector: if it is zero a random vector (a constant one if vrandom is false) is generated. On output Vect is the normalized eigenvector and Nitermax is the number of iterations done. The eigenvalue is real also for complex hermitian matrices.

The values of iverbose, threshold, ncheck and vrandom are stored in module variables and keep their value in the following calls, also those that do not pass them. The default values are false, 1d-12, 10 and true.

The parallel versions mpi_lanczos_eigh_d and mpi_lanczos_eigh_c are available only when compiling with _MPI. They take the MPI communicator MpiComm as first argument. The vectors are distributed over the processes in chunks: MatVec acts on the local chunk of size Nloc and Vect holds the local chunk. The eigenvalue is the same on all the processes. Nothing is done if MpiComm is MPI_COMM_NULL.

Parameters:

  • matvec [external] – vout=H*vin on vectors of size Nloc (local chunk if MPI)

  • egs [real] – lowest eigenvalue

  • vect (•) [real/complex] – starting vector, local chunk if MPI; on output the eigenvector

  • nitermax [integer] – maximum number of iterations; on output the iterations done

  • mpicomm [integer] – MPI communicator

Options:

  • iverbose [logical] – if T print the convergence information (default F)

  • threshold [real] – convergence threshold on the energy change and on b (default 1d-12)

  • ncheck [integer] – iterations done before the convergence test starts (default 10)

  • vrandom [logical] – if T (default) random start vector when vect is 0, else constant

  • numop [integer] – number of matrix-vector products

interface  sf_sp_linalg/sp_lanc_tridiag(matvec, vin, alanc, blanc, mpicomm[, threshold])

This subroutine tri-diagonalizes a large sparse real symmetric (lanczos_tridiag_d) or complex hermitian (lanczos_tridiag_c) matrix \(H\), defined only through its action on a vector, vout = H*vin, by the user routine MatVec. It applies the plain Lanczos iteration, without re-orthogonalization, starting from the non-zero vector vin.

The coefficients of the tridiagonal matrix are returned in alanc (diagonal) and blanc (off-diagonal), real also in the complex case, with blanc(i+1) connecting the elements i and i+1, so that blanc(1) is not set. The number of iterations is size(alanc), and the iteration stops earlier if the Lanczos coefficient b falls below threshold (default 1d-12, stored in a module variable that keeps its value in the following calls), the remaining elements being left unchanged. In the serial versions vin is overwritten: it is normalized and then advanced up to the last Lanczos vector. The program stops if the norm of vin is zero.

The parallel versions mpi_lanczos_tridiag_d and mpi_lanczos_tridiag_c are available only when compiling with _MPI. They take the MPI communicator MpiComm as first argument. The vectors are distributed over the processes in chunks: MatVec acts on the local chunk and vin holds the local chunk, and it is NOT modified. The coefficients are the same on all the processes. Nothing is done if MpiComm is MPI_COMM_NULL.

Parameters:

  • matvec [external] – vout=H*vin on vectors of size Nloc (local chunk if MPI)

  • vin (•) [real/complex, inout] – vout=H*vin on vectors of size Nloc (local chunk if MPI)

  • alanc (•) [real, inout/in,required] – diagonal of the tridiagonal matrix, size = number of iterations

  • blanc (size(alanc)) [real, inout] – off-diagonal: blanc(i+1) couples i and i+1, blanc(1) is not set

  • mpicomm [integer] – MPI communicator

Options:

threshold [real] – stop if the Lanczos coefficient b is below this value (default 1d-12)

interface  sf_sp_linalg/sp_dvdson_eigh(matvec, eval, evec[, nblock, nitermax, tol])

This subroutine computes the lowest Neigen=size(eval) eigenpairs of a large sparse real symmetric matrix \(H\), defined only through its action on a vector, vout = H*vin, by the user routine MatVec, using the Davidson method (routine DVDSON). The eigenvalues are returned in eval in ascending order, and the eigenvectors in the columns of evec, which has shape [Ns,Neigen]. The program stops if evec does not have the shape [Ns,Neigen], or if the Davidson routine returns an error code.

The first column of evec on input is taken as the diagonal of \(H\), used by the method as preconditioner. If it is identically zero, the diagonal is computed by applying MatVec to the Ns unit vectors, which costs Ns matrix-vector products. The size of the expanding basis is min(Neigen+40,Ns), the convergence thresholds on the residuals and on the coefficients are fixed at 1d-15, and Tol sets the threshold on the eigenvalues. The method stops after Nitermax iterations. Only the serial version is available.

Parameters:

  • matvec [external] – vout=H*vin on vectors of size Nloc

  • eval (•) [real, inout] – lowest eigenvalues in ascending order: [Neigen]

  • evec (•, •) [real, inout] – in: 1st column = diagonal of H (computed if 0); out: eigenvectors [Ns,Neigen]

Options:

  • nblock [integer] – number of eigenpairs targeted in each iteration (default Neigen)

  • nitermax [integer] – maximum number of iterations (default max(200,40*Neigen))

  • tol [real] – convergence threshold on the eigenvalues (default 1d-15)Dimensions: