test_2 Subroutine

subroutine test_2() 

 >>> 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.]))
 >>>

Arguments

None

Calls

proc~~test_2~~CallsGraph proc~test_2 main::test_2 proc~solve lusol_ez_module::solve proc~test_2->proc~solve proc~lu1fac lusol::lu1fac proc~solve->proc~lu1fac proc~lu6sol lusol::lu6sol proc~solve->proc~lu6sol proc~lu1fad lusol::lu1fad proc~lu1fac->proc~lu1fad proc~lu1or1 lusol::lu1or1 proc~lu1fac->proc~lu1or1 proc~lu1or2 lusol::lu1or2 proc~lu1fac->proc~lu1or2 proc~lu1or3 lusol::lu1or3 proc~lu1fac->proc~lu1or3 proc~lu1or4 lusol::lu1or4 proc~lu1fac->proc~lu1or4 proc~lu1pq1 lusol::lu1pq1 proc~lu1fac->proc~lu1pq1 proc~lu1slk lusol::lu1slk proc~lu1fac->proc~lu1slk proc~lu6chk lusol::lu6chk proc~lu1fac->proc~lu6chk proc~lu6l lusol::lu6L proc~lu6sol->proc~lu6l proc~lu6ld lusol::lu6LD proc~lu6sol->proc~lu6ld proc~lu6lt lusol::lu6Lt proc~lu6sol->proc~lu6lt proc~lu6u lusol::lu6U proc~lu6sol->proc~lu6u proc~lu6ut lusol::lu6Ut proc~lu6sol->proc~lu6ut proc~hbuild lusol::Hbuild proc~lu1fad->proc~hbuild proc~hchange lusol::Hchange proc~lu1fad->proc~hchange proc~hdelete lusol::Hdelete proc~lu1fad->proc~hdelete proc~lu1ful lusol::lu1ful proc~lu1fad->proc~lu1ful proc~lu1gau lusol::lu1gau proc~lu1fad->proc~lu1gau proc~lu1mar lusol::lu1mar proc~lu1fad->proc~lu1mar proc~lu1mrp lusol::lu1mRP proc~lu1fad->proc~lu1mrp proc~lu1msp lusol::lu1mSP proc~lu1fad->proc~lu1msp proc~lu1mxc lusol::lu1mxc proc~lu1fad->proc~lu1mxc proc~lu1mxr lusol::lu1mxr proc~lu1fad->proc~lu1mxr proc~lu1pen lusol::lu1pen proc~lu1fad->proc~lu1pen proc~lu1pq2 lusol::lu1pq2 proc~lu1fad->proc~lu1pq2 proc~lu1pq3 lusol::lu1pq3 proc~lu1fad->proc~lu1pq3 proc~lu1rec lusol::lu1rec proc~lu1fad->proc~lu1rec proc~hinsert lusol::Hinsert proc~hbuild->proc~hinsert proc~hdown lusol::Hdown proc~hchange->proc~hdown proc~hup lusol::Hup proc~hchange->proc~hup proc~hdelete->proc~hchange proc~lu1dcp lusol::lu1DCP proc~lu1ful->proc~lu1dcp proc~lu1dpp lusol::lu1DPP proc~lu1ful->proc~lu1dpp proc~hinsert->proc~hup proc~jdamax lusol::jdamax proc~lu1dcp->proc~jdamax proc~lu1dpp->proc~jdamax

Called by

proc~~test_2~~CalledByGraph proc~test_2 main::test_2 program~main~2 main program~main~2->proc~test_2

Source Code

    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