EVSL  1.1.0 EigenValues Slicing Library
LapPLanN_MatFree.f90
Go to the documentation of this file.
1 program driver
2
3  use csrmatrix ! This is where our matrix class is stored
4
5  implicit none
6
7  ! Variable declarations
8  ! n - The number of rows in the matrix
9  ! nx, ny, nz - The dimension of the mesh to generate the laplacean from
10  ! nslices - The number of slices to divide the spectrum into
11  ! Mdeg - Polynomial degree
12  ! nvec - Number of sample vectors
13  ! ev_int - Number of eigenvalues per slice
14  ! mlan - Dimension of krylov subspace
15  ! nev - Approximate number of eigenvalues wanted
16  integer :: n, nx, ny, nz, nslices, mdeg, nvec, ev_int, mlan, nev
17
18  ! Find the eigenvalues in the interval [a, b]
19  ! a - lower bound of the interval
20  ! b - upper bound of the interval
21  ! lmax - largest eigenvalue
22  ! lmin - smallest eigenvalue
23  ! thresh_int, thresh_ext - Polynomial variables
24  double precision :: a, b, lmax, lmin, tol, thresh_int, thresh_ext
25  ! sli - The spectrum slices
26  double precision, dimension(:), pointer :: sli
27  double precision, dimension(4) :: xintv
28
29  double precision, dimension(:), pointer :: eigval, eigvec
30  ! Matrix peices
31  ! We build the matrix using the arrays vals, ia, ja and will then
32  ! use the csrmat type to point to these arrays.
33  type(csrmat) :: mat
34  double precision, dimension(:), pointer :: vals
35  integer, dimension(:), pointer :: ia
36  integer, dimension(:), pointer :: ja
37
38  ! Store the polynomial as an 8 byte integer
39  ! This will never be used on the fortran side, only passed to EVSL
40  integer*8 :: pol
41
42  ! Loop varialbe declarations
43  integer :: i, j, k
44
45  ! DEBUG Variables
46  integer :: s, e
47
48  ! Variables for reading command line arguments
49  character (len=32) :: arg
50  integer :: readerr
51
52  ! Variables for using five point gen from SPARSKIT
53  double precision, dimension(:), pointer :: rhs
54  double precision, dimension(6) :: al
55  integer, dimension(:), pointer :: iau
56  integer :: mode
57
58  ! Read in command line arguments to set important values.
59  ! The user can pass the phrase help to get the program to
60  ! print a usage statement then terminate
61
62  !Set default values
63  nx = 16
64  ny = 16
65  nz = 20
66  a = .40d0
67  b = .80d0
68  nslices = 4
69
70  !Crude but works.
71  ! This loop and if statements process the command line arguments.
72  ! Type ./LapPLanN.out help for usage information
73  do i = 1, iargc()
74  call getarg(i, arg)
75  arg = trim(arg)
76
77  if(arg(1:2) == 'nx') then
78  read(arg(4:), *, iostat = readerr) nx
79  elseif(arg(1:2) == 'ny') then
80  read(arg(4:), *, iostat = readerr) ny
81  elseif(arg(1:2) == 'nz') then
82  read(arg(4:), *, iostat = readerr) nz
83  elseif(arg(1:1) == 'a') then
84  read(arg(3:), *, iostat = readerr) a
85  elseif(arg(1:1) == 'b') then
86  read(arg(3:), *, iostat = readerr) b
87  elseif(arg(:7) == 'nslices') then
88  read(arg(9:), *, iostat = readerr) nslices
89  elseif(arg == 'help') then
90  write(*,*) 'Usage: ./testL.ex nx=[int] ny=[int] nz=[int] a=[double] b=[double] nslices=[int]'
91  stop
92  endif
93  if(readerr /= 0) then
94  write(*,*) 'There was an error while reading argument: ', arg
95  stop
96  endif
97  enddo
98
99  ! Initialize eigenvalue bounds set by hand
100  lmin = 0.0d0
101  if(nz == 1) then
102  lmax = 8.0d0
103  else
104  lmax = 12.0d0
105  endif
106  xintv(1) = a
107  xintv(2) = b
108  xintv(3) = lmin
109  xintv(4) = lmax
110  tol = 1e-08
111
112  ! Use SPARSKIT to generate the 2D/3D Laplacean matrix
113  ! We change the grid size to account for the boundaries that
114  ! SPARSKIT uses but not used by the LapGen tests in EVSL
115  nx = nx+2
116  if(ny > 1) then
117  ny = ny+2
118  if(nz > 1) then
119  nz = nz+2
120  endif
121  endif
122  n = nx*ny*nz
123
124  ! allocate our csr matrix
125  allocate(vals(n*7)) !Size of number of nonzeros
126  allocate(ja(n*7)) !Size of number of nonzeros
127  allocate(ia(n+1)) !Size of number of rows + 1
128
129  ! Allocate sparskit things
130  allocate(rhs(n*7)) ! Righthand side of size n
131  allocate(iau(n)) ! iau of size n
132  al(1) = 0.0d0; al(2) = 0.0d0;
133  al(3) = 0.0d0; al(4) = 0.0d0;
134  al(5) = 0.0d0; al(6) = 0.0d0;
135  mode = 0
136  call gen57pt(nx,ny,nz,al,mode,n,vals,ja,ia,iau,rhs)
137
138  ! Cleanup extra sparskit information
139  deallocate(rhs)
140  deallocate(iau)
141
142  ! Set up the csrmat object to point to the necessary arrays
143  mat%ia => ia
144  mat%ja => ja
145  mat%a => vals
146  mat%nrows = n
147  mat%ncols = n
148  mat%nnz = ia(n+1)
149
150  ! This section of the code will run the EVSL code.
151  ! Initialize the EVSL global data
152  call evsl_start_f90()
153
154  ! Since we are storing the matrix on the Fortran side we need to provide
155  ! EVSL with the matvec routine, and a pointer to our matrix data.
156  call evsl_setamv_f90(mat%nrows, csrmatvec, mat)
157
158  ! kmpdos in EVSL for the DOS for dividing the spectrum
159  ! Set up necessary variables for kpmdos
160  mdeg = 300;
161  nvec = 60;
162  allocate(sli(nslices+1))
163  ! Call EVSL kpmdos and spslicer
164  call evsl_kpm_spslicer_f90(mdeg, nvec, xintv, nslices, sli, ev_int)
165
166  ! For each slice call ChebLanr
167  do i = 1, nslices
168  ! Prepare parameters for this slice
169  xintv(1) = sli(i)
170  xintv(2) = sli(i+1)
171  thresh_int = .5
172  thresh_ext = .15
173
174  ! Call the EVSL function to create the polynomial
175  call evsl_find_pol_f90(xintv, thresh_int, thresh_ext, pol)
176
177  ! Necessary paramters
178  nev = ev_int + 2
179  mlan = max(4*nev, 100)
180  mlan = min(mlan, n)
181
182  ! Call the EVSL cheblannr function to find the eigenvalues in the slice
183  call evsl_cheblannr_f90(xintv, mlan, tol, pol)
184
185  ! Extract the number of eigenvalues found from the EVSL global data
186  call evsl_get_nev_f90(nev)
187
188  ! Allocate storage for the eigenvalue and vectors found from cheblannr
189  allocate(eigval(nev))
190  allocate(eigvec(nev*size(ia))) ! number of eigen values * number of rows
191
192  ! Extract the arrays of eigenvalues and eigenvectors from the EVSL global data
193  call evsl_copy_result_f90(eigval, eigvec)
194  write(*,*) nev, ' Eigs in this slice'
195
196  ! Here you can do something with the eigenvalues that were found
197  ! The eigenvalues are stored in eigval and eigenvectors are in eigvec
198
199  ! Be sure to deallocate the polynomial stored by EVSL
200  call evsl_free_pol_f90(pol)
201  deallocate(eigval)
202  deallocate(eigvec)
203  enddo
204  deallocate(sli)
205
206  call evsl_finish_f90()
207
208  ! Necessary Cleanup
209  deallocate(vals)
210  deallocate(ja)
211  deallocate(ia)
212 end program driver
subroutine gen57pt(nx, ny, nz, al, mode, n, a, ja, ia, iau, rhs)
Definition: genmat.f:20
subroutine csrmatvec(x, y, mat)
Definition: csr_module.f90:21
program driver
Definition: LapPLanN.f90:1
#define min(a, b)
Definition: def.h:62
#define max(a, b)
Definition: def.h:56