linpack.f90 Source File

module linpack
!! Selected routines from LINPACK.

   use daskr_kinds, only: rk, zero, one
   implicit none

contains

   pure subroutine dgbfa(abd, lda, n, ml, mu, ipvt, info)
   !! [[dgbfa]] factors a real band matrix by elimination.
      !
      !  Discussion:
      !
      !    DGBFA is usually called by DGBCO, but it can be called
      !    directly with a saving in time if RCOND is not needed.
      !
      !  Licensing:
      !
      !    This code is distributed under the GNU LGPL license.
      !
      !  Modified:
      !
      !    17 May 2005
      !
      !  Author:
      !
      !    FORTRAN90 version by John Burkardt.
      !
      !  Reference:
      !
      !    Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart,
      !    LINPACK User's Guide,
      !    SIAM, 1979,
      !    ISBN13: 978-0-898711-72-1,
      !    LC: QA214.L56.
      !
      !  Parameters:
      !
      !    Input/output, real ( kind = 8 ) ABD(LDA,N).  On input, the matrix in band
      !    storage.  The columns of the matrix are stored in the columns of ABD
      !    and the diagonals of the matrix are stored in rows ML+1 through
      !    2*ML+MU+1 of ABD.  On output, an upper triangular matrix in band storage
      !    and the multipliers which were used to obtain it.  The factorization
      !    can be written A = L*U where L is a product of permutation and unit lower
      !    triangular matrices and U is upper triangular.
      !
      !    Input, integer ( kind = 4 ) LDA, the leading dimension of the array ABD.
      !    2*ML + MU + 1 <= LDA is required.
      !
      !    Input, integer ( kind = 4 ) N, the order of the matrix.
      !
      !    Input, integer ( kind = 4 ) ML, MU, the number of diagonals below and above
      !    the main diagonal.  0 <= ML < N, 0 <= MU < N.
      !
      !    Output, integer ( kind = 4 ) IPVT(N), the pivot indices.
      !
      !    Output, integer ( kind = 4 ) INFO, error flag.
      !    0, normal value.
      !    K, if U(K,K) == 0.0.  This is not an error condition for this
      !      subroutine, but it does indicate that DGBSL will divide by zero if
      !      called.  Use RCOND in DGBCO for a reliable indication of singularity.

      use blas_interfaces, only: idamax, daxpy, dscal

      real(rk), intent(inout) :: abd(lda, n)
      integer, intent(in) :: lda
      integer, intent(in) :: n
      integer, intent(in) :: ml
      integer, intent(in) :: mu
      integer, intent(out) :: ipvt(n)
      integer, intent(out) :: info

      integer :: i
      integer :: i0
      integer :: j
      integer :: j0
      integer :: j1
      integer :: ju
      integer :: jz
      integer :: k
      integer :: l
      integer :: lm
      integer :: m
      integer :: mm
      real(rk) :: t

      m = ml + mu + 1
      info = 0
      !
      !  Zero initial fill-in columns.
      !
      j0 = mu + 2
      j1 = min(n, m) - 1

      do jz = j0, j1
         i0 = m + 1 - jz
         do i = i0, ml
            abd(i, jz) = zero
         end do
      end do

      jz = j1
      ju = 0
      !
      !  Gaussian elimination with partial pivoting.
      !
      do k = 1, n - 1
         !
         !  Zero out the next fill-in column.
         !
         jz = jz + 1
         if (jz <= n) then
            abd(1:ml, jz) = zero
         end if
         !
         !  Find L = pivot index.
         !
         lm = min(ml, n - k)
         l = idamax(lm + 1, abd(m, k), 1) + m - 1
         ipvt(k) = l + k - m
         !
         !  Zero pivot implies this column already triangularized.
         !
         if (abd(l, k) == zero) then

            info = k
            !
            !  Interchange if necessary.
            !
         else

            if (l /= m) then
               t = abd(l, k)
               abd(l, k) = abd(m, k)
               abd(m, k) = t
            end if
            !
            !  Compute multipliers.
            !
            t = -one/abd(m, k)
            call dscal(lm, t, abd(m + 1, k), 1)
            !
            !  Row elimination with column indexing.
            !
            ju = min(max(ju, mu + ipvt(k)), n)
            mm = m

            do j = k + 1, ju
               l = l - 1
               mm = mm - 1
               t = abd(l, j)
               if (l /= mm) then
                  abd(l, j) = abd(mm, j)
                  abd(mm, j) = t
               end if
               call daxpy(lm, t, abd(m + 1, k), 1, abd(mm + 1, j), 1)
            end do

         end if

      end do

      ipvt(n) = n

      if (abd(m, n) == zero) then
         info = n
      end if

   end subroutine dgbfa

   pure subroutine dgbsl(abd, lda, n, ml, mu, ipvt, b, job)
   !! [[dgbsl]] solves a real banded system factored by [[dgbco]] or [[dgbfa]].
      !
      !  Discussion:
      !
      !    DGBSL can solve either A * X = B  or  A' * X = B.
      !
      !    A division by zero will occur if the input factor contains a
      !    zero on the diagonal.  Technically this indicates singularity
      !    but it is often caused by improper arguments or improper
      !    setting of LDA.  It will not occur if the subroutines are
      !    called correctly and if DGBCO has set 0.0 < RCOND
      !    or DGBFA has set INFO == 0.
      !
      !    To compute inverse(A) * C  where C is a matrix with P columns:
      !
      !      call dgbco ( abd, lda, n, ml, mu, ipvt, rcond, z )
      !
      !      if ( rcond is too small ) then
      !        exit
      !      end if
      !
      !      do j = 1, p
      !        call dgbsl ( abd, lda, n, ml, mu, ipvt, c(1,j), 0 )
      !      end do
      !
      !  Licensing:
      !
      !    This code is distributed under the GNU LGPL license.
      !
      !  Modified:
      !
      !    17 May 2005
      !
      !  Author:
      !
      !    FORTRAN90 version by John Burkardt.
      !
      !  Reference:
      !
      !    Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart,
      !    LINPACK User's Guide,
      !    SIAM, 1979,
      !    ISBN13: 978-0-898711-72-1,
      !    LC: QA214.L56.
      !
      !  Parameters:
      !
      !    Input, real ( kind = 8 ) ABD(LDA,N), the output from DGBCO or DGBFA.
      !
      !    Input, integer ( kind = 4 ) LDA, the leading dimension of the array ABD.
      !
      !    Input, integer ( kind = 4 ) N, the order of the matrix.
      !
      !    Input, integer ( kind = 4 ) ML, MU, the number of diagonals below and above
      !    the main diagonal.  0 <= ML < N, 0 <= MU < N.
      !
      !    Input, integer ( kind = 4 ) IPVT(N), the pivot vector from DGBCO or DGBFA.
      !
      !    Input/output, real ( kind = 8 ) B(N).  On input, the right hand side.
      !    On output, the solution.
      !
      !    Input, integer ( kind = 4 ) JOB, job choice.
      !    0, solve A*X=B.
      !    nonzero, solve A'*X=B.

      use blas_interfaces, only: daxpy, ddot

      real(rk), intent(in) :: abd(lda, n)
      integer, intent(in) :: lda
      integer, intent(in) :: n
      integer, intent(in) :: ml
      integer, intent(in) :: mu
      integer, intent(in) :: ipvt(n)
      real(rk), intent(inout) :: b(n)
      integer, intent(in) :: job

      integer :: k
      integer :: l
      integer :: la
      integer :: lb
      integer :: lm
      integer :: m
      real(rk) :: t

      m = mu + ml + 1
      !
      !  JOB = 0, Solve A * x = b.
      !
      !  First solve L * y = b.
      !
      if (job == 0) then

         if (0 < ml) then

            do k = 1, n - 1
               lm = min(ml, n - k)
               l = ipvt(k)
               t = b(l)
               if (l /= k) then
                  b(l) = b(k)
                  b(k) = t
               end if
               call daxpy(lm, t, abd(m + 1, k), 1, b(k + 1), 1)
            end do

         end if
         !
         !  Now solve U * x = y.
         !
         do k = n, 1, -1
            b(k) = b(k)/abd(m, k)
            lm = min(k, m) - 1
            la = m - lm
            lb = k - lm
            t = -b(k)
            call daxpy(lm, t, abd(la, k), 1, b(lb), 1)
         end do
         !
         !  JOB nonzero, solve A' * x = b.
         !
         !  First solve U' * y = b.
         !
      else

         do k = 1, n
            lm = min(k, m) - 1
            la = m - lm
            lb = k - lm
            t = ddot(lm, abd(la, k), 1, b(lb), 1)
            b(k) = (b(k) - t)/abd(m, k)
         end do
         !
         !  Now solve L' * x = y.
         !
         if (0 < ml) then

            do k = n - 1, 1, -1
               lm = min(ml, n - k)
               b(k) = b(k) + ddot(lm, abd(m + 1, k), 1, b(k + 1), 1)
               l = ipvt(k)
               if (l /= k) then
                  t = b(l)
                  b(l) = b(k)
                  b(k) = t
               end if
            end do

         end if

      end if

   end subroutine dgbsl

   pure subroutine dgefa(a, lda, n, ipvt, info)
   !! [[dgefa]] factors a real general matrix.
      !
      !  Licensing:
      !
      !    This code is distributed under the GNU LGPL license.
      !
      !  Modified:
      !
      !    07 March 2001
      !
      !  Author:
      !
      !    FORTRAN90 version by John Burkardt.
      !
      !  Reference:
      !
      !    Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart,
      !    LINPACK User's Guide,
      !    SIAM, 1979,
      !    ISBN13: 978-0-898711-72-1,
      !    LC: QA214.L56.
      !
      !  Parameters:
      !
      !    Input/output, real ( kind = 8 ) A(LDA,N).
      !    On intput, the matrix to be factored.
      !    On output, an upper triangular matrix and the multipliers used to obtain
      !    it.  The factorization can be written A=L*U, where L is a product of
      !    permutation and unit lower triangular matrices, and U is upper triangular.
      !
      !    Input, integer ( kind = 4 ) LDA, the leading dimension of A.
      !
      !    Input, integer ( kind = 4 ) N, the order of the matrix A.
      !
      !    Output, integer ( kind = 4 ) IPVT(N), the pivot indices.
      !
      !    Output, integer ( kind = 4 ) INFO, singularity indicator.
      !    0, normal value.
      !    K, if U(K,K) == 0.  This is not an error condition for this subroutine,
      !    but it does indicate that DGESL or DGEDI will divide by zero if called.
      !    Use RCOND in DGECO for a reliable indication of singularity.

      use blas_interfaces, only: idamax

      real(rk), intent(inout) :: a(lda, n)
      integer, intent(in) :: lda
      integer, intent(in) :: n
      integer, intent(out) :: ipvt(n)
      integer, intent(out) :: info

      integer :: j
      integer :: k
      integer :: l
      real(rk) :: t
      !
      !  Gaussian elimination with partial pivoting.
      !
      info = 0

      do k = 1, n - 1
         !
         !  Find L = pivot index.
         !
         l = idamax(n - k + 1, a(k, k), 1) + k - 1
         ipvt(k) = l
         !
         !  Zero pivot implies this column already triangularized.
         !
         if (a(l, k) == zero) then
            info = k
            cycle
         end if
         !
         !  Interchange if necessary.
         !
         if (l /= k) then
            t = a(l, k)
            a(l, k) = a(k, k)
            a(k, k) = t
         end if
         !
         !  Compute multipliers.
         !
         t = -one/a(k, k)
         a(k + 1:n, k) = a(k + 1:n, k)*t
         !
         !  Row elimination with column indexing.
         !
         do j = k + 1, n
            t = a(l, j)
            if (l /= k) then
               a(l, j) = a(k, j)
               a(k, j) = t
            end if
            a(k + 1:n, j) = a(k + 1:n, j) + t*a(k + 1:n, k)
         end do

      end do

      ipvt(n) = n

      if (a(n, n) == zero) then
         info = n
      end if

   end subroutine dgefa

   pure subroutine dgesl(a, lda, n, ipvt, b, job)
   !! [[dgesl]] solves a real general linear system A * X = B.
      !
      !  Discussion:
      !
      !    DGESL can solve either of the systems A * X = B or A' * X = B.
      !
      !    The system matrix must have been factored by DGECO or DGEFA.
      !
      !    A division by zero will occur if the input factor contains a
      !    zero on the diagonal.  Technically this indicates singularity
      !    but it is often caused by improper arguments or improper
      !    setting of LDA.  It will not occur if the subroutines are
      !    called correctly and if DGECO has set 0.0 < RCOND
      !    or DGEFA has set INFO == 0.
      !
      !  Licensing:
      !
      !    This code is distributed under the GNU LGPL license.
      !
      !  Modified:
      !
      !    07 March 2001
      !
      !  Author:
      !
      !    FORTRAN90 version by John Burkardt.
      !
      !  Reference:
      !
      !    Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart,
      !    LINPACK User's Guide,
      !    SIAM, 1979,
      !    ISBN13: 978-0-898711-72-1,
      !    LC: QA214.L56.
      !
      !  Parameters:
      !
      !    Input, real ( kind = 8 ) A(LDA,N), the output from DGECO or DGEFA.
      !
      !    Input, integer ( kind = 4 ) LDA, the leading dimension of A.
      !
      !    Input, integer ( kind = 4 ) N, the order of the matrix A.
      !
      !    Input, integer ( kind = 4 ) IPVT(N), the pivot vector from DGECO or DGEFA.
      !
      !    Input/output, real ( kind = 8 ) B(N).
      !    On input, the right hand side vector.
      !    On output, the solution vector.
      !
      !    Input, integer ( kind = 4 ) JOB.
      !    0, solve A * X = B;
      !    nonzero, solve A' * X = B.

      real(rk), intent(in) :: a(lda, n)
      integer, intent(in) :: lda
      integer, intent(in) :: n
      integer, intent(in) :: ipvt(n)
      real(rk), intent(inout) :: b(n)
      integer, intent(in) :: job

      integer :: k
      integer :: l
      real(rk) :: t
      !
      !  Solve A * X = B.
      !
      if (job == 0) then

         do k = 1, n - 1

            l = ipvt(k)
            t = b(l)

            if (l /= k) then
               b(l) = b(k)
               b(k) = t
            end if

            b(k + 1:n) = b(k + 1:n) + t*a(k + 1:n, k)

         end do

         do k = n, 1, -1
            b(k) = b(k)/a(k, k)
            t = -b(k)
            b(1:k - 1) = b(1:k - 1) + t*a(1:k - 1, k)
         end do

      else
         !
         !  Solve A' * X = B.
         !
         do k = 1, n
            t = dot_product(a(1:k - 1, k), b(1:k - 1))
            b(k) = (b(k) - t)/a(k, k)
         end do

         do k = n - 1, 1, -1

            b(k) = b(k) + dot_product(a(k + 1:n, k), b(k + 1:n))
            l = ipvt(k)

            if (l /= k) then
               t = b(l)
               b(l) = b(k)
               b(k) = t
            end if

         end do

      end if

   end subroutine dgesl

end module linpack