# Solving Linear Algebra Problem
## Linear equation
Linear equations try to solve following kind of equation
\begin{equation}
\textbf{A}\textbf{x} = \textbf{b}
\tag{1}
\end{equation}
where $\textbf{A}$ is a NxN matrix, $\textbf{x}$ is a Nx1 and $\textbf{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 $\textbf{x}$ immediately. Second method is that we can find the inverse matrix of $\textbf{A}$ , then multiplied by $\textbf{b}$. Finding inverse matrix in LAPACK require two subroutine, $dgetrf$ and $dgetri$.
### Example 1
In the first example we first solve $\textbf{x}$ and then verified the result by multiply $\textbf{x}$ by the matrix $\textbf{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 $\textbf{A}\textbf{x}$ and
$\textbf{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 $\textbf{A}$ is a symmetric matrix the later method has a higher accuracy despite a slower speed.
### Code
```fortran=
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 $\textbf{A}$. Following print the result of five eigenvalue and five eigenvector. By definition the relation between eigenvalue and eigenvector is given by
\begin{equation}
\textbf{A} * v(j) = \lambda(j) * v(j)
\tag{2}
\end{equation}
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{equation}
P(\lambda) = \dfrac{2}{\pi R^2}\sqrt{R^2-x^2}
\tag{3}
\end{equation}
### Code
```fortran=
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.
```fortran=
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
```
Sign in with Wallet
Connect another wallet