EVSL  1.1.0
EigenValues Slicing Library
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Macros
LapPLanR.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  ! max_it - Maximum number of iterations before a restart occurs
15  integer :: n, nx, ny, nz, nslices, mdeg, nvec, ev_int, mlan, nev, max_it
16 
17  ! Find the eigenvalues in the interval [a, b]
18  ! a - lower bound of the interval
19  ! b - upper bound of the interval
20  ! lmax - largest eigenvalue
21  ! lmin - smallest eigenvalue
22  ! thresh_int, thresh_ext - Polynomial variables
23  double precision :: a, b, lmax, lmin, tol, thresh_int, thresh_ext
24  ! sli - The spectrum slices
25  double precision, dimension(:), pointer :: sli
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  ! Set the A matrix in EVSL global data to point to the arrays built here
147  call evsl_arr2csr_f90(n, ia, ja, vals, csr)
148 
149  call evsl_seta_csr_f90(csr)
150 
151  ! kmpdos in EVSL for the DOS for dividing the spectrum
152  ! Set up necessary variabls for kpmdos
153  mdeg = 300;
154  nvec = 60;
155  allocate(sli(nslices+1))
156  ! Call EVSL kpmdos and spslicer
157  call evsl_kpm_spslicer_f90(mdeg, nvec, xintv, nslices, sli, ev_int)
158 
159  ! For each slice call ChebLatr
160  do i = 1, nslices
161  ! Prepare parameters for this slice
162  xintv(1) = sli(i)
163  xintv(2) = sli(i+1)
164  thresh_int = .5
165  thresh_ext = .15
166 
167  ! Call the EVSL function to create the polynomial
168  call evsl_find_pol_f90(xintv, thresh_int, thresh_ext, pol)
169 
170  ! Necessary paramters
171  nev = ev_int + 2
172  mlan = max(4*nev, 100)
173  mlan = min(mlan, n)
174  max_it = 3*mlan
175 
176  ! Call the EVSL cheblannr function to find the eigenvalues in the slice
177  call evsl_cheblantr_f90(mlan, nev, xintv, max_it, 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  write(*,*) nev, ' Eigs in this slice'
189 
190  ! Here you can do something with the eigenvalues that were found
191  ! The eigenvalues are stored in eigval and eigenvectors are in eigvec
192 
193  ! Be sure to deallocate the polynomial stored by EVSL
194  call evsl_free_pol_f90(pol)
195  deallocate(eigval)
196  deallocate(eigvec)
197  enddo
198  deallocate(sli)
199 
200  call evsl_free_csr_f90(csr)
201 
202  call evsl_finish_f90()
203 
204  ! Necessary Cleanup
205  deallocate(vals)
206  deallocate(ja)
207  deallocate(ia)
208 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