Merge branch 'master' into reconnection

sources/algebra.F90

@@ -51,7 +51,7 @@ module algebra
   public :: quadratic, quadratic_normalized
   public :: cubic, cubic_normalized
   public :: quartic
-  public :: tridiag, invert
+  public :: tridiag, pentadiag, invert
 
 !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 !
@@ -1148,6 +1148,87 @@ module algebra
 !
 !===============================================================================
 !
+! subroutine PENTADIAG:
+! --------------------
+!
+!   Subroutine solves the system of pentadiagonal linear equations using
+!   PTRANS-I algorithm.
+!
+!   Arguments:
+!
+!     n    - the size of the matrix;
+!     e(:) - the vector of the outer lower off-diagonal values;
+!     c(:) - the vector of the inner lower off-diagonal values;
+!     d(:) - the vector of the diagonal values;
+!     a(:) - the vector of the inner upper off-diagonal values;
+!     b(:) - the vector of the outer upper off-diagonal values;
+!     y(:) - the vector of the right side values;
+!     x(:) - the solution vector;
+!
+!   References:
+!
+!     [1] Askar, S. S., and Karawia, A. A.,
+!         "On Solving Pentadiagonal Linear Systems via Transformations",
+!         Mathematical Problems in Engineering,
+!         2015, Volume 2015, Article ID 232456, 9 pages
+!         http://dx.doi.org/10.1155/2015/232456
+!
+!===============================================================================
+!
+  subroutine pentadiag(n, e, c, d, a, b, y, x)
+
+    implicit none
+
+    integer                   , intent(in)  :: n
+    real(kind=8), dimension(n), intent(in)  :: e, c, d, a, b, y
+    real(kind=8), dimension(n), intent(out) :: x
+
+    integer                    :: i
+    real(kind=8), dimension(n) :: al, bt, gm, mu, z
+!
+!-------------------------------------------------------------------------------
+!
+    mu(1) = d(1)
+    al(1) = a(1) / mu(1)
+    bt(1) = b(1) / mu(1)
+    z (1) = y(1) / mu(1)
+
+    gm(2) = c(2)
+    mu(2) = d(2) - al(1) * gm(2)
+    al(2) = (a(2) - bt(1) * gm(2)) / mu(2)
+    bt(2) = b(2) / mu(2)
+    z (2) = (y(2) - z(1) * gm(2)) / mu(2)
+
+    do i = 3, n-2
+      gm(i) = c(i) - al(i-2) * e(i)
+      mu(i) = d(i) - bt(i-2) * e(i) - al(i-1) * gm(i)
+      al(i) = (a(i) - bt(i-1) * gm(i)) / mu(i)
+      bt(i) = b(i) / mu(i)
+      z (i) = (y(i) - z(i-2) * e(i) - z(i-1) * gm(i)) / mu(i)
+    end do
+
+    gm(n-1) = c(n-1) - al(n-3) * e(n-1)
+    mu(n-1) = d(n-1) - bt(n-3) * e(n-1) - al(n-2) * gm(n-1)
+    al(n-1) = (a(n-1) - bt(n-2) * gm(n-1)) / mu(n-1)
+
+    gm(n) = c(n) - al(n-2) * e(n)
+    mu(n) = d(n) - bt(n-3) * e(n) - al(n-1) * gm(n)
+
+    z(n-1) = (y(n-1) - z(n-2) * e(n-1) - z(n-2) * gm(n-1)) / mu(n-1)
+    z(n  ) = (y(n) - z(n-1) * e(n) - z(n-1) * gm(n)) / mu(n)
+
+    x(n  ) = z(n)
+    x(n-1) = z(n-1) - al(n-1) * x(n)
+    do i = n - 2, 1, -1
+      x(i) = z(i) - al(i) * x(i+1) - bt(i) * x(i+2)
+    end do
+
+!-------------------------------------------------------------------------------
+!
+  end subroutine pentadiag
+!
+!===============================================================================
+!
 ! subroutine INVERT:
 ! -----------------
 !