EVSL  1.1.0
EigenValues Slicing Library
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Macros
LapRLanR.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 - 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  ! max_it - Maximum number of iterations before a restart occurs
15  integer :: n, nx, ny, nz, nslices, mdeg, nvec, ev_int, mlan, nev, max_its
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  double precision :: a, b, lmax, lmin, tol
23  ! sli - The spectrum slices
24  double precision, dimension(:), pointer :: sli
25  double precision, dimension(4) :: xintv
26 
27  double precision, dimension(:), pointer :: eigval, eigvec
28  ! Matrix peices
29  double precision, dimension(:), pointer :: vals, valst
30  integer, dimension(:), pointer :: ia, iat
31  integer, dimension(:), pointer :: ja, jat
32 
33  integer*8 :: rat, asigbss, csr, dummy
34 
35  ! Loop varialbe declarations
36  integer :: i, j, k, zero
37 
38  ! DEBUG Variables
39  integer :: s, e
40 
41  ! Variables for reading command line arguments
42  character (len=32) :: arg
43  integer :: readerr
44 
45  ! Variables for using five point gen from SPARSKIT
46  double precision, dimension(:), pointer :: rhs
47  double precision, dimension(6) :: al
48  integer, dimension(:), pointer :: iau
49  integer :: mode
50 
51  ! Read in command line arguments to set important values.
52  ! The user can pass the phrase help to get the program to
53  ! print a usage statement then terminate
54 
55  zero = 0
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(valst(n*7))
121  allocate(ja(n*7)) !Size of number of nonzeros
122  allocate(jat(n*7))
123  allocate(ia(n+1)) !Size of number of rows + 1
124  allocate(iat(n+1))
125 
126  ! Allocate sparskit things
127  allocate(rhs(n*7)) ! Righthand side of size n
128  allocate(iau(n)) ! iau of size n
129  al(1) = 0.0d0; al(2) = 0.0d0;
130  al(3) = 0.0d0; al(4) = 0.0d0;
131  al(5) = 0.0d0; al(6) = 0.0d0;
132  mode = 0
133  call gen57pt(nx,ny,nz,al,mode,n,vals,ja,ia,iau,rhs)
134 
135  ! Conver to csc and back to csr to sort the indexing of the columns
136  call csrcsc(n, 1, 1, vals, ja, ia, valst, jat, iat)
137  call csrcsc(n, 1, 1, valst, jat, iat, vals, ja, ia)
138 
139  ! Since we're using this array with C we need to accomodate for C indexing
140  ia = ia - 1
141  ja = ja - 1
142  ! Cleanup extra sparskit information
143  deallocate(rhs)
144  deallocate(iau)
145  deallocate(valst)
146  deallocate(jat)
147  deallocate(iat)
148 
149  ! This section of the code will run the EVSL code.
150  ! This file is not utilizing the matrix free format and we'll pass
151  ! a CSR matrix in
152 
153  ! Initialize the EVSL global data
154  call evsl_start_f90()
155 
156  ! Set the A matrix in EVSL global data to point to the arrays built here
157  call evsl_arr2csr_f90(n, ia, ja, vals, csr)
158 
159  call evsl_seta_csr_f90(csr)
160 
161  ! kmpdos in EVSL for the DOS for dividing the spectrum
162  ! Set up necessary variables for kpmdos
163  mdeg = 300;
164  nvec = 60;
165  allocate(sli(nslices+1))
166  ! Call EVSL kpmdos and spslicer
167  call evsl_kpm_spslicer_f90(mdeg, nvec, xintv, nslices, sli, ev_int)
168 
169  ! For each slice call RatLanTr
170  do i = 1, nslices
171  ! Prepare parameters for this slice
172  xintv(1) = sli(i)
173  xintv(2) = sli(i+1)
174 
175  ! Call EVSL to create our rational filter
176  call evsl_find_rat_f90(xintv, rat)
177 
178  ! Get the factorizations of A-\sigma*B
179  ! matrix B is not present
180  call setup_asigmabsol_direct_f90(csr, zero, dummy, rat, asigbss)
181  ! Set the direct solver function
182  call set_asigmabsol_direct_f90(rat, asigbss)
183 
184  ! Necessary paramters
185  nev = ev_int + 2
186  mlan = max(4*nev, 100)
187  max_its = 3*mlan
188 
189  ! Call the EVSL ratlantr function to find the eigenvalues in the slice
190  call evsl_ratlantr_f90(mlan, nev, xintv, max_its, tol, rat)
191 
192  ! Extract the number of eigenvalues found from the EVSL global data
193  call evsl_get_nev_f90(nev)
194 
195  ! Allocate storage for the eigenvalue and vectors found from cheblannr
196  allocate(eigval(nev))
197  allocate(eigvec(nev*size(ia))) ! number of eigen values * number of rows
198 
199  ! Extract the arrays of eigenvalues and eigenvectors from the EVSL global data
200  call evsl_copy_result_f90(eigval, eigvec)
201  write(*,*) nev, ' Eigs in this slice'
202 
203  ! Here you can do something with the eigenvalues that were found
204  ! The eigenvalues are stored in eigval and eigenvectors are in eigvec
205 
206  ! Be sure to deallocate the rational filter as well as the solve data
207  call free_asigmabsol_direct_f90(rat, asigbss)
208  ! this must be done after FREE_ASIGMABSOL_DIRECT_F90
209  call evsl_free_rat_f90(rat)
210  deallocate(eigval)
211  deallocate(eigvec)
212  enddo
213  deallocate(sli)
214 
215  call evsl_free_csr_f90(csr)
216 
217  call evsl_finish_f90()
218 
219  ! Necessary Cleanup
220  deallocate(vals)
221  deallocate(ja)
222  deallocate(ia)
223 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 max(a, b)
Definition: def.h:56