Solving Linear Algebra Problem

Linear equation

Linear equations try to solve following kind of equation

$$\begin{array}{}\text{(1)}& \mathbf{\text{A}}\mathbf{\text{x}}=\mathbf{\text{b}}\end{array}$$
where

$\mathbf{\text{A}}$ is a NxN matrix,

$\mathbf{\text{x}}$ is a Nx1 and

$\mathbf{\text{b}}$ is a Nx1 column. There are two ways to solve it by using linear package, LAPACK. First is using the subroutine

$dgesv$ which solve

$\mathbf{\text{x}}$ immediately. Second method is that we can find the inverse matrix of

$\mathbf{\text{A}}$ , then multiplied by

$\mathbf{\text{b}}$. Finding inverse matrix in LAPACK require two subroutine,

$dgetrf$ and

$dgetri$.

Example 1

In the first example we first solve

$\mathbf{\text{x}}$ and then verified the result by multiply

$\mathbf{\text{x}}$ by the matrix

$\mathbf{\text{A}}$. I then solve this equation for

${10}^6$ times to compare its run time to another method. Finally I computer the error of solution by calculate the difference between

$\mathbf{\text{A}}\mathbf{\text{x}}$ and

$\mathbf{\text{b}}$ of two method.

As you can see both run time and accuracy are superior when using subroutine

$dgesv$.

Example 2

However, when the matrix

$\mathbf{\text{A}}$ is a symmetric matrix the later method has a higher accuracy despite a slower speed.

Code

program main
    ! Example program that solves the matrix
    ! equation A * x = B using LAPACK.
    implicit none

    double precision :: a(5, 5)  ! Matrix A.
    double precision :: b(5)     ! Matrix B.
    ! Pivot indices (list of swap operations).
    double precision :: pivot(5) 
    integer :: rc       ! Return code.
    double precision :: aa(5, 5), bb(5), aaa(5, 5)
    double precision :: work2(5)
    double precision :: start, finish
    integer :: i

    a = reshape([4.d0 , 1.d0 , 1.d0 , 0.d0 , 1.d0,  &
                 -1.d0, 3.d0 , 1.d0 , 1.d0 , 0.d0,  &
                 2.d0 , 1.d0 , 5.d0 , -1.d0, -1.d0, &
                 -1.d0, -1.d0, -1.d0, 4.d0 , 0.d0,  &
                 0.d0 , 2.d0 , -1.d0, 1.d0 , 4.d0], &
                 [5, 5])
    b = [6.d0, 6.d0, 6.d0, 6.d0, 6.d0]

    ! method one
    call cpu_time(start)
    do i = 1, 1000000
        aa = a
        bb = b
        call dgesv(5, 1, aa, 5, pivot, bb, 5, rc)
    end do
    call cpu_time(finish)

    if (rc == 0) then
        write(*, '(a)') repeat('=', 50)
        print *, 'Using dgesv:'
        print '(a, 5(f0.4, ", "))', 'x = ', bb
    else
        print '(a, i0)', 'Error: ', rc
    end if
    print '(a, 5(f0.4, ", "))', 'Ax = ', matmul(a, bb)
    write(*, *) 'cost(run 10^6 times)', finish-start, 's'
    write(*, *) 'accuracy: ', sum(matmul(a, bb)-b)

    ! method two
    write(*, '(a)') repeat('=', 50)
    print *, 'Using dgetrf+dgetri:'

    call cpu_time(start)
    do i = 1, 1000000
        aaa = a
        call dgetrf(5, 5, aaa, 5, pivot, rc)
        call dgetri(5, aaa, 5, pivot, work2, 5, rc)
    end do
    call cpu_time(finish)

    if (rc == 0) then
        print '(a, 5(f0.4, ", "))', 'x = ', matmul(aaa, b)
        print '(a, 5(f0.4, ", "))', 'Ax = ', matmul(a, matmul(aaa, b))
    else
        print '(a, i0)', 'Error: ', rc
    end if
    write(*, *) 'cost(run 10^6 times)', finish-start, 's'
    write(*, *) 'accuracy: ', sum(matmul(a, matmul(aaa, b))-b)
    write(*, '(a)') repeat('=', 50)
end program main

Eigen-system

Using

$dgeev$ we can solve the eigenvalue and eigenvector of the matrix

$\mathbf{\text{A}}$. Following print the result of five eigenvalue and five eigenvector. By definition the relation between eigenvalue and eigenvector is given by

$$\begin{array}{}\text{(2)}& \mathbf{\text{A}}\ast v(j)=\lambda (j)\ast v(j)\end{array}$$
So I use this relation to test the accuracy of my result.

Example 1

Example 2

Same procedure for finding another eigenvalue and eigenvector of second matrix.

Wigner's semi-circle and Girko's Circular Laws in random matrix

Real-Symmetric NxN matrix

I create a symmetric matrix with the size of

$5000\times 5000$ that fill with uniform distribution in

$[-0.5,0.5]$. Then I solve the eigenvalue of this system, since this matrix is symmetric, all the eigenvalue are real. I plot the probability distribution of eigenvalue, the distribution when the size of the array become infinite is given by

$$\begin{array}{}\text{(3)}& P(\lambda )={\displaystyle \frac{2}{\pi R^2}}\sqrt{R^2-x^2}\end{array}$$

Code

program main
    ! find the eigenvalue and eigenvector of asymmetric matrix A
    ! equation A * x = B using LAPACK.
    implicit none

    integer, parameter :: n=5
    double precision :: a(n, n) 
    integer :: rc, lwork
    double precision :: aa(n, n)
    double precision :: wr(n), wi(n), work(1000), vl(n, n), vr(n, n)
    integer :: i, j
    double complex :: lambda(n), v(n, n)

    a = reshape([4.d0 , 1.d0 , 1.d0 , 0.d0 , 1.d0,  &
                 -1.d0, 3.d0 , 1.d0 , 1.d0 , 0.d0,  &
                 2.d0 , 1.d0 , 5.d0 , -1.d0, -1.d0, &
                 -1.d0, -1.d0, -1.d0, 4.d0 , 0.d0,  &
                 0.d0 , 2.d0 , -1.d0, 1.d0 , 4.d0], &
                 [5, 5])
    aa = a

    lwork = -1
    CALL DGEEV('N', 'V', n, aa, n, wr, wi, vl, n, &
                vr, n, work, lwork, rc )
    lwork = MIN( 1000, INT( work( 1 ) ) )
    CALL DGEEV('N', 'V', n, aa, n, wr, wi, vl, n, &
                vr, n, work, lwork, rc )

    ! write eigenvalue
    write(*, *) 'Eigenvalue'
    do i = 1, n
        write(*, *) i, wr(i), wi(i)
    end do
    write(*, *) repeat('=', 70)

    lambda = 0.d0
    do i = 1, n
        lambda(i) = complex(wr(i), wi(i))
    end do

    ! write eigenvector
    write(*, *) 'Eigenvector'
    do i = 1, n
        j = 1
        do while( j.le.n )
            if( wi(j).eq.0.d0) then
                write(*, 9998, advance='no') vr( i, j )
                v(i, j) = complex(vr(i, j), 0.d0)
                j = j + 1
            else
                write(*, 9999, advance='no') vr( i, j ), vr( i, j+1 )
                write(*,9999,advance='no') vr( i, j ), -vr( i, j+1 )
                v(i, j) = complex(vr(i, j), vr(i, j+1))
                v(i, j+1) = complex(vr(i, j), -vr(i, j+1))
                j = j + 2
            end if
        end do
        write(*,*)
    end do
    write(*, *) repeat('=', 70)

    write(*, *) 'check: A * v(j) - lambda(j) * v(j) ?= 0'
    do j = 1, 5
        print 9999, sum(matmul(a, v(:, j)) - lambda(j) * v(:, j))
    end do
    write(*, *) repeat('=', 70)

9998 FORMAT( 11(:,1X,F6.2) )
9999 FORMAT( 11(:,1X,'(',F6.2,',',F6.2,')') )
end program main

Girko’s circle law

If we use more general asymmetric N-by-N matrix the eigenvalue would be a complex number. I can the plot distribution of real part and imaginary part in the complex plane. I then scale the data by

$\sqrt{12/n}$, so it seats on a unit circle. As you can see by increasing the size of the matrix, the dot are filling more in the unit circle.

Linear boundary-value problem ODE

The last question are trying to solve a boundary-value ODE, using Finite-Difference method we can easily transform this ODE into a problem of linear algebra in which I can use previous result to solve it. I divide the interval into different number, and we can see that the method converge to the exact solution very fast. It only require few point to converges to the exact solution to this ODE.

program main
    ! find the eigenvalue and eigenvector of asymmetric matrix A
    ! equation A * x = B using LAPACK.
    implicit none

    integer, parameter :: n=5000
    integer, parameter :: lmax = 100000
    double precision :: a(n, n) 
    integer :: rc, lwork
    double precision :: aa(n, n)
    double precision :: wr(n), wi(n), work(lmax), vl(n, n), vr(n, n)
    integer :: i, j
    double complex :: lambda(n)
    double precision :: rand

    call random_seed()
    do i = 1, n
        do j = 1, i
            call random_number(rand)
            ! symmetric matrix
            a(i, j) = rand - 0.5d0
            a(j, i) = rand - 0.5d0
        end do
    end do
    aa = a

    lwork = -1
    CALL DGEEV('N', 'N', n, aa, n, wr, wi, vl, n, &
                vr, n, work, lwork, rc )
    lwork = MIN( lmax, INT( work( 1 ) ) )
    CALL DGEEV('N', 'N', n, aa, n, wr, wi, vl, n, &
                vr, n, work, lwork, rc )
!
    ! write eigenvalue
    ! write(*, *) 'Eigenvalue'
    ! do i = 1, n
        ! write(*, *) i, wr(i), wi(i)
    ! end do
    ! write(*, *) repeat('=', 70)

    open(66, file='a.dat', status='unknown')
    lambda = 0.d0
    do i = 1, n
        lambda(i) = complex(wr(i), wi(i))
        write(66, *) wr(i), wi(i)
    end do

    ! write eigenvector
    ! write(*, *) 'Eigenvector'
    ! do i = 1, n
        ! j = 1
        ! do while( j.le.n )
            ! if( wi(j).eq.0.d0) then
                ! write(*, 9998, advance='no') vr( i, j )
                ! v(i, j) = complex(vr(i, j), 0.d0)
                ! j = j + 1
            ! else
                ! write(*, 9999, advance='no') vr( i, j ), vr( i, j+1 )
                ! write(*,9999,advance='no') vr( i, j ), -vr( i, j+1 )
                ! v(i, j) = complex(vr(i, j), vr(i, j+1))
                ! v(i, j+1) = complex(vr(i, j), -vr(i, j+1))
                ! j = j + 2
            ! end if
        ! end do
        ! write(*,*)
    ! end do
    ! write(*, *) repeat('=', 70)
!
    ! write(*, *) 'check: A * v(j) - lambda(j) * v(j) ?= 0'
    ! do j = 1, 5
        ! print 9999, sum(matmul(a, v(:, j)) - lambda(j) * v(:, j))
    ! end do
    ! write(*, *) repeat('=', 70)

! 9998 FORMAT( 11(:,1X,F6.2) )
! 9999 FORMAT( 11(:,1X,'(',F6.2,',',F6.2,')') )
end program main
\end{minted}
\subsection*{4.f95}
\begin{minted}[frame=lines, linenos, fontsize=\large]
{fortran}
program FD
    ! solve y'' = p(x)y' + q(x)y + r(x)
    ! solve Aw = b
    implicit none
    integer, parameter :: N=101
    double precision, parameter :: pi=4.d0*atan(1.d0)
    double precision, external :: p, q, r
    double precision :: A(N, N), w(N), b(N)
    double precision :: xi, xf, h, wi, wf
    integer :: i
    double precision :: pivot(N), work2(N)
    integer :: rc
    ! if N is 1
    ! double precision :: A(N, N), w(N), b(N)
    
    ! bd condition
    wi = -.3d0
    wf = -.1d0
    xi = 0.d0
    xf = pi / 2.d0
    h = (xf - xi) / (dble(N) + 1.d0)

    ! A matrix
    ! first row
    A(1, 1) = 2.d0 + h**2 * q(xi+h)
    A(1, 2) = -1.d0 + .5d0 * h * p(xi+h)
    do i = 2, N-1
        A(i, i-1) = -1.d0 - .5d0 * h * p(xi+dble(i)*h)
        A(i, i) = 2.d0 + h**2 * q(xi+dble(i)*h)
        A(i, i+1) = -1.d0 + .5d0 * h * p(xi+dble(i)*h)
    end do
    ! last row
    A(N, N-1) = -1.d0 - .5d0 * h * p(xi+dble(N)*h)
    A(N, N) = 2.d0 + h**2 * q(xi+dble(N)*h)

    ! b matrix
    b(1) = -h**2 * r(xi+h) + (1.d0 + .5d0 * h * p(xi+h)) * wi
    do i = 2, N-1
        b(i) = -h**2 * r(xi+dble(i)*h)
    end do
    b(N) = -h**2 * r(xi+dble(N)*h) + (1.d0 + .5d0 * h * p(xi+dble(N)*h)) * wf

    ! find the inverse of A
    call dgetrf(N, N, A, N, pivot, rc)
    call dgetri(N, A, N, pivot, work2, N, rc)

    ! w = A^{-1}b
    w = matmul(A, b)

    ! case: if N = 1
    ! A(1, 1) = 2.d0 + h**2 * q(xi+h)
    ! b(1) = -h**2 * r(xi+h) + (1.d0 + .5d0 * h * p(xi+h)) * wi + (1.d0 - .5d0 * h * p(xi+h)) * wf
    ! w(1) = b(1) / A(1, 1)



    ! output
    open(66, file='d.dat', status='unknown')
    write(66, *) wi
    do i = 1, N
        write(66, *) w(i)
    end do
    write(66, *) wf
    close(66)
end program FD

function p(x)
    implicit none
    double precision :: p, x
    p = 1.d0
end function

function q(x)
    implicit none
    double precision :: q, x
    q = 2.d0
end function

function r(x)
    implicit none
    double precision :: r, x
    r = dcos(x)
end function