FILE ./indexer.f90 -------------- MODULE indexer ============== ! Indexer class. Useful to flatten a multi-dimensional indexed expansion ! onto a flat data structure. implicit none character(*), parameter, private :: version = '$Id: doc.txt,v 1.1 2003/12/18 21:25:49 pletzer Exp $' type indexer_obj integer counter ! flat array index integer ntot ! flat array total size integer ndims ! no of dimensions ! the following contain the sizes and starting indices ! in each dimension integer, pointer :: sizes(:), start(:) end type indexer_obj contains SUBROUTINE indexer_init(this, sizes, start, ier) SUBROUTINE indexer_free(this, ier) SUBROUTINE indexer_next(this, ier) SUBROUTINE indexer_get(this, indices, ier) SUBROUTINE indexer_save(this, filename, ier) SUBROUTINE indexer_load(this, filename, ier) SUBROUTINE indexer_error(this, nunit, ier) SUBROUTINE indexer_test(this, ier) == in module indexer == ------------------------------------------------ SUBROUTINE indexer_init(this, sizes, start, ier) ================================================ ! Constructor. type(indexer_obj) :: this integer, intent(inout) :: sizes(:) ! An array containing the sizes integer, intent(in) :: start(:) integer, intent(out) :: ier ------------------------------------------------ ---------------------------------- SUBROUTINE indexer_free(this, ier) ================================== ! Destructor type(indexer_obj) :: this integer, intent(out) :: ier ---------------------------------- ---------------------------------- SUBROUTINE indexer_next(this, ier) ================================== ! Increment counter type(indexer_obj) :: this integer, intent(out) :: ier ---------------------------------- ------------------------------------------ SUBROUTINE indexer_get(this, indices, ier) ========================================== ! Return sub-indices type(indexer_obj) :: this integer, intent(out) :: indices(:) integer, intent(out) :: ier integer j1, j2, periodicity(this % ndims) ------------------------------------------ -------------------------------------------- SUBROUTINE indexer_save(this, filename, ier) ============================================ USE ezcdf ! Save state in netCDF file 'filename' type(indexer_obj) :: this character(*), intent(in) :: filename integer, intent(out) :: ier -------------------------------------------- -------------------------------------------- SUBROUTINE indexer_load(this, filename, ier) ============================================ USE ezcdf ! Load state from netCDF file 'filename' type(indexer_obj) :: this character(*), intent(in) :: filename integer, intent(out) :: ier -------------------------------------------- ------------------------------------------ SUBROUTINE indexer_error(this, nunit, ier) ========================================== ! Print error message type(indexer_obj) :: this integer, intent(in) :: nunit ! output channel number integer, intent(in) :: ier ------------------------------------------ ---------------------------------- SUBROUTINE indexer_test(this, ier) ================================== ! Test unit type(indexer_obj) :: this integer, intent(out) :: ier ---------------------------------- FILE ./airy.f90 ----------- MODULE airy =========== ! Module for Airy test problem use supralu_mod use indexer use Gabor use kzgrid use xgrid use operators use plotmtv implicit none integer, parameter, private :: r8 = selected_real_kind(12,100) character(*), parameter, private :: version = '$Id: doc.txt,v 1.1 2003/12/18 21:25:49 pletzer Exp $' integer, parameter :: Nd = 1 ! space dimension integer, parameter :: Nc = 1 ! no of components type airy_obj integer :: Nz ! No of Gabor intervals integer :: Nj ! Fourier modes -Nj .. + Nj real(r8) :: sampling ! sampling rate > 1 real(r8) :: alpha ! equation parameter real(r8) :: tol ! sparsity tolerance complex(r8) :: bc0, bc1 ! Boundary conditions real(r8) :: xmin, xmax ! domain boundary real(r8) :: eps real(r8) :: sigma integer :: Nx type(sparseC16_obj) :: amat type(indexer_obj) :: iter2 type(kzgrid_obj) :: kzg complex(r8), dimension(:), pointer :: xvec, bvec end type airy_obj contains SUBROUTINE airy_init(this, Nz, Nj, sigma_dx2, sampling) SUBROUTINE airy_free(this) SUBROUTINE airy_set_domain(this, xmin, xmax) SUBROUTINE airy_set_tol(this, tol) SUBROUTINE airy_set_equation(this, alpha) SUBROUTINE airy_set_bcs(this, bc0, bc1) SUBROUTINE airy_assemble(this) SUBROUTINE airy_solve(this) SUBROUTINE airy_get_solution(this, x, y) SUBROUTINE airy_plot(this, x, y) SUBROUTINE airy_test(this) == in module airy == ------------------------------------------------------- SUBROUTINE airy_init(this, Nz, Nj, sigma_dx2, sampling) ======================================================= type(airy_obj) :: this integer, intent(in) :: Nz ! No of Gabor intervals integer, intent(in) :: Nj ! Fourier modes -Nj .. + Nj real(r8), intent(in) :: sigma_dx2 ! inverse Gabor extent normalized to grid size real(r8), intent(in) :: sampling ! sampling rate > 1 ------------------------------------------------------- -------------------------- SUBROUTINE airy_free(this) ========================== type(airy_obj) :: this -------------------------- -------------------------------------------- SUBROUTINE airy_set_domain(this, xmin, xmax) ============================================ type(airy_obj) :: this real(r8), intent(in) :: xmin, xmax ! domain boundary -------------------------------------------- ---------------------------------- SUBROUTINE airy_set_tol(this, tol) ================================== type(airy_obj) :: this real(r8), intent(in) :: tol ! sparsity tolerance ! A low value of tol will make the matrices less sparse ---------------------------------- ----------------------------------------- SUBROUTINE airy_set_equation(this, alpha) ========================================= type(airy_obj) :: this real(r8), intent(in) :: alpha ! equation parameter ----------------------------------------- --------------------------------------- SUBROUTINE airy_set_bcs(this, bc0, bc1) ======================================= type(airy_obj) :: this complex(r8), intent(in) :: bc0, bc1 ! Dirichlet boundary conditions --------------------------------------- ------------------------------ SUBROUTINE airy_assemble(this) ============================== type(airy_obj) :: this ! In order to exercize iterators, the matrix ! is assembled using the general, multi-dimensional ! iterator based method, although for a 1-d ! equation it might be more efficient to use do-loops. ------------------------------ --------------------------- SUBROUTINE airy_solve(this) =========================== type(airy_obj) :: this --------------------------- ---------------------------------------- SUBROUTINE airy_get_solution(this, x, y) ======================================== type(airy_obj) :: this real(r8), intent(in) :: x(:) ! abscissae complex(r8), intent(out) :: y(:) ! ordinates ---------------------------------------- -------------------------------- SUBROUTINE airy_plot(this, x, y) ================================ type(airy_obj) :: this real(r8), intent(in) :: x(:) ! abscissae complex(r8), intent(in) :: y(:) ! ordinates -------------------------------- -------------------------- SUBROUTINE airy_test(this) ========================== ! Test module type(airy_obj) :: this -------------------------- FILE ./tester.f90 ------------- MODULE tester ============= == in module tester == FILE ./coupled.f90 -------------- MODULE coupled ============== ! 4-th order coupled wave test problem use supralu_mod use indexer use Gabor use kzgrid use xgrid use operators use plotmtv implicit none integer, parameter, private :: r8 = selected_real_kind(12,100) character(*), parameter, private :: version = '$Id: doc.txt,v 1.1 2003/12/18 21:25:49 pletzer Exp $' integer, parameter :: Nd = 1 ! space dimension integer, parameter :: Nc = 1 ! no of components type coupled_obj integer :: Nz ! No of Gabor intervals integer :: Nj ! Fourier modes -Nj .. + Nj real(r8) :: sampling ! sampling rate > 1 real(r8) :: k1, k2, alpha, beta1, beta2 ! equation parameters real(r8) :: tol ! sparsity tolerance complex(r8) :: bc0, bc1, bcp0, bcp1 ! Boundary conditions real(r8) :: xmin, xmax ! domain boundary real(r8) :: eps real(r8) :: sigma integer :: Nx type(sparseC16_obj) :: amat type(indexer_obj) :: iter2 type(kzgrid_obj) :: kzg complex(r8), dimension(:), pointer :: xvec, bvec end type coupled_obj contains SUBROUTINE coupled_init(this, Nz, Nj, sigma_dx2, sampling) SUBROUTINE coupled_free(this) FUNCTION coupled_k1sq_function(this, x) FUNCTION coupled_k2sq_function(this, x) FUNCTION coupled_k2sq_prime_function(this, x) SUBROUTINE coupled_set_domain(this, xmin, xmax) SUBROUTINE coupled_set_tol(this, tol) SUBROUTINE coupled_set_equation(this, alpha, beta1, beta2, k1, k2) SUBROUTINE coupled_set_bcs(this, bc, bcp) SUBROUTINE coupled_assemble(this) SUBROUTINE coupled_solve(this) SUBROUTINE coupled_get_solution(this, x, y) SUBROUTINE coupled_plot(this, x, y, label) SUBROUTINE coupled_test(this) == in module coupled == ---------------------------------------------------------- SUBROUTINE coupled_init(this, Nz, Nj, sigma_dx2, sampling) ========================================================== type(coupled_obj) :: this integer, intent(in) :: Nz ! No of Gabor intervals integer, intent(in) :: Nj ! Fourier modes -Nj .. + Nj real(r8), intent(in) :: sigma_dx2 ! inverse Gabor extent normalized to grid size real(r8), intent(in) :: sampling ! sampling rate > 1 ---------------------------------------------------------- ----------------------------- SUBROUTINE coupled_free(this) ============================= type(coupled_obj) :: this ----------------------------- --------------------------------------- FUNCTION coupled_k1sq_function(this, x) ======================================= type(coupled_obj) :: this real(r8), intent(in) :: x --------------------------------------- --------------------------------------- FUNCTION coupled_k2sq_function(this, x) ======================================= type(coupled_obj) :: this real(r8), intent(in) :: x --------------------------------------- --------------------------------------------- FUNCTION coupled_k2sq_prime_function(this, x) ============================================= type(coupled_obj) :: this real(r8), intent(in) :: x --------------------------------------------- ----------------------------------------------- SUBROUTINE coupled_set_domain(this, xmin, xmax) =============================================== type(coupled_obj) :: this real(r8), intent(in) :: xmin, xmax ! domain boundary ----------------------------------------------- ------------------------------------- SUBROUTINE coupled_set_tol(this, tol) ===================================== type(coupled_obj) :: this real(r8), intent(in) :: tol ! sparsity tolerance ! A low value of tol will make the matrices less sparse ------------------------------------- ------------------------------------------------------------------ SUBROUTINE coupled_set_equation(this, alpha, beta1, beta2, k1, k2) ================================================================== type(coupled_obj) :: this real(r8), intent(in) :: alpha, beta1, beta2, k1, k2 ! equation parameters ------------------------------------------------------------------ ----------------------------------------- SUBROUTINE coupled_set_bcs(this, bc, bcp) ========================================= type(coupled_obj) :: this complex(r8), intent(in) :: bc(2), bcp(2) ! boundary conditions ----------------------------------------- --------------------------------- SUBROUTINE coupled_assemble(this) ================================= type(coupled_obj) :: this ! In order to exercize iterators, the matrix ! is assembled using the general, multi-dimensional ! iterator based method, although for a 1-d ! equation it might be more efficient to use do-loops. --------------------------------- ------------------------------ SUBROUTINE coupled_solve(this) ============================== type(coupled_obj) :: this ------------------------------ ------------------------------------------- SUBROUTINE coupled_get_solution(this, x, y) =========================================== type(coupled_obj) :: this real(r8), intent(in) :: x(:) ! abscissae complex(r8), intent(out) :: y(:) ! ordinates ------------------------------------------- ------------------------------------------ SUBROUTINE coupled_plot(this, x, y, label) ========================================== type(coupled_obj) :: this real(r8), intent(in) :: x(:) ! abscissae complex(r8), intent(in) :: y(:) ! ordinates character(*), intent(in) :: label ------------------------------------------ ----------------------------- SUBROUTINE coupled_test(this) ============================= ! Test module type(coupled_obj) :: this ----------------------------- FILE ./gabor.f90 ------------ MODULE gabor ============ integer, parameter, private :: r8 = selected_real_kind(12,100) complex(r8), parameter, private :: imag = (0._r8, 1._r8) character(*), parameter, private :: version = '$Id: doc.txt,v 1.1 2003/12/18 21:25:49 pletzer Exp $' contains FUNCTION Gabor0(x, z, ak, sigma) FUNCTION GaborN(x, z, ak, sigma) == in module gabor == -------------------------------- FUNCTION Gabor0(x, z, ak, sigma) ================================ ! Scalar Gabor function real(r8), intent(in) :: x ! collocation point real(r8), intent(in) :: z ! envelop position real(r8), intent(in) :: ak ! wave number real(r8), intent(in) :: sigma ! inverse extent -------------------------------- -------------------------------- FUNCTION GaborN(x, z, ak, sigma) ================================ ! N-dimensional Gabor function real(r8), intent(in) :: x(:) ! collocation point real(r8), intent(in) :: z(:) ! envelop position real(r8), intent(in) :: ak(:) ! wave number real(r8), intent(in) :: sigma(:) ! inverse extent -------------------------------- FILE ./ktae.f90 ----------- MODULE ktae =========== ! Module for Eigenvalue Bessel test problem use supralu_mod use indexer use Gabor use kzgrid use xgrid use operators use plotmtv implicit none integer, parameter, private :: r8 = selected_real_kind(12,100) character(*), parameter, private :: version = '$Id: doc.txt,v 1.1 2003/12/18 21:25:49 pletzer Exp $' integer, parameter :: Nd = 1 ! space dimension integer, parameter :: Nc = 1 ! no of components type ktae_obj integer :: Nz ! No of Gabor intervals integer :: Nj ! Fourier modes -Nj .. + Nj real(r8) :: sampling ! sampling rate > 1 real(r8) :: tol ! sparsity tolerance real(r8) :: xmin, xmax ! domain boundary ! equation parameters integer :: mpol, ntor ! poloidal/toroidal mode number real(r8) :: q0, q1 ! axis/edgesafety factor values real(r8) :: rhobar, rho_delta ! gyroradius real(r8) :: eps real(r8) :: sigma integer :: Nx type(sparseC16_obj) :: amat, bmat type(indexer_obj) :: iter2 type(kzgrid_obj) :: kzg end type ktae_obj contains SUBROUTINE ktae_init(this, Nz, Nj, sigma_dx2, sampling) SUBROUTINE ktae_free(this) SUBROUTINE ktae_set_domain(this, xmin, xmax) SUBROUTINE ktae_set_tol(this, tol) SUBROUTINE ktae_assemble(this) SUBROUTINE ktae_solve(this, xvec, lambda, nmax) SUBROUTINE ktae_get_solution(this, xvec, x, y) SUBROUTINE ktae_plot(this, x, y, title) SUBROUTINE ktae_test(this) FUNCTION q(this, x) FUNCTION q_prime(this, x) FUNCTION K_square(this, x, m) FUNCTION K_square_prime(this, x, m) == in module ktae == ------------------------------------------------------- SUBROUTINE ktae_init(this, Nz, Nj, sigma_dx2, sampling) ======================================================= type(ktae_obj) :: this integer, intent(in) :: Nz ! No of Gabor intervals integer, intent(in) :: Nj ! Fourier modes -Nj .. + Nj real(r8), intent(in) :: sigma_dx2 ! inverse Gabor extent normalized to grid size real(r8), intent(in) :: sampling ! sampling rate > 1 ------------------------------------------------------- -------------------------- SUBROUTINE ktae_free(this) ========================== type(ktae_obj) :: this -------------------------- -------------------------------------------- SUBROUTINE ktae_set_domain(this, xmin, xmax) ============================================ type(ktae_obj) :: this real(r8), intent(in) :: xmin, xmax ! domain boundary -------------------------------------------- ---------------------------------- SUBROUTINE ktae_set_tol(this, tol) ================================== type(ktae_obj) :: this real(r8), intent(in) :: tol ! sparsity tolerance ! A low value of tol will make the matrices less sparse ---------------------------------- ------------------------------ SUBROUTINE ktae_assemble(this) ============================== type(ktae_obj) :: this ! In order to exercize iterators, the matrix ! is assembled using the general, multi-dimensional ! iterator based method, although for a 1-d ! equation it might be more efficient to use do-loops. ------------------------------ ----------------------------------------------- SUBROUTINE ktae_solve(this, xvec, lambda, nmax) =============================================== type(ktae_obj) :: this complex(r8), intent(inout) :: xvec(:) ! initial eigenvector guess on input/eigenvector on output complex(r8), intent(inout) :: lambda ! initial eigenvalue guess on input/eigenvalue on output integer, intent(inout) :: nmax ! max number of iterations on input/no of iterations on output ----------------------------------------------- ---------------------------------------------- SUBROUTINE ktae_get_solution(this, xvec, x, y) ============================================== type(ktae_obj) :: this complex(r8), intent(in) :: xvec(:) ! eigenvector returned by solve method real(r8), intent(in) :: x(:) ! abscissae complex(r8), intent(out) :: y(:) ! ordinates ---------------------------------------------- --------------------------------------- SUBROUTINE ktae_plot(this, x, y, title) ======================================= type(ktae_obj) :: this real(r8), intent(in) :: x(:) ! abscissae complex(r8), intent(in) :: y(:) ! ordinates character(*), intent(in) :: title --------------------------------------- -------------------------- SUBROUTINE ktae_test(this) ========================== USE supralu_mod ! Test module type(ktae_obj) :: this -------------------------- ------------------- FUNCTION q(this, x) =================== ! safety factor profile type(ktae_obj) :: this real(r8), intent(in) :: x ------------------- ------------------------- FUNCTION q_prime(this, x) ========================= ! derivative of safety factor profile type(ktae_obj) :: this real(r8), intent(in) :: x ------------------------- ----------------------------- FUNCTION K_square(this, x, m) ============================= type(ktae_obj) :: this real(r8), intent(in) :: x integer, intent(in) :: m ----------------------------- ----------------------------------- FUNCTION K_square_prime(this, x, m) =================================== type(ktae_obj) :: this real(r8), intent(in) :: x integer, intent(in) :: m ----------------------------------- FILE ./xkzgrid.f90 ------------- MODULE kzgrid ============= ! Handles the conversion of indices returned by the indexer into ! ak (=wavenumber) and z (=Gabor positons). implicit none integer, parameter, private :: r8 = selected_real_kind(12,100) real(r8), parameter, private :: twopi = 6.2831853071795865_r8 character(*), parameter, private :: version = '$Id: doc.txt,v 1.1 2003/12/18 21:25:49 pletzer Exp $' type kzgrid_obj integer ndims integer, pointer :: L(:), J(:) real(r8), pointer :: xmin(:), xmax(:) real(r8), pointer :: delta_k(:), delta_z(:) real(r8), pointer :: sampling(:) integer, pointer :: map_kz_indices(:) end type kzgrid_obj contains SUBROUTINE kzgrid_init(this, xmin, xmax, L, J, sampling, map_kz_indices, ier) SUBROUTINE kzgrid_free(this, ier) SUBROUTINE kzgrid_get(this, indices, ak, z, ier) SUBROUTINE kzgrid_test(this, ier) == in module kzgrid == ----------------------------------------------------------------------------- SUBROUTINE kzgrid_init(this, xmin, xmax, L, J, sampling, map_kz_indices, ier) ============================================================================= ! Constructor type(kzgrid_obj) :: this ! min/max points in the domain real(r8), intent(in) :: xmin(:), xmax(:) integer, intent(in) :: L(:) ! L+1 envelopes integer, intent(in) :: J(:) ! 2*J+1 Fourier modes real(r8), intent(in) :: sampling(:) ! sampling rate should be > 1 ! map_kz_indices maps the set of indices returned by the indexer into ! the indices of k and z (in that order). This provides the flexibility ! the reorder the memory layout for best performance. ! ! Example: in 3D, a possible layout could be map_kz_indices = (/2,3,4, 5,6,7/), ! ie the 2nd, 3rd, and 4th indices give the k-vector while indices 5, 6, and 7 ! yield the Gabor positions. In this case index 1 would be associated ! with field components, for instance. If one rather prefers to iterate first over ! the k-vector and then over field components, set map_kz_indices = (/1,2,3, 5,6,7/). integer, intent(in) :: map_kz_indices(:) integer, intent(out) :: ier ----------------------------------------------------------------------------- --------------------------------- SUBROUTINE kzgrid_free(this, ier) ================================= ! Destructor type(kzgrid_obj) :: this integer, intent(out) :: ier --------------------------------- ------------------------------------------------ SUBROUTINE kzgrid_get(this, indices, ak, z, ier) ================================================ ! Given a set of indices, return the components of ! x, k, and z. type(kzgrid_obj) :: this integer, intent(in) :: indices(:) ! from indexer real(r8), intent(out) :: ak(:), z(:) ! Wavenumber and Gabor position integer, intent(out) :: ier ------------------------------------------------ --------------------------------- SUBROUTINE kzgrid_test(this, ier) ================================= ! Tester type(kzgrid_obj) :: this integer, intent(out) :: ier ! computational domain --------------------------------- ------------ MODULE xgrid ============ ! Handles the conversion of indices returned by the indexer into ! x (= collocation) positions. implicit none integer, parameter, private :: r8 = selected_real_kind(12,100) real(r8), parameter, private :: twopi = 6.2831853071795865_r8 character(*), parameter, private :: version = '$Id: doc.txt,v 1.1 2003/12/18 21:25:49 pletzer Exp $' type xgrid_obj integer ndims integer, pointer :: M(:) ! no of intervals real(r8), pointer :: xmin(:), xmax(:) real(r8), pointer :: delta_x(:) integer, pointer :: map_x_indices(:) end type xgrid_obj contains SUBROUTINE xgrid_init(this, xmin, xmax, M, map_x_indices, ier) SUBROUTINE xgrid_free(this, ier) SUBROUTINE xgrid_get(this, indices, x, ier) SUBROUTINE xgrid_test(this, ier) == in module xgrid == -------------------------------------------------------------- SUBROUTINE xgrid_init(this, xmin, xmax, M, map_x_indices, ier) ============================================================== ! Constructor type(xgrid_obj) :: this ! min/max points in the domain real(r8), intent(in) :: xmin(:), xmax(:) integer, intent(in) :: M(:) ! M+1 collocation points ! map_x_indices maps the set of indices returned by the indexer into ! the indices of collocation points. This provides the flexibility ! the reorder the memory layout for best performance. ! ! Example: in 3D, a possible layout could be map_x_indices = (/2,3,4/), ! ie the 2nd, 3rd and 4th indices give the x-point. In this case index 1 ! would be associated with field components, for instance. If one rather prefers ! to iterate first over x-components and then over field components, set ! map_x_indices = (/1,2,3/). integer, intent(in) :: map_x_indices(:) integer, intent(out) :: ier -------------------------------------------------------------- -------------------------------- SUBROUTINE xgrid_free(this, ier) ================================ ! Destructor type(xgrid_obj) :: this integer, intent(out) :: ier -------------------------------- ------------------------------------------- SUBROUTINE xgrid_get(this, indices, x, ier) =========================================== ! Given a set of indices, return the components of ! x. type(xgrid_obj) :: this integer, intent(in) :: indices(:) ! from indexer real(r8), intent(out) :: x(:) ! collocation position integer, intent(out) :: ier ------------------------------------------- -------------------------------- SUBROUTINE xgrid_test(this, ier) ================================ ! Tester type(xgrid_obj) :: this integer, intent(out) :: ier ! computational domain -------------------------------- FILE ./order4.f90 ------------- MODULE order4 ============= ! 4-th order test problem use supralu_mod use indexer use Gabor use kzgrid use xgrid use operators use plotmtv implicit none integer, parameter, private :: r8 = selected_real_kind(12,100) character(*), parameter, private :: version = '$Id: doc.txt,v 1.1 2003/12/18 21:25:49 pletzer Exp $' integer, parameter :: Nd = 1 ! space dimension integer, parameter :: Nc = 1 ! no of components type order4_obj integer :: Nz ! No of Gabor intervals integer :: Nj ! Fourier modes -Nj .. + Nj real(r8) :: sampling ! sampling rate > 1 real(r8) :: k1, k2, alpha ! equation parameters real(r8) :: tol ! sparsity tolerance complex(r8) :: bc0, bc1, bcp0, bcp1 ! Boundary conditions real(r8) :: xmin, xmax ! domain boundary real(r8) :: eps real(r8) :: sigma integer :: Nx type(sparseC16_obj) :: amat type(indexer_obj) :: iter2 type(kzgrid_obj) :: kzg complex(r8), dimension(:), pointer :: xvec, bvec end type order4_obj contains SUBROUTINE order4_init(this, Nz, Nj, sigma_dx2, sampling) SUBROUTINE order4_free(this) SUBROUTINE order4_set_domain(this, xmin, xmax) SUBROUTINE order4_set_tol(this, tol) SUBROUTINE order4_set_equation(this, alpha, k1, k2) SUBROUTINE order4_set_bcs(this, bc, bcp) SUBROUTINE order4_assemble(this) SUBROUTINE order4_solve(this) SUBROUTINE order4_get_solution(this, x, y) SUBROUTINE order4_plot(this, x, y, label) SUBROUTINE order4_test(this) == in module order4 == --------------------------------------------------------- SUBROUTINE order4_init(this, Nz, Nj, sigma_dx2, sampling) ========================================================= type(order4_obj) :: this integer, intent(in) :: Nz ! No of Gabor intervals integer, intent(in) :: Nj ! Fourier modes -Nj .. + Nj real(r8), intent(in) :: sigma_dx2 ! inverse Gabor extent normalized to grid size real(r8), intent(in) :: sampling ! sampling rate > 1 --------------------------------------------------------- ---------------------------- SUBROUTINE order4_free(this) ============================ type(order4_obj) :: this ---------------------------- ---------------------------------------------- SUBROUTINE order4_set_domain(this, xmin, xmax) ============================================== type(order4_obj) :: this real(r8), intent(in) :: xmin, xmax ! domain boundary ---------------------------------------------- ------------------------------------ SUBROUTINE order4_set_tol(this, tol) ==================================== type(order4_obj) :: this real(r8), intent(in) :: tol ! sparsity tolerance ! A low value of tol will make the matrices less sparse ------------------------------------ --------------------------------------------------- SUBROUTINE order4_set_equation(this, alpha, k1, k2) =================================================== type(order4_obj) :: this real(r8), intent(in) :: alpha, k1, k2 ! equation parameters --------------------------------------------------- ---------------------------------------- SUBROUTINE order4_set_bcs(this, bc, bcp) ======================================== type(order4_obj) :: this complex(r8), intent(in) :: bc(2), bcp(2) ! boundary conditions ---------------------------------------- -------------------------------- SUBROUTINE order4_assemble(this) ================================ type(order4_obj) :: this ! In order to exercize iterators, the matrix ! is assembled using the general, multi-dimensional ! iterator based method, although for a 1-d ! equation it might be more efficient to use do-loops. -------------------------------- ----------------------------- SUBROUTINE order4_solve(this) ============================= type(order4_obj) :: this ----------------------------- ------------------------------------------ SUBROUTINE order4_get_solution(this, x, y) ========================================== type(order4_obj) :: this real(r8), intent(in) :: x(:) ! abscissae complex(r8), intent(out) :: y(:) ! ordinates ------------------------------------------ ----------------------------------------- SUBROUTINE order4_plot(this, x, y, label) ========================================= type(order4_obj) :: this real(r8), intent(in) :: x(:) ! abscissae complex(r8), intent(in) :: y(:) ! ordinates character(*), intent(in) :: label ----------------------------------------- ---------------------------- SUBROUTINE order4_test(this) ============================ ! Test module type(order4_obj) :: this ---------------------------- FILE ./operators.f90 ---------------- MODULE operators ================ ! Implements various operators implicit none integer, parameter, private :: r8 = selected_real_kind(12,100) complex(r8), parameter, private :: imag = (0._r8, 1._r8) character(*), parameter, private :: version = '$Id: doc.txt,v 1.1 2003/12/18 21:25:49 pletzer Exp $' contains SUBROUTINE second_derivative(x, ak, z, sigma, res, ier) SUBROUTINE derivative(order, x, ak, z, sigma, res, ier) SUBROUTINE curl_curl(ia, ib, x, ak, z, sigma, res, ier) == in module operators == ------------------------------------------------------- SUBROUTINE second_derivative(x, ak, z, sigma, res, ier) ======================================================= ! Returns the derivative of order 'order' / Gabor real(r8), intent(in) :: x, ak, z, sigma complex(r8), intent(out) :: res integer, intent(out) :: ier ------------------------------------------------------- ------------------------------------------------------- SUBROUTINE derivative(order, x, ak, z, sigma, res, ier) ======================================================= integer, intent(in) :: order real(r8), intent(in) :: x, ak, z, sigma complex(r8), intent(out) :: res integer, intent(out) :: ier ------------------------------------------------------- ------------------------------------------------------- SUBROUTINE curl_curl(ia, ib, x, ak, z, sigma, res, ier) ======================================================= integer, intent(in) :: ia, ib ! field components real(r8), intent(in) :: x(:) ! collocation position real(r8), intent(in) :: ak(:) ! Gabor wave-vector real(r8), intent(in) :: z(:) ! Gabor envelop position real(r8), intent(in) :: sigma(:) ! inverse Gabor extents complex(r8), intent(out) :: res ! result integer, intent(out) :: ier complex(r8) :: i_ktilde(size(x)) ------------------------------------------------------- FILE ./driver.f90 ------------- MODULE driver ============= == in module driver == FILE ./supralu_mod.f90 ------------------ MODULE supralu_mod ================== ! Column-compressed-storage based sparse matrix container implicit none integer, parameter, private :: r8 = selected_real_kind(12,100) character(*), parameter, private :: version = '$Id: doc.txt,v 1.1 2003/12/18 21:25:49 pletzer Exp $' type sparseC16_obj integer, pointer :: irow(:) integer, pointer :: jcol_ptr(:) complex(r8), pointer :: values(:) integer :: size, size_ptr, nnz, col_counter integer :: base ! SuperLU stuff ! column permutation specification ! permc_spec = 0: natural ordering ! permc_spec = 1: minimum degree on structure of A'*A ! permc_spec = 2: minimum degree on structure of A'+A ! permc_spec = 3: approximate minimum degree for unsymmetric matrices integer :: permc_spec end type sparseC16_obj type sparseR8_obj integer, pointer :: irow(:) integer, pointer :: jcol_ptr(:) real(r8), pointer :: values(:) integer :: size, size_ptr, nnz, col_counter integer :: base ! SuperLU stuff ! column permutation specification ! permc_spec = 0: natural ordering ! permc_spec = 1: minimum degree on structure of A'*A ! permc_spec = 2: minimum degree on structure of A'+A ! permc_spec = 3: approximate minimum degree for unsymmetric matrices integer :: permc_spec end type sparseR8_obj contains SUBROUTINE sparseC16_init(this, rank, size, base, ier) SUBROUTINE sparseC16_free(this, ier) SUBROUTINE sparseC16_copy(this, that, ier) SUBROUTINE sparseC16_compare_shapes(this, that, ier) SUBROUTINE sparseC16_set_next(this, row, newcol_flag, value, ier) SUBROUTINE sparseC16_get_size(this, nrows, ncols, ier) SUBROUTINE sparseC16_solve(this, rhs, tol, ier) SUBROUTINE sparseC16_det(this, res_mantissa, res_exponent, ier) SUBROUTINE sparseC16_eigen(Amat, Bmat, x, lambda, tol, nmax, ier) SUBROUTINE sparseC16_A_dot_x(this, x, res, ier) SUBROUTINE sparseC16_save(this, filename, ier) SUBROUTINE sparseC16_load(this, filename, ier) SUBROUTINE sparseC16_test(this, ier) SUBROUTINE sparseC16_error(this, nunit, ier) SUBROUTINE sparseR8_init(this, rank, size, base, ier) SUBROUTINE sparseR8_free(this, ier) SUBROUTINE sparseR8_copy(this, that, ier) SUBROUTINE sparseR8_compare_shapes(this, that, ier) SUBROUTINE sparseR8_set_next(this, row, newcol_flag, value, ier) SUBROUTINE sparseR8_get_size(this, nrows, ncols, ier) SUBROUTINE sparseR8_solve(this, rhs, tol, ier) SUBROUTINE sparseR8_det(this, res_mantissa, res_exponent, ier) SUBROUTINE sparseR8_A_dot_x(this, x, res, ier) SUBROUTINE sparseR8_save(this, filename, ier) SUBROUTINE sparseR8_load(this, filename, ier) SUBROUTINE sparseR8_test(this, ier) SUBROUTINE sparseR8_error(this, nunit, ier) == in module supralu_mod == ------------------------------------------------------ SUBROUTINE sparseC16_init(this, rank, size, base, ier) ====================================================== ! Constructor. ! ! rank: matrix rank (or approximation if not known) ! size: >= no of non-zero values to store. If ! size turns out to be too tight, fresh arrays will be ! reallocated and the data will be copied, however this may ! can cause performance degradation. ! ! base: indexing base (0 like C or 1 like Fortran) type(sparseC16_obj) :: this integer, intent(in) :: rank integer, intent(in) :: size integer, intent(in) :: base integer, intent(out) :: ier ------------------------------------------------------ ------------------------------------ SUBROUTINE sparseC16_free(this, ier) ==================================== ! Destructor type(sparseC16_obj) :: this integer, intent(out) :: ier ------------------------------------ ------------------------------------------ SUBROUTINE sparseC16_copy(this, that, ier) ========================================== ! Copy this into that. that should be unitialized. type(sparseC16_obj) :: this, that integer, intent(out) :: ier ------------------------------------------ ---------------------------------------------------- SUBROUTINE sparseC16_compare_shapes(this, that, ier) ==================================================== ! Compare two sparse matrices type(sparseC16_obj) :: this, that integer, intent(out) :: ier ---------------------------------------------------- ----------------------------------------------------------------- SUBROUTINE sparseC16_set_next(this, row, newcol_flag, value, ier) ================================================================= ! Set next value. Must iterate down each column ! running from left to right. ! type(sparseC16_obj) :: this integer, intent(in) :: row logical, intent(in) :: newcol_flag ! .TRUE. if new column, .FALSE. otherwise complex(r8), intent(in) :: value integer, intent(out) :: ier ----------------------------------------------------------------- ------------------------------------------------------ SUBROUTINE sparseC16_get_size(this, nrows, ncols, ier) ====================================================== ! Return no of rows and columns type(sparseC16_obj) :: this integer, intent(out) :: nrows, ncols integer, intent(out) :: ier ------------------------------------------------------ ----------------------------------------------- SUBROUTINE sparseC16_solve(this, rhs, tol, ier) =============================================== ! Solve matrix system type(sparseC16_obj) :: this complex(r8), intent(inout) :: rhs(:) ! right-hand side vector, solution vector upon return real(r8), intent(in) :: tol ! tolerance integer, intent(out) :: ier ----------------------------------------------- --------------------------------------------------------------- SUBROUTINE sparseC16_det(this, res_mantissa, res_exponent, ier) =============================================================== ! Return the determinant of A in the form of ! ! res_mantissa * 2**res_exponent type(sparseC16_obj) :: this complex(r8), intent(out) :: res_mantissa integer, intent(out) :: res_exponent integer, intent(out) :: ier --------------------------------------------------------------- ----------------------------------------------------------------- SUBROUTINE sparseC16_eigen(Amat, Bmat, x, lambda, tol, nmax, ier) ================================================================= ! Solve A * x = lambda B * x by inverse iterations type(sparseC16_obj) :: Amat, Bmat ! Must have the same sparsity pattern! complex(r8), intent(inout) :: lambda ! initial guess/computed eigenvalue complex(r8), intent(inout) :: x(:) ! initial egenvector guess/computed eigenvector real(r8), intent(inout) :: tol ! tolerance integer, intent(inout) :: nmax ! max number of iterations integer, intent(out) :: ier ----------------------------------------------------------------- ----------------------------------------------- SUBROUTINE sparseC16_A_dot_x(this, x, res, ier) =============================================== ! matrix . x => res type(sparseC16_obj) :: this complex(r8), intent(in) :: x(:) ! vector complex(r8), intent(out) :: res(:) ! result integer, intent(out) :: ier ----------------------------------------------- ---------------------------------------------- SUBROUTINE sparseC16_save(this, filename, ier) ============================================== USE ezcdf ! Save state in netCDF file 'filename' type(sparseC16_obj) :: this character(*), intent(in) :: filename integer, intent(out) :: ier ---------------------------------------------- ---------------------------------------------- SUBROUTINE sparseC16_load(this, filename, ier) ============================================== USE ezcdf ! Load state from netCDF file 'filename' type(sparseC16_obj) :: this character(*), intent(in) :: filename integer, intent(out) :: ier ---------------------------------------------- ------------------------------------ SUBROUTINE sparseC16_test(this, ier) ==================================== ! Test unit type(sparseC16_obj) :: this, that, amat, bmat integer, intent(out) :: ier ------------------------------------ -------------------------------------------- SUBROUTINE sparseC16_error(this, nunit, ier) ============================================ ! Error handler type(sparseC16_obj) :: this integer, intent(in) :: nunit ! output channel integer, intent(in) :: ier -------------------------------------------- ----------------------------------------------------- SUBROUTINE sparseR8_init(this, rank, size, base, ier) ===================================================== ! Constructor. ! ! rank: matrix rank (or approximation if not known) ! size: >= no of non-zero values to store. If ! size turns out to be too tight, fresh arrays will be ! reallocated and the data will be copied, however this may ! can cause performance degradation. ! ! base: indexing base (0 like C or 1 like Fortran) type(sparseR8_obj) :: this integer, intent(in) :: rank integer, intent(in) :: size integer, intent(in) :: base integer, intent(out) :: ier ----------------------------------------------------- ----------------------------------- SUBROUTINE sparseR8_free(this, ier) =================================== ! Destructor type(sparseR8_obj) :: this integer, intent(out) :: ier ----------------------------------- ----------------------------------------- SUBROUTINE sparseR8_copy(this, that, ier) ========================================= ! Copy this into that. that should be unitialized. type(sparseR8_obj) :: this, that integer, intent(out) :: ier ----------------------------------------- --------------------------------------------------- SUBROUTINE sparseR8_compare_shapes(this, that, ier) =================================================== ! Compare two sparse matrices type(sparseR8_obj) :: this, that integer, intent(out) :: ier --------------------------------------------------- ---------------------------------------------------------------- SUBROUTINE sparseR8_set_next(this, row, newcol_flag, value, ier) ================================================================ ! Set next value. Must iterate down each column ! running from left to right. ! type(sparseR8_obj) :: this integer, intent(in) :: row logical, intent(in) :: newcol_flag ! .TRUE. if new column, .FALSE. otherwise real(r8), intent(in) :: value integer, intent(out) :: ier ---------------------------------------------------------------- ----------------------------------------------------- SUBROUTINE sparseR8_get_size(this, nrows, ncols, ier) ===================================================== ! Return no of rows and columns type(sparseR8_obj) :: this integer, intent(out) :: nrows, ncols integer, intent(out) :: ier ----------------------------------------------------- ---------------------------------------------- SUBROUTINE sparseR8_solve(this, rhs, tol, ier) ============================================== ! Solve matrix system type(sparseR8_obj) :: this real(r8), intent(inout) :: rhs(:) ! right-hand side vector, solution vector upon return real(r8), intent(in) :: tol ! tolerance integer, intent(out) :: ier ---------------------------------------------- -------------------------------------------------------------- SUBROUTINE sparseR8_det(this, res_mantissa, res_exponent, ier) ============================================================== ! Return the determinant of A in the form of ! ! res_mantissa * 2**res_exponent type(sparseR8_obj) :: this real(r8), intent(out) :: res_mantissa integer, intent(out) :: res_exponent integer, intent(out) :: ier -------------------------------------------------------------- ---------------------------------------------- SUBROUTINE sparseR8_A_dot_x(this, x, res, ier) ============================================== ! matrix . x => res type(sparseR8_obj) :: this real(r8), intent(in) :: x(:) ! vector real(r8), intent(out) :: res(:) ! result integer, intent(out) :: ier ---------------------------------------------- --------------------------------------------- SUBROUTINE sparseR8_save(this, filename, ier) ============================================= USE ezcdf ! Save state in netCDF file 'filename' type(sparseR8_obj) :: this character(*), intent(in) :: filename integer, intent(out) :: ier --------------------------------------------- --------------------------------------------- SUBROUTINE sparseR8_load(this, filename, ier) ============================================= USE ezcdf ! Load state from netCDF file 'filename' type(sparseR8_obj) :: this character(*), intent(in) :: filename integer, intent(out) :: ier --------------------------------------------- ----------------------------------- SUBROUTINE sparseR8_test(this, ier) =================================== ! Test unit type(sparseR8_obj) :: this, that, amat, bmat integer, intent(out) :: ier ----------------------------------- ------------------------------------------- SUBROUTINE sparseR8_error(this, nunit, ier) =========================================== ! Error handler type(sparseR8_obj) :: this integer, intent(in) :: nunit ! output channel integer, intent(in) :: ier ------------------------------------------- FILE ./plotmtv.f90 -------------- MODULE plotmtv ============== ! Interface to plotmtv implicit none integer, parameter, private :: r8 = selected_real_kind(12,100) character(*), parameter, private :: version = '$Id: doc.txt,v 1.1 2003/12/18 21:25:49 pletzer Exp $' type plotmtv_obj integer :: plotid ! id number logical :: hold ! set this to true for multiplots logical :: save ! set this to true if file needs to be saved logical :: freshstart ! true if new file is to be created logical :: display ! true to dosplay plot character(80) :: toplabel, xlabel, ylabel, zlabel, linelabel integer :: linewidth ! 1-32 integer :: markersize ! 1-30 integer :: fontsize ! 1-30 end type plotmtv_obj contains SUBROUTINE plotmtv_init(this) SUBROUTINE plotmtv_free(this) SUBROUTINE plotmtv_plot(this, x, y, line, marker) SUBROUTINE plotmtv_test(this) SUBROUTINE plotmtv_set_line(line, code) SUBROUTINE plotmtv_set_color(color, code) SUBROUTINE plotmtv_set_marker(marker, code) == in module plotmtv == ----------------------------- SUBROUTINE plotmtv_init(this) ============================= type(plotmtv_obj) :: this ----------------------------- ----------------------------- SUBROUTINE plotmtv_free(this) ============================= type(plotmtv_obj) :: this ----------------------------- ------------------------------------------------- SUBROUTINE plotmtv_plot(this, x, y, line, marker) ================================================= ! Dump x and y into file tableN.mtv where N is the plotid. Fire ! plotmtv if display = true. Multiple line plots are obtained ! by calling plot in sequence and setting hold = true. type(plotmtv_obj) :: this real(r8), intent(in) :: x(:), y(:) character(*) :: line, marker ------------------------------------------------- ----------------------------- SUBROUTINE plotmtv_test(this) ============================= type(plotmtv_obj) :: this ----------------------------- --------------------------------------- SUBROUTINE plotmtv_set_line(line, code) ======================================= character, intent(in) :: line integer, intent(out) :: code --------------------------------------- ----------------------------------------- SUBROUTINE plotmtv_set_color(color, code) ========================================= character, intent(in) :: color integer, intent(out) :: code ----------------------------------------- ------------------------------------------- SUBROUTINE plotmtv_set_marker(marker, code) =========================================== character, intent(in) :: marker integer, intent(out) :: code ------------------------------------------- FILE ./singular.f90 --------------- MODULE singular =============== ! Module for Singular test problem use supralu_mod use indexer use Gabor use kzgrid use xgrid use operators use plotmtv implicit none integer, parameter, private :: r8 = selected_real_kind(12,100) character(*), parameter, private :: version = '$Id: doc.txt,v 1.1 2003/12/18 21:25:49 pletzer Exp $' integer, parameter :: Nd = 1 ! space dimension integer, parameter :: Nc = 1 ! no of components type singular_obj integer :: Nz ! No of Gabor intervals integer :: Nj ! Fourier modes -Nj .. + Nj real(r8) :: sampling ! sampling rate > 1 real(r8) :: alpha ! eqn parameter real(r8) :: tol ! sparsity tolerance complex(r8) :: bc0, bc1 ! Boundary conditions real(r8) :: xmin, xmax ! domain boundary real(r8) :: eps real(r8) :: sigma integer :: Nx type(sparseC16_obj) :: amat type(indexer_obj) :: iter2 type(kzgrid_obj) :: kzg complex(r8), dimension(:), pointer :: xvec, bvec end type singular_obj contains SUBROUTINE singular_init(this, Nz, Nj, sigma_dx2, sampling) SUBROUTINE singular_free(this) SUBROUTINE singular_set_domain(this, xmin, xmax) SUBROUTINE singular_set_tol(this, tol) SUBROUTINE singular_set_bcs(this, bc0, bc1) SUBROUTINE singular_assemble(this) SUBROUTINE singular_solve(this) SUBROUTINE singular_get_solution(this, x, y) SUBROUTINE singular_plot(this, x, y) SUBROUTINE singular_test(this) == in module singular == ----------------------------------------------------------- SUBROUTINE singular_init(this, Nz, Nj, sigma_dx2, sampling) =========================================================== type(singular_obj) :: this integer, intent(in) :: Nz ! No of Gabor intervals integer, intent(in) :: Nj ! Fourier modes -Nj .. + Nj real(r8), intent(in) :: sigma_dx2 ! inverse Gabor extent normalized to grid size real(r8), intent(in) :: sampling ! sampling rate > 1 ----------------------------------------------------------- ------------------------------ SUBROUTINE singular_free(this) ============================== type(singular_obj) :: this ------------------------------ ------------------------------------------------ SUBROUTINE singular_set_domain(this, xmin, xmax) ================================================ type(singular_obj) :: this real(r8), intent(in) :: xmin, xmax ! domain boundary ------------------------------------------------ -------------------------------------- SUBROUTINE singular_set_tol(this, tol) ====================================== type(singular_obj) :: this real(r8), intent(in) :: tol ! sparsity tolerance ! A low value of tol will make the matrices less sparse -------------------------------------- ------------------------------------------- SUBROUTINE singular_set_bcs(this, bc0, bc1) =========================================== type(singular_obj) :: this complex(r8), intent(in) :: bc0, bc1 ! Dirichlet boundary conditions ------------------------------------------- ---------------------------------- SUBROUTINE singular_assemble(this) ================================== type(singular_obj) :: this ! In order to exercize iterators, the matrix ! is assembled using the general, multi-dimensional ! iterator based method, although for a 1-d ! equation it might be more efficient to use do-loops. ---------------------------------- ------------------------------- SUBROUTINE singular_solve(this) =============================== type(singular_obj) :: this ------------------------------- -------------------------------------------- SUBROUTINE singular_get_solution(this, x, y) ============================================ type(singular_obj) :: this real(r8), intent(in) :: x(:) ! abscissae complex(r8), intent(out) :: y(:) ! ordinates -------------------------------------------- ------------------------------------ SUBROUTINE singular_plot(this, x, y) ==================================== type(singular_obj) :: this real(r8), intent(in) :: x(:) ! abscissae complex(r8), intent(in) :: y(:) ! ordinates ------------------------------------ ------------------------------ SUBROUTINE singular_test(this) ============================== ! Test module type(singular_obj) :: this ------------------------------ FILE ./eigenwave.f90 ---------------- MODULE eigenwave ================ ! Module for Eigenvalue wave test problem use supralu_mod use indexer use Gabor use kzgrid use xgrid use operators use plotmtv implicit none integer, parameter, private :: r8 = selected_real_kind(12,100) real(r8), parameter, private :: twopi = 6.2831853071795865_r8 character(*), parameter, private :: version = '$Id: doc.txt,v 1.1 2003/12/18 21:25:49 pletzer Exp $' integer, parameter :: Nd = 1 ! space dimension integer, parameter :: Nc = 1 ! no of components type eigenwave_obj integer :: Nz ! No of Gabor intervals integer :: Nj ! Fourier modes -Nj .. + Nj real(r8) :: sampling ! sampling rate > 1 real(r8) :: tol ! sparsity tolerance real(r8) :: xmin, xmax ! domain boundary real(r8) :: eps real(r8) :: sigma integer :: Nx type(sparseC16_obj) :: amat, bmat type(indexer_obj) :: iter2 type(kzgrid_obj) :: kzg end type eigenwave_obj contains SUBROUTINE eigenwave_init(this, Nz, Nj, sigma_dx2, sampling) SUBROUTINE eigenwave_free(this) SUBROUTINE eigenwave_set_domain(this, xmin, xmax) SUBROUTINE eigenwave_set_tol(this, tol) SUBROUTINE eigenwave_assemble(this) SUBROUTINE eigenwave_solve(this, xvec, lambda, nmax) SUBROUTINE eigenwave_get_solution(this, xvec, x, y) SUBROUTINE eigenwave_plot(this, x, y, title) SUBROUTINE eigenwave_test(this) == in module eigenwave == ------------------------------------------------------------ SUBROUTINE eigenwave_init(this, Nz, Nj, sigma_dx2, sampling) ============================================================ type(eigenwave_obj) :: this integer, intent(in) :: Nz ! No of Gabor intervals integer, intent(in) :: Nj ! Fourier modes -Nj .. + Nj real(r8), intent(in) :: sigma_dx2 ! inverse Gabor extent normalized to grid size real(r8), intent(in) :: sampling ! sampling rate > 1 ------------------------------------------------------------ ------------------------------- SUBROUTINE eigenwave_free(this) =============================== type(eigenwave_obj) :: this ------------------------------- ------------------------------------------------- SUBROUTINE eigenwave_set_domain(this, xmin, xmax) ================================================= type(eigenwave_obj) :: this real(r8), intent(in) :: xmin, xmax ! domain boundary ------------------------------------------------- --------------------------------------- SUBROUTINE eigenwave_set_tol(this, tol) ======================================= type(eigenwave_obj) :: this real(r8), intent(in) :: tol ! sparsity tolerance ! A low value of tol will make the matrices less sparse --------------------------------------- ----------------------------------- SUBROUTINE eigenwave_assemble(this) =================================== type(eigenwave_obj) :: this ! In order to exercize iterators, the matrix ! is assembled using the general, multi-dimensional ! iterator based method, although for a 1-d ! equation it might be more efficient to use do-loops. ----------------------------------- ---------------------------------------------------- SUBROUTINE eigenwave_solve(this, xvec, lambda, nmax) ==================================================== type(eigenwave_obj) :: this complex(r8), intent(inout) :: xvec(:) ! initial egenvector guess on input/eigenvector on output complex(r8), intent(inout) :: lambda ! initial egenvalue guess on input/eigenvalue on output integer, intent(inout) :: nmax ! max number of iterations on input/no of iterations on output ---------------------------------------------------- --------------------------------------------------- SUBROUTINE eigenwave_get_solution(this, xvec, x, y) =================================================== type(eigenwave_obj) :: this complex(r8), intent(in) :: xvec(:) ! eigenvector returned by solve method real(r8), intent(in) :: x(:) ! abscissae complex(r8), intent(out) :: y(:) ! ordinates --------------------------------------------------- -------------------------------------------- SUBROUTINE eigenwave_plot(this, x, y, title) ============================================ type(eigenwave_obj) :: this real(r8), intent(in) :: x(:) ! abscissae complex(r8), intent(in) :: y(:) ! ordinates character(*), intent(in) :: title -------------------------------------------- ------------------------------- SUBROUTINE eigenwave_test(this) =============================== USE supralu_mod ! Test module type(eigenwave_obj) :: this ------------------------------- FILE ./eigenbessel.f90 ------------------ MODULE eigenbessel ================== ! Module for Eigenvalue Bessel test problem use supralu_mod use indexer use Gabor use kzgrid use xgrid use operators use plotmtv implicit none integer, parameter, private :: r8 = selected_real_kind(12,100) character(*), parameter, private :: version = '$Id: doc.txt,v 1.1 2003/12/18 21:25:49 pletzer Exp $' integer, parameter :: Nd = 1 ! space dimension integer, parameter :: Nc = 1 ! no of components type eigenbessel_obj integer :: Nz ! No of Gabor intervals integer :: Nj ! Fourier modes -Nj .. + Nj real(r8) :: sampling ! sampling rate > 1 real(r8) :: tol ! sparsity tolerance real(r8) :: xmin, xmax ! domain boundary ! equation parameters integer :: mpol real(r8) :: eps real(r8) :: sigma integer :: Nx type(sparseC16_obj) :: amat, bmat type(indexer_obj) :: iter2 type(kzgrid_obj) :: kzg end type eigenbessel_obj contains SUBROUTINE eigenbessel_init(this, Nz, Nj, sigma_dx2, sampling) SUBROUTINE eigenbessel_free(this) SUBROUTINE eigenbessel_set_domain(this, xmin, xmax) SUBROUTINE eigenbessel_set_tol(this, tol) SUBROUTINE eigenbessel_assemble(this) SUBROUTINE eigenbessel_solve(this, xvec, lambda, nmax) SUBROUTINE eigenbessel_get_solution(this, xvec, x, y) SUBROUTINE eigenbessel_plot(this, x, y, title) SUBROUTINE eigenbessel_test(this) == in module eigenbessel == -------------------------------------------------------------- SUBROUTINE eigenbessel_init(this, Nz, Nj, sigma_dx2, sampling) ============================================================== type(eigenbessel_obj) :: this integer, intent(in) :: Nz ! No of Gabor intervals integer, intent(in) :: Nj ! Fourier modes -Nj .. + Nj real(r8), intent(in) :: sigma_dx2 ! inverse Gabor extent normalized to grid size real(r8), intent(in) :: sampling ! sampling rate > 1 -------------------------------------------------------------- --------------------------------- SUBROUTINE eigenbessel_free(this) ================================= type(eigenbessel_obj) :: this --------------------------------- --------------------------------------------------- SUBROUTINE eigenbessel_set_domain(this, xmin, xmax) =================================================== type(eigenbessel_obj) :: this real(r8), intent(in) :: xmin, xmax ! domain boundary --------------------------------------------------- ----------------------------------------- SUBROUTINE eigenbessel_set_tol(this, tol) ========================================= type(eigenbessel_obj) :: this real(r8), intent(in) :: tol ! sparsity tolerance ! A low value of tol will make the matrices less sparse ----------------------------------------- ------------------------------------- SUBROUTINE eigenbessel_assemble(this) ===================================== type(eigenbessel_obj) :: this ! In order to exercize iterators, the matrix ! is assembled using the general, multi-dimensional ! iterator based method, although for a 1-d ! equation it might be more efficient to use do-loops. ------------------------------------- ------------------------------------------------------ SUBROUTINE eigenbessel_solve(this, xvec, lambda, nmax) ====================================================== type(eigenbessel_obj) :: this complex(r8), intent(inout) :: xvec(:) ! initial egenvector guess on input/eigenvector on output complex(r8), intent(inout) :: lambda ! initial egenvalue guess on input/eigenvalue on output integer, intent(inout) :: nmax ! max number of iterations on input/no of iterations on output ------------------------------------------------------ ----------------------------------------------------- SUBROUTINE eigenbessel_get_solution(this, xvec, x, y) ===================================================== type(eigenbessel_obj) :: this complex(r8), intent(in) :: xvec(:) ! eigenvector returned by solve method real(r8), intent(in) :: x(:) ! abscissae complex(r8), intent(out) :: y(:) ! ordinates ----------------------------------------------------- ---------------------------------------------- SUBROUTINE eigenbessel_plot(this, x, y, title) ============================================== type(eigenbessel_obj) :: this real(r8), intent(in) :: x(:) ! abscissae complex(r8), intent(in) :: y(:) ! ordinates character(*), intent(in) :: title ---------------------------------------------- --------------------------------- SUBROUTINE eigenbessel_test(this) ================================= USE supralu_mod ! Test module type(eigenbessel_obj) :: this ---------------------------------