sparsetest.f90

module auxsubs

contains

  subroutine vecinit(n, b)
    ! initialize source vector b
    ! source: simple quadrupol
    
    use kinds
    implicit none    

    real(kind=REAL8), dimension(:,:), pointer :: b
    integer                                   :: n
    real(kind=REAL8)                          :: rho, dist, offset
    integer                                   :: i1, i2

    allocate(b(n*n,1))
    
    rho = n*n                                ! total charge is constant

    dist = 0.4
    offset = (1.0 - dist)/2.0
    i1 = floor(offset*n) + 1
    i2 = floor((1.0 - offset)*n)

    b(n*i1 + i1, 1) =  rho                      ! b(i1, i1)
    b(n*i2 + i2, 1) =  rho                      ! b(i2, i2)
    b(n*i1 + i2, 1) = -rho                      ! b(i1, i2)
    b(n*i2 + i1, 1) = -rho                      ! b(i2, i1)

    return
  end subroutine vecinit

  subroutine matinit(n, cols, rows, values)
    ! initialize sparse matrix for poisson equation
    
    use kinds
    implicit none    
    real(kind=REAL8), dimension(:), pointer   :: values
    integer, dimension(:), pointer            :: cols, rows
    integer                                   :: n
    integer                                   :: nz, i, offset

    nz = 3*n*n - n - 1

    allocate(cols(nz))
    allocate(rows(nz))
    allocate(values(nz))
    
    offset = 0

    ! diagonal
    do i=1, n*n
       rows(offset + i) = i
       cols(offset + i) = i
       values(offset + i) = -4.0D0
    end do
    offset = offset + n*n

    ! lower subdiagonal    
    do i=2, n*n
       rows(offset + i -1) = i
       cols(offset + i -1) = i-1
       values(offset + i -1) = 1.0D0
    end do
    offset = offset + n*n - 1
    
    ! lower sub-n-diagonal    
    do i=n+1, n*n
       rows(offset + i - n) = i
       cols(offset + i - n) = i - n
       values(offset + i - n) = 1.0D0
    end do
    offset = offset + n*n - n
    
    return
  end subroutine matinit

  subroutine solve(cols, rows, values, b)
    ! solve sparse system

    use kinds
    implicit none    

    integer, dimension(:), pointer             :: rows, cols
    real(kind=REAL8), dimension(:), pointer    :: values
    real(kind=REAL8), dimension(:,:), pointer  :: b
    integer                                    :: n, nz, info, i, j, k
    real(kind=REAL8)                           :: cond, value
    real(kind=REAL8), dimension(150)           :: global
    integer, dimension(3)                      :: inertia

    
    n  = size(b, 1)
    nz = size(rows, 1)

    ! init sparse lib
    call dslein(n, 1, 6, global, info)
    
    ! enter matrix structure
    do k=1, nz
       i = rows(k)
       j = cols(k)
       call dslei1(i, j, global, info)
    end do
    call dsleif(global, info)
    
    ! reorder and factorize symbolically
    call dsleor(10, global, info)
    
    ! enter matrix values
    do k=1, nz
       i = rows(k)
       j = cols(k)
       value = values(k)
       call dslev1(i, j, value, global, info)
    end do
    
    ! factorize
    call dslefa(0.1D0, inertia, global, info)
    
    ! solve for given rhs
    call dslesl(1, b, n, global, info)
    
    return
  end subroutine solve


  subroutine show_result(b)
    ! display the values
    
    use kinds
    implicit none    
    real(kind=REAL8), dimension(:,:), pointer :: b
    integer                                   :: n, i

    n  = size(b, 1)
    do i=1, n
       print '(f16.3)', b(i,1)
    enddo

    return
  end subroutine show_result
end module auxsubs


program sparsetest

  ! example program solving ax=b, a sparse,  via veclib routine dslefs
  ! set number of threads with MLIB_NUMBER_OF_THREADS
  
  use kinds
  use auxsubs
  implicit none
  
  integer, dimension(:), pointer             :: rows, cols
  real(kind=REAL8), dimension(:), pointer    :: values
  real(kind=REAL8), dimension(:,:), pointer  :: b
  integer                                    :: n
  character                                  :: c

  ! get size of problem
  print*, 'enter (linear) size:'
  read*, n

  ! initialize system matrix
  call matinit(n, cols, rows, values)

  ! initialize source vector
  call vecinit(n, b)

  ! solve the equation
  call solve(cols, rows, values, b)

  ! display the result
  print*, 'dump result (y/n)?'
  read*, c
  if ((c == 'y') .or. (c == 'Y')) then
     call show_result(b)
  endif
  
end program sparsetest