


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
---------------------------------
