lusol_test.f90 Source File

!***************************************************************************************************
!>
!  Main program for EZ test.

    program main

    use lusol_precision, wp => rp
    use lusol_ez_module,     only: solve

    implicit none

    call test_1()
    call test_2()

    contains

    subroutine test_1()

        ! define a 3x3 dense system to solve:
        integer,parameter :: m = 3 !! number of rows in `A` matrix
        integer,parameter :: n = 3 !! number of columns in `A` matrix
        real(wp),dimension(m),parameter :: b = real([1,2,3],wp)
        integer,dimension(m*n),parameter :: icol = [1,1,1,2,2,2,3,3,3]
        integer,dimension(m*n),parameter :: irow = [1,2,3,1,2,3,1,2,3]
        real(wp),dimension(m*n),parameter :: a = real([1,4,7,2,5,88,3,66,9],wp)
        real(wp),dimension(m,1),parameter :: b_vec = reshape(b,[m,1])
        real(wp),dimension(m,n),parameter :: a_mat = reshape(a,[3,3])

        real(wp),dimension(n) :: x
        real(wp),dimension(n,1) :: x_vec
        integer :: istat

        write(*,*) ''
        write(*,*) '***** test 1'
        write(*,*) ''

        call solve(n,m,m*n,irow,icol,a,b,x,istat)

        x_vec(1:3,1) = x

        write(*,*) ''
        write(*,*) 'istop = ', istat
        write(*,*) ''
        write(*,'(1P,A,*(E16.6))') 'x       = ', x
        write(*,'(1P,A,*(E16.6))') 'A*x     = ', matmul(a_mat, x_vec)
        write(*,'(1P,A,*(E16.6))') 'A*x - b = ', matmul(a_mat, x_vec) - b_vec

        if (any(abs(matmul(a_mat, x_vec) - b_vec) > 1.0e-12)) error stop 'TEST FAILED'

    end subroutine test_1

    subroutine test_2()

        ! another test case (with n>m). compare to scipy.

        !!```
        !! >>> a
        !! array([[  4. ,   5. ,  66. ,   0.1],
        !!        [  1. ,  -3. ,   8. ,  -9. ],
        !!        [ 11. ,   3. , -87. ,   2. ]])
        !! >>> b
        !! array([1, 2, 3])
        !! >>> scipy.sparse.linalg.lsqr(a, b)
        !! (array([ 0.26437473,  0.04901579, -0.00426183, -0.21297414]), 1, 3, 5.5785963493386424e-12, 5.5785963493386424e-12, 110.70234866523838, 15.316189089999897, 6.119548932941366e-10, 0.343034538979173, array([0., 0., 0., 0.]))
        !! >>>
        !!```

        integer,parameter :: m = 3 !! number of rows in `A` matrix
        integer,parameter :: n = 4 !! number of columns in `A` matrix
        real(wp),dimension(m),parameter :: b = real([1,2,3],wp)
        integer,dimension(m*n),parameter :: icol = [1,1,1,2,2,2,3,3,3,4,4,4]
        integer,dimension(m*n),parameter :: irow = [1,2,3,1,2,3,1,2,3,1,2,3]
        real(wp),dimension(m*n),parameter :: a = [4.1_wp,1.1_wp,11.1_wp,&
                                                  5.1_wp,-3.1_wp,3.1_wp,&
                                                  66.1_wp,8.1_wp,-87.1_wp,&
                                                  0.1_wp,-9.1_wp,2.1_wp]
        real(wp),dimension(m,1),parameter :: b_vec = reshape(b,[m,1])
        real(wp),dimension(m,n),parameter :: a_mat = reshape(a,[3,4])

        real(wp),dimension(n) :: x
        real(wp),dimension(n,1) :: x_vec
        integer :: istat

        write(*,*) ''
        write(*,*) '***** test 2'
        write(*,*) ''

        call solve(n,m,m*n,irow,icol,a,b,x,istat)

        x_vec(1:4,1) = x

        write(*,*) ''
        write(*,*) 'istop = ', istat
        write(*,*) ''
        write(*,'(1P,A,*(E16.6))') 'x       = ', x
        write(*,'(1P,A,*(E16.6))') 'A*x     = ', matmul(a_mat, x_vec)
        write(*,'(1P,A,*(E16.6))') 'A*x - b = ', matmul(a_mat, x_vec) - b_vec

        if (any(abs(matmul(a_mat, x_vec) - b_vec) > 1.0e-12)) error stop 'TEST FAILED'

    end subroutine test_2


    end program main
!***************************************************************************************************