EVSL  1.1.0 EigenValues Slicing Library
LapRLanN.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 dimensino of the mesh to generate the laplacean from
8  ! nslices - The number of slices to divide the spectrum into
9  ! Mdeg - Degree for kpmdos
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  double precision :: a, b, lmax, lmin, tol
22  ! sli - The spectrum slices
23  double precision, dimension(:), pointer :: sli
24  double precision, dimension(4) :: xintv
25
26  double precision, dimension(:), pointer :: eigval, eigvec
27  ! Matrix peices
28  double precision, dimension(:), pointer :: vals, valst
29  integer, dimension(:), pointer :: ia, iat
30  integer, dimension(:), pointer :: ja, jat
31
32  integer*8 :: rat, asigbss, csr, dummy
33
34  ! Loop varialbe declarations
35  integer :: i, j, k, zero
36
37  ! DEBUG Variables
38  integer :: s, e
39
40  ! Variables for reading command line arguments
41  character (len=32) :: arg
43
44  ! Variables for using five point gen from SPARSKIT
45  double precision, dimension(:), pointer :: rhs
46  double precision, dimension(6) :: al
47  integer, dimension(:), pointer :: iau
48  integer :: mode
49
50  ! Read in command line arguments to set important values.
51  ! The user can pass the phrase help to get the program to
52  ! print a usage statement then terminate
53
54  zero = 0
55  !Set default values
56  nx = 16
57  ny = 16
58  nz = 20
59  a = .40d0
60  b = .80d0
61  nslices = 4
62
63  !Crude but works.
64  ! This loop and if statements process the command line arguments.
65  ! Type ./LapPLanN.out help for usage information
66  do i = 1, iargc()
67  call getarg(i, arg)
68  arg = trim(arg)
69
70  if(arg(1:2) == 'nx') then
72  elseif(arg(1:2) == 'ny') then
74  elseif(arg(1:2) == 'nz') then
76  elseif(arg(1:1) == 'a') then
78  elseif(arg(1:1) == 'b') then
80  elseif(arg(:7) == 'nslices') then
82  elseif(arg == 'help') then
83  write(*,*) 'Usage: ./testL.ex nx=[int] ny=[int] nz=[int] a=[double] b=[double] nslices=[int]'
84  stop
85  endif
87  write(*,*) 'There was an error while reading argument: ', arg
88  stop
89  endif
90  enddo
91
92  ! Initialize eigenvalue bounds set by hand
93  lmin = 0.0d0
94  if(nz == 1) then
95  lmax = 8.0d0
96  else
97  lmax = 12.0d0
98  endif
99  xintv(1) = a
100  xintv(2) = b
101  xintv(3) = lmin
102  xintv(4) = lmax
103  tol = 1e-08
104
105  ! Use SPARSKIT to generate the 2D/3D Laplacean matrix
106  ! We change the grid size to account for the boundaries that
107  ! SPARSKIT uses but not used by the LapGen tests in EVSL
108  nx = nx+2
109  if(ny > 1) then
110  ny = ny+2
111  if(nz > 1) then
112  nz = nz+2
113  endif
114  endif
115  n = nx*ny*nz
116
117  ! allocate our csr matrix
118  allocate(vals(n*7)) !Size of number of nonzeros
119  allocate(valst(n*7))
120  allocate(ja(n*7)) !Size of number of nonzeros
121  allocate(jat(n*7))
122  allocate(ia(n+1)) !Size of number of rows + 1
123  allocate(iat(n+1))
124
125  ! Allocate sparskit things
126  allocate(rhs(n*7)) ! Righthand side of size n
127  allocate(iau(n)) ! iau of size n
128  al(1) = 0.0d0; al(2) = 0.0d0;
129  al(3) = 0.0d0; al(4) = 0.0d0;
130  al(5) = 0.0d0; al(6) = 0.0d0;
131  mode = 0
132  call gen57pt(nx,ny,nz,al,mode,n,vals,ja,ia,iau,rhs)
133
134  ! Conver to csc and back to csr to sort the indexing of the columns
135  call csrcsc(n, 1, 1, vals, ja, ia, valst, jat, iat)
136  call csrcsc(n, 1, 1, valst, jat, iat, vals, ja, ia)
137
138  ! Since we're using this array with C we need to accomodate for C indexing
139  ia = ia - 1
140  ja = ja - 1
141  ! Cleanup extra sparskit information
142  deallocate(rhs)
143  deallocate(iau)
144  deallocate(valst)
145  deallocate(jat)
146  deallocate(iat)
147
148  ! This section of the code will run the EVSL code.
149  ! This file is not utilizing the matrix free format and we'll pass
150  ! a CSR matrix in
151
152  ! Initialize the EVSL global data
153  call evsl_start_f90()
154
155  ! Set the A matrix in EVSL global data to point to the arrays built here
156  call evsl_arr2csr_f90(n, ia, ja, vals, csr)
157
158  call evsl_seta_csr_f90(csr)
159
160  ! kmpdos in EVSL for the DOS for dividing the spectrum
161  ! Set up necessary variables for kpmdos
162  mdeg = 300;
163  nvec = 60;
164  allocate(sli(nslices+1))
165  ! Call EVSL kpmdos and spslicer
166  call evsl_kpm_spslicer_f90(mdeg, nvec, xintv, nslices, sli, ev_int)
167
168  ! For each slice call RatLanNr
169  do i = 1, nslices
170  ! Prepare parameters for this slice
171  xintv(1) = sli(i)
172  xintv(2) = sli(i+1)
173
174  ! Call EVSL to create our rational filter
175  call evsl_find_rat_f90(xintv, rat)
176
177  ! Get the factorizations of A-\sigma*B
178  ! matrix B is not present
179  call setup_asigmabsol_direct_f90(csr, zero, dummy, rat, asigbss)
180  ! Set the direct solver function
181  call set_asigmabsol_direct_f90(rat, asigbss)
182
183  ! Necessary paramters
184  nev = ev_int + 2
185  mlan = max(4*nev, 100)
186  mlan = min(mlan, n)
187
188  ! Call the EVSL ratlannr function to find the eigenvalues in the slice
189  call evsl_ratlannr_f90(xintv, mlan, tol, rat)
190
191  ! Extract the number of eigenvalues found from the EVSL global data
192  call evsl_get_nev_f90(nev)
193
194  ! Allocate storage for the eigenvalue and vectors found from cheblannr
195  allocate(eigval(nev))
196  allocate(eigvec(nev*size(ia))) ! number of eigen values * number of rows
197
198  ! Extract the arrays of eigenvalues and eigenvectors from the EVSL global data
199  call evsl_copy_result_f90(eigval, eigvec)
200  write(*,*) nev, ' Eigs in this slice'
201
202  ! Here you can do something with the eigenvalues that were found
203  ! The eigenvalues are stored in eigval and eigenvectors are in eigvec
204
205  ! Be sure to deallocate the rational filter as well as the solve data
206  call free_asigmabsol_direct_f90(rat, asigbss)
207  ! this must be done after FREE_ASIGMABSOL_DIRECT_F90
208  call evsl_free_rat_f90(rat)
209  deallocate(eigval)
210  deallocate(eigvec)
211  enddo
212  deallocate(sli)
213
214  call evsl_free_csr_f90(csr)
215
216  call evsl_finish_f90()
217
218  ! Necessary Cleanup
219  deallocate(vals)
220  deallocate(ja)
221  deallocate(ia)
222 end program driver
subroutine gen57pt(nx, ny, nz, al, mode, n, a, ja, ia, iau, rhs)
Definition: genmat.f:20
void csrcsc(int OUTINDEX, const int nrow, const int ncol, int job, double *a, int *ja, int *ia, double *ao, int *jao, int *iao)
convert csr to csc Assume input csr is 0-based index output csc 0/1 index specified by OUTINDEX * ...
Definition: spmat.c:17
program driver
Definition: LapPLanN.f90:1
#define min(a, b)
Definition: def.h:62
#define max(a, b)
Definition: def.h:56