program main !******************************************************************************* ! !! SSSIMP is a simple program to call ARPACK for a symmetric eigenproblem. ! ! Discussion: ! ! This example program is intended to illustrate the ! simplest case of using ARPACK in considerable detail. ! This code may be used to understand basic usage of ARPACK ! and as a template for creating an interface to ARPACK. ! ! This code shows how to use ARPACK to find a few eigenvalues ! LAMBDA and corresponding eigenvectors X for the standard ! eigenvalue problem: ! ! A * X = LAMBDA * X ! ! where A is an N by N real symmetric matrix. ! ! The main points illustrated here are: ! ! 1) How to declare sufficient memory to find NEV ! eigenvalues of largest magnitude. Other options ! are available. ! ! 2) Illustration of the reverse communication interface ! needed to utilize the top level ARPACK routine SSAUPD ! that computes the quantities needed to construct ! the desired eigenvalues and eigenvectors (if requested). ! ! 3) How to extract the desired eigenvalues and eigenvectors ! using the ARPACK routine SSEUPD. ! ! The only things that must be supplied in order to use this ! routine on your problem is: ! ! * to change the array dimensions appropriately, ! * to specify WHICH eigenvalues you want to compute ! * to supply a matrix-vector product ! w <- A * v ! in place of the call to AV( ) below. ! ! Once usage of this routine is understood, you may wish to explore ! the other available options to improve convergence, to solve generalized ! problems, etc. Look at the file ex-sym.doc in DOCUMENTS directory. ! ! ! For this problem, we want to solve A*x = lambda*x in regular mode. ! ! A is derived from the central difference discretization ! of the 2-dimensional Laplacian on the unit square with ! zero Dirichlet boundary condition. ! ! OP = A and B = I. ! Assume "call av (n,x,y)" computes y = A*x ! Use mode 1 of SSAUPD. ! ! Author: ! ! Richard Lehoucq, Danny Sorensen, Chao Yang, ! Department of Computational and Applied Mathematics, ! Rice University, ! Houston, Texas. ! ! Storage: ! ! The maximum dimensions for all arrays are set here to accommodate ! a problem size of N <= MAXN ! ! NEV is the number of eigenvalues requested. ! See specifications for ARPACK usage below. ! ! 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. ! ! You must set: ! ! MAXN: Maximum dimension of the A allowed. ! MAXNEV: Maximum NEV allowed. ! MAXNCV: Maximum NCV allowed. ! implicit none integer, parameter :: maxn = 256 integer, parameter :: maxnev = 10 integer, parameter :: maxncv = 25 integer, parameter :: ldv = maxn intrinsic abs real ax(maxn) character bmat real d(maxncv,2) integer ido integer ierr integer info integer iparam(11) integer ipntr(11) integer ishfts integer j integer lworkl integer maxitr integer mode1 integer n integer nconv integer ncv integer nev integer nx real resid(maxn) logical rvec external saxpy logical select(maxncv) real sigma real, external :: snrm2 real tol real v(ldv,maxncv) character ( len = 2 ) which real workl(maxncv*(maxncv+8)) real workd(3*maxn) real, parameter :: zero = 0.0E+00 ! ! The following include statement and assignments control trace output ! from the internal actions of ARPACK. See debug.doc in the ! DOCUMENTS directory for usage. ! ! Initially, the most useful information will be a breakdown of ! time spent in the various stages of computation given by setting ! msaupd = 1. ! include 'debug.h' call timestamp ( ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SSSIMP:' write ( *, '(a)' ) ' FORTRAN90 version' write ( *, '(a)' ) ' A simple ARPACK calling program.' write ( *, '(a)' ) ' This program defines an eigenproblem for a' write ( *, '(a)' ) ' symmetric matrix.' ndigit = -3 logfil = 6 msgets = 0 msaitr = 0 msapps = 0 msaupd = 1 msaup2 = 0 mseigt = 0 mseupd = 0 ! ! Set dimensions for this problem. ! nx = 10 n = nx * nx ! ! Specifications for ARPACK usage are set below: ! ! 1) NEV = 4 asks for 4 eigenvalues to be computed. ! ! 2) NCV = 20 sets the length of the Arnoldi factorization. ! ! 3) This is a standard problem(indicated by bmat = 'I') ! ! 4) Ask for the NEV eigenvalues of largest magnitude ! (indicated by which = 'LM') ! ! See documentation in SSAUPD for the other options SM, LA, SA, LI, SI. ! ! NEV and NCV must satisfy the following conditions: ! ! NEV <= MAXNEV ! NEV + 1 <= NCV <= MAXNCV ! nev = 4 ncv = 20 bmat = 'I' which = 'LM' if ( maxn < n ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SSSIMP - Fatal error!' write ( *, '(a)' ) ' N is greater than MAXN ' stop else if ( maxnev < nev ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SSSIMP - Fatal error!' write ( *, '(a)' ) ' NEV is greater than MAXNEV ' stop else if ( maxncv < ncv ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SSSIMP - Fatal error!' write ( *, '(a)' ) ' NCV is greater than MAXNCV ' stop end if ! ! Specification of stopping rules and initial ! conditions before calling SSAUPD ! ! TOL determines the stopping criterion. Expect ! abs(lambdaC - lambdaT) < TOL*abs(lambdaC) ! computed true ! If TOL <= 0, then TOL <- macheps (machine precision) is used. ! ! IDO is the REVERSE COMMUNICATION parameter ! used to specify actions to be taken on return ! from SSAUPD. (See usage below.) ! It MUST initially be set to 0 before the first ! call to SSAUPD. ! ! INFO on entry specifies starting vector information ! and on return indicates error codes ! Initially, setting INFO=0 indicates that a ! random starting vector is requested to ! start the ARNOLDI iteration. Setting INFO to ! a nonzero value on the initial call is used ! if you want to specify your own starting ! vector. (This vector must be placed in RESID.) ! ! The work array WORKL is used in SSAUPD as workspace. Its dimension ! LWORKL is set as illustrated below. ! lworkl = ncv * ( ncv + 8 ) tol = zero info = 0 ido = 0 ! ! Specification of Algorithm Mode: ! ! This program uses the exact shift strategy ! (indicated by setting PARAM(1) = 1). ! ! IPARAM(3) specifies the maximum number of Arnoldi iterations allowed. ! ! Mode 1 of SSAUPD is used (IPARAM(7) = 1). ! ! All these options can be changed by the user. For details see the ! documentation in SSAUPD. ! ishfts = 1 maxitr = 300 mode1 = 1 iparam(1) = ishfts iparam(3) = maxitr iparam(7) = mode1 ! ! MAIN LOOP (Reverse communication loop) ! ! Repeatedly call SSAUPD and take actions indicated by parameter ! IDO until convergence is indicated or MAXITR is exceeded. ! do call ssaupd ( ido, bmat, n, which, nev, tol, resid, & ncv, v, ldv, iparam, ipntr, workd, workl, & lworkl, info ) if ( ido /= -1 .and. ido /= 1 ) then exit end if ! ! Perform matrix vector multiplication ! ! y <--- OP*x ! ! The user supplies a matrix-vector multiplication routine that takes ! workd(ipntr(1)) as the input, and return the result to workd(ipntr(2)). ! call av ( nx, workd(ipntr(1)), workd(ipntr(2)) ) end do ! ! Either we have convergence or there is an error. ! if ( info < 0 ) then ! ! Error message. Check the documentation in SSAUPD. ! write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SSSIMP - Fatal error!' write ( *, '(a,i6)' ) ' Error with SSAUPD, INFO = ', info write ( *, '(a)' ) ' Check documentation in SSAUPD.' else ! ! No fatal errors occurred. ! Post-Process using SSEUPD. ! ! Computed eigenvalues may be extracted. ! ! Eigenvectors may be also computed now if ! desired. (indicated by rvec = .true.) ! ! The routine SSEUPD now called to do this ! post processing (Other modes may require ! more complicated post processing than mode1.) ! rvec = .true. call sseupd ( rvec, 'All', select, d, v, ldv, sigma, & bmat, n, which, nev, tol, resid, ncv, v, ldv, & iparam, ipntr, workd, workl, lworkl, ierr ) ! ! Eigenvalues are returned in the first column of the two dimensional ! array D and the corresponding eigenvectors are returned in the first ! NCONV (=IPARAM(5)) 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 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SSSIMP - Fatal error!' write ( *, '(a,i6)' ) ' Error with SSEUPD, IERR = ', ierr write ( *, '(a)' ) ' Check the documentation of SSEUPD.' ! ! Compute the residual norm ! ! || A*x - lambda*x || ! ! for the NCONV accurately computed eigenvalues and ! eigenvectors. (iparam(5) indicates how many are ! accurate to the requested tolerance) ! else nconv = iparam(5) do j = 1, nconv call av ( nx, v(1,j), ax ) call saxpy ( n, -d(j,1), v(1,j), 1, ax, 1 ) d(j,2) = snrm2 ( n, ax, 1) d(j,2) = d(j,2) / abs ( d(j,1) ) end do ! ! Display computed residuals. ! call smout ( 6, nconv, 2, d, maxncv, -6, & 'Ritz values and relative residuals' ) end if ! ! Print additional convergence information. ! if ( info == 1) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Maximum number of iterations reached.' else if ( info == 3) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' No shifts could be applied during implicit' & // ' Arnoldi update, try increasing NCV.' end if write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SSSIMP:' write ( *, '(a)' ) '====== ' write ( *, '(a)' ) ' ' write ( *, '(a,i6)' ) ' Size of the matrix is ', n write ( *, '(a,i6)' ) ' The number of Ritz values requested is ', nev write ( *, '(a,i6)' ) & ' The number of Arnoldi vectors generated (NCV) is ', ncv write ( *, '(a)' ) ' What portion of the spectrum: ' // which write ( *, '(a,i6)' ) & ' The number of converged Ritz values is ', nconv write ( *, '(a,i6)' ) & ' The number of Implicit Arnoldi update iterations taken is ', iparam(3) write ( *, '(a,i6)' ) ' The number of OP*x is ', iparam(9) write ( *, '(a,g14.6)' ) ' The convergence criterion is ', tol end if write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SSSIMP:' write ( *, '(a)' ) ' Normal end of execution.' write ( *, '(a)' ) ' ' call timestamp ( ) stop end subroutine av ( nx, v, w ) !******************************************************************************* ! !! AV computes w <- A * V where A is a discretized Laplacian. ! ! Discussion: ! ! The matrix used is the 2 dimensional discrete Laplacian on unit ! square with zero Dirichlet boundary condition. ! ! Compute w <--- OP * v, where OP is the nx*nx by nx*nx block ! tridiagonal matrix ! ! | T -I | ! |-I T -I | ! | -I T | ! | ... -I | ! | -I T | ! ! Parameters: ! ! Input, integer NX, the length of the vectors. ! ! Input, real V(NX), the vector to be operated on by A. ! ! Output, real W(NX), the result of A*V. ! implicit none integer nx real h2 integer j integer lo integer n2 real, parameter :: one = 1.0E+00 real v(nx*nx) real w(nx*nx) call tv ( nx, v(1), w(1) ) call saxpy ( nx, -one, v(nx+1), 1, w(1), 1 ) do j = 2, nx-1 lo = (j-1) * nx call tv ( nx, v(lo+1), w(lo+1) ) call saxpy ( nx, -one, v(lo-nx+1), 1, w(lo+1), 1 ) call saxpy ( nx, -one, v(lo+nx+1), 1, w(lo+1), 1 ) end do lo = (nx-1) * nx call tv ( nx, v(lo+1), w(lo+1) ) call saxpy ( nx, -one, v(lo-nx+1), 1, w(lo+1), 1 ) ! ! Scale the vector W by (1/H^2), where H is the mesh size. ! n2 = nx * nx h2 = one / real ( ( nx + 1 ) * ( nx + 1 ) ) call sscal ( n2, one/h2, w, 1 ) return end subroutine tv ( nx, x, y ) !******************************************************************************* ! !! TV computes y <-- T*x. ! ! Discussion: ! ! The full matrix A is the 2 dimensional discrete Laplacian on unit ! square with zero Dirichlet boundary condition. ! ! A can be represented as the nx*nx by nx*nx block tridiagonal matrix ! ! | T -I | ! |-I T -I | ! | -I T | ! | ... -I| ! | -I T| ! ! It is the job of this routine to compute products of the form T*x. ! ! Parameters: ! ! Input, integer NX, the length of the vectors. ! ! Input, real X(NX), the vector to be multiplied by T. ! ! Output, real Y(NX), the product T*X. ! implicit none integer nx integer j real dd real dl real du real x(nx) real y(nx) ! ! Compute the matrix vector multiplication y<---T*x ! where T is a nx by nx tridiagonal matrix with DD on the ! diagonal, DL on the subdiagonal, and DU on the superdiagonal. ! dd = 4.0E+00 dl = -1.0E+00 du = -1.0E+00 y(1) = dd * x(1) + du * x(2) do j = 2, nx-1 y(j) = dl * x(j-1) + dd * x(j) + du * x(j+1) end do y(nx) = dl * x(nx-1) + dd * x(nx) return end