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