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