View difference between Paste ID: 4ZXMYJKJ and g7GkWMwa

-module global_vars
+subroutine LU(n, m, A, P, Y)
-    integer :: N = 100, l = 1
+    integer, intent(in) :: n, m
-    real :: R = 15.0, V_0 = 4.0, a = 5.0, lk = 0.05, rk = 10.0, dk = 0.05, k
+    integer :: i, j, k, pivot, temp
-    save
+    double precision :: pivotValue
-end module
+    integer, dimension(:), intent(inout) :: P
-
+    double precision, dimension(:), allocatable :: x
-function V(x)
+    double precision, dimension(:, :), intent(inout) :: A, Y
-    use global_vars
+    allocate(x(n))
     do i = 1, n
-    real :: V
+        !searching for max
-    real, intent(in) :: x
+        pivotValue = 0.0
-    if (x < a) then
+        pivot = -1
-        V = -V_0
+        do k = i, n
-    else
+            if (abs(A(k, i)) > pivotValue) then
-        V = 0.0
+                pivotValue = abs(A(k, i))
                 pivot = k
-    return
+            end if
-end function
+
-
+        !is matrix invertible?
-function F(x)
+        if (abs(pivotValue) < 1.0/10000) then
-    use global_vars
+            print *, "matrix is invertiblle"
         end if
-    real :: F, V
+        !row swap
-    real, intent(in) :: x
+        !call dswap(1, P(i), 0, P(pivot), 0)
-    F = V(x) + l * (l + 1) / x**2 - k**2
+        temp = P(i)
-    return
+        P(i) = P(pivot)
-end function
+        P(pivot) = temp
-
+        call dswap(n, A(i,1:n), 1, A(pivot, 1:n), 1)
-recursive function jl(x, l) result(ans)
+        do j = i + 1, n
-    real, intent(in) :: x
+            A(j, i) = A(j, i) / A(i, i)
-    real :: ans
+            do k = i + 1, n
-    integer, intent(in) :: l
+                A(j, k) = A(j, k) - A(j, i) * A(i, k)
-    if (l == 0) then
+            end do
-        ans = sin(x) / x
+
-    else if (l == 1) then
+
-        ans = sin(x) / x**2 - cos(x) / x
+    !A = L + U - I
-    else
+    !Gauss
-        ans = (2 * (l - 1) + 1) / x * jl(x, l - 1) - jl(x, l - 2)
+    do k = 1, m
         !forward
-end function
+        do i = 1, n
             x(i) = Y(P(i), k)
-function jl_hat(x, l)
+            x(i) = x(i) - ddot(i-1, A(i,1:i-1), 1, x(1:i-1), 1)
-    real, intent(in) :: x
+            !do j = 1, i - 1
-    integer, intent(in) :: l
+            !    x(i) = x(i) - A(i, j) * x(j)
-    real :: jl, jl_hat
+            !end do
-    jl_hat = jl(x, l) * x
+
-    return
+        !backward
-end function
+        do i = n, 1, -1
             if (i < n) then
-subroutine Simpson(h, f, N, ans)
+                x(i) = x(i) - ddot(n-i, x(i+1:n), 1, A(i,i+1:n), 1)
             end if
-    integer :: i, N
+            !do j = i + 1, n
-    real :: h
+            !    x(i) = x(i) - A(i, j) * x(j)
-    complex :: f(N), ans
+            !end do
-    ans = 0.0
+            x(i) = x(i) / A(i, i)
-    do i = 2, N - 1, 2
+
-        ans = ans + f(i - 1) + 4.0 * f(i) + f(i + 1)
+        !answer
         do i = 1, n
-    ans = ans * h / 3.0
+            Y(i, k) = x(i)
         end do
-
+
-program scattering
+    deallocate(x)
-    use global_vars
+
 
-    integer :: i, j, out = 42
+program lu_decomposition
-    real :: F, V, h, jl_hat,  pi = 3.1415926
+
-    complex :: Al, delta
+    integer :: n, m, i, inp, outp
-    integer, allocatable :: ipiv(:)
+    integer, dimension(:), allocatable :: P, sP
-    complex, allocatable :: X(:),  mA(:, :)
+    double precision, dimension(:,:), allocatable :: A, sA, tempY, Y, sY
-    allocate(X(N))
+    interface
-    allocate(ipiv(N))
+        subroutine LU(n, m, A, P, Y)
-    allocate(mA(N, N))
+            integer, intent(in) :: n, m
-    h = R / N
+            integer :: i, j, k, pivot, temp
-    k = lk
+            double precision :: pivotValue
-    open(unit = out, file = "output.txt")
+            integer, dimension(:), intent(inout) :: P
-    do while (k <= rk)
+            double precision, dimension(:), allocatable :: x
-        !mA = -d^2/dr^2 + F(r) - k^2
+            double precision, dimension(:, :), intent(inout) :: A, Y
-        mA(:, :) = 0.0
+        end subroutine LU
-        do i = 1, N
+    end interface
-            if(i > 1) mA(i, i - 1) = 12.0 / h**2  - F(h * (i - 1))
+    !input
-            mA(i, i) = -24.0 / h**2 - 10.0 * F(h * i)
+    inp = 42
-            if(i < N) mA(i, i + 1) = 12.0 / h**2 - F(h * (i + 1))
+    outp = 43
     open(unit = inp, file = "input.txt")
-        mA(N, N) = mA(N, N) + exp((0.0, 1.0) * k * h) * (12.0 / h**2 - F(h * (N + 1)))
+    read(inp, *) n, m
-        !X = -V(r)jl(kr)
+    allocate(P(n))
-        do i = 1, N
+    do i = 1, n
-            X(i) = 10.0 * V(h * i) * jl_hat(k * i * h, l)
+        P(i) = i
-            if(i > 1) X(i) = X(i) + V(h * (i - 1)) * jl_hat(k * (i - 1) * h, l)
+
-            if(i < n) X(i) = X(i) + V(h * (i + 1)) * jl_hat(k * (i + 1) * h, l)
+    allocate(sP(n))
     sP = P
-        !X = chi_l
+    allocate(A(n, n))
-        call cgesv(N, 1, mA, N, ipiv, X, N, i)
+    allocate(sA(n, n))
-        !X = jl(kr)V(r)f_l(r, k) = jl(kr)V(r)(jl(kr) + chi_l(r, k))
+    allocate(Y(n, m))
-        do i =1,n
+    allocate(sY(n, m))
-            X(i) = jl_hat(k * h * i, l) * V(h * i) * (jl_hat(k * h * i, l) + X(i))
+    allocate(tempY(m, n))
     read (inp, *) A
-        call Simpson(h, X, N, Al)
+    A = transpose(A)
-        Al = -1.0 * Al / k**2
+    sA = A
-        delta = log(1 + (0.0, 2.0) * k * Al) / (0.0, 2.0)
+    read (inp, *) tempY
-        write(out, *) k, real(delta)
+    Y = transpose(tempY)
-        k = k + dk
+    sY = Y
     close(inp)
-    close(out)
+    open(unit = outp, file = "output.txt")
     !LU decomposition
     call LU(n, m, A, P, Y)
     !LU output
     write(outp, *) "A"
     do i = 1, n
         write(outp, *) A(i, :)
     end do
     write(outp, *) "P"
     write(outp, *) P
     write(outp, *) "Y"
     do i = 1, n
         write(outp, "(1000F10.3)") Y(i, :)
     end do
     !sgesv
     call dgesv(n, m, sA, n, sP, sY, n, i)
     !sgesv output
     if (i .eq. 0) then
         print *, "successful sgesv execution"
     end if
     write(outp, *) "A by sgesv"
     do i = 1, n
         write(outp, *) sA(i, :)
     end do
     write(outp, *) "P by sgesv"
     write(outp, *) sP
     write(outp, *) "Y by sgesv"
     do i = 1, n
         write(outp, "(1000F10.3)") sY(i, :)
     end do
     close(outp)
     !memory
     deallocate(P)
     deallocate(sP)
     deallocate(A)
     deallocate(Y)
     deallocate(sY)
 end program