Line data    Source code

       1             : !/*************************************************************************
       2             : ! *
       3             : ! *  Project
       4             : ! *                         _____ _____  __  __ _____
       5             : ! *                        / ____|  __ \|  \/  |  __ \
       6             : ! *  ___  _ __   ___ _ __ | |  __| |__) | \  / | |__) |
       7             : ! * / _ \| '_ \ / _ \ '_ \| | |_ |  ___/| |\/| |  ___/
       8             : ! *| (_) | |_) |  __/ | | | |__| | |    | |  | | |
       9             : ! * \___/| .__/ \___|_| |_|\_____|_|    |_|  |_|_|
      10             : ! *      | |
      11             : ! *      |_|
      12             : ! *
      13             : ! * Copyright (C) Akiel Aries, <akiel@akiel.org>, et al.
      14             : ! *
      15             : ! * This software is licensed as described in the file LICENSE, which
      16             : ! * you should have received as part of this distribution. The terms
      17             : ! * among other details are referenced in the official documentation
      18             : ! * seen here : https://akielaries.github.io/openGPMP/ along with
      19             : ! * important files seen in this project.
      20             : ! *
      21             : ! * You may opt to use, copy, modify, merge, publish, distribute
      22             : ! * and/or sell copies of the Software, and permit persons to whom
      23             : ! * the Software is furnished to do so, under the terms of the
      24             : ! * LICENSE file. As this is an Open Source effort, all implementations
      25             : ! * must be of the same methodology.
      26             : ! *
      27             : ! *
      28             : ! *
      29             : ! * This software is distributed on an AS IS basis, WITHOUT
      30             : ! * WARRANTY OF ANY KIND, either express or implied.
      31             : ! *
      32             : ! ************************************************************************/
      33             : ! linsys.f90
      34             : 
      35           0 : subroutine solveLinearSystem(A, b, x, n)
      36             :     implicit none
      37             :     integer, intent(in) :: n
      38             :     real, dimension(n, n), intent(inout) :: A
      39             :     real, dimension(n), intent(inout) :: b
      40             :     real, dimension(n), intent(out) :: x
      41             :     real :: factor
      42             :     integer :: i, j, k
      43             : 
      44             :     ! Forward elimination
      45           0 :     do k = 1, n - 1
      46           0 :         do i = k + 1, n
      47           0 :             factor = A(i, k)/A(k, k)
      48           0 :             A(i, k + 1:n) = A(i, k + 1:n) - factor*A(k, k + 1:n)
      49           0 :             b(i) = b(i) - factor*b(k)
      50             :         end do
      51             :     end do
      52             : 
      53             :     ! Back substitution
      54           0 :     x(n) = b(n)/A(n, n)
      55           0 :     do i = n - 1, 1, -1
      56           0 :         x(i) = (b(i) - dot_product(A(i, i + 1:n), x(i + 1:n)))/A(i, i)
      57             :     end do
      58             : 
      59           0 : end subroutine solveLinearSystem
      60             : 
      61             : !program main
      62             : !    implicit none
      63             : !    integer, parameter :: n = 3
      64             : !    real :: A(n, n), b(n), x(n)
      65             : !    integer :: i
      66             : !
      67             : ! Initialize matrix A and vector b
      68             : !    A = reshape([1.0, 2.0, -1.0, 3.0, 1.0, 1.0, 2.0, -1.0, 2.0], [n, n])
      69             : !    b = [3.0, 4.0, 5.0]
      70             : 
      71             : ! Solve the linear system
      72             : !    call solveLinearSystem(A, b, x, n)
      73             : 
      74             : ! Output the solution
      75             : !    write (*, *) "Solution:"
      76             : !    do i = 1, n
      77             : !        write (*, *) "x(", i, ") = ", x(i)
      78             : !    end do
      79             :