tools_asa007.F90

Go to the documentation of this file.

# 1 "/Users/miesch/JEDI/code/working_copy/internal/mpas-bundle/saber/src/saber/external/tools_asa007.fypp"
# 1 "/Users/miesch/JEDI/code/working_copy/internal/mpas-bundle/saber/src/saber/external/../instrumentation.fypp" 1
# 1 "/Users/miesch/JEDI/code/working_copy/internal/mpas-bundle/saber/src/saber/external/../subr_list.fypp" 1
!----------------------------------------------------------------------
! Header: subr_list
!> Subroutines/functions list
! Author: Benjamin Menetrier
! Licensing: this code is distributed under the CeCILL-C license
! Copyright 2015-... UCAR, CERFACS, METEO-FRANCE and IRIT
!----------------------------------------------------------------------
 
# 926 "/Users/miesch/JEDI/code/working_copy/internal/mpas-bundle/saber/src/saber/external/../subr_list.fypp"
# 2 "/Users/miesch/JEDI/code/working_copy/internal/mpas-bundle/saber/src/saber/external/../instrumentation.fypp" 2
!----------------------------------------------------------------------
! Header: instrumentation
!> Instrumentation functions
! Author: Benjamin Menetrier
! Licensing: this code is distributed under the CeCILL-C license
! Copyright 2015-... UCAR, CERFACS, METEO-FRANCE and IRIT
!----------------------------------------------------------------------
 
# 112 "/Users/miesch/JEDI/code/working_copy/internal/mpas-bundle/saber/src/saber/external/../instrumentation.fypp"
# 2 "/Users/miesch/JEDI/code/working_copy/internal/mpas-bundle/saber/src/saber/external/tools_asa007.fypp" 2
!----------------------------------------------------------------------
! Module: tools_asa007
!> Inverse of symmetric positive definite matrix routines
! Source: https://people.sc.fsu.edu/~jburkardt/f_src/asa007/asa007.html
! Author: Michael Healy
! Original licensing: none
! Modified by Alan Miller
! Fortran 90 version by John Burkardt
! Modified by Benjamin Menetrier for BUMP
! Licensing: this code is distributed under the CeCILL-C license
! Copyright 2015-... UCAR, CERFACS, METEO-FRANCE and IRIT
!----------------------------------------------------------------------
module tools_asa007
 
use tools_const, only: zero,one
use tools_kinds, only: kind_real
use tools_repro, only: inf,infeq
use type_mpl, only: mpl_type
 
 
implicit none
 
real(kind_real),parameter :: eta = 1.0e-9_kind_real !< Small parameter for the Cholesky decomposition
 
interface cholesky
   module procedure asa007_cholesky
end interface
interface syminv
   module procedure asa007_syminv
end interface
 
private
public :: cholesky,syminv
 
contains
 
!----------------------------------------------------------------------
! Subroutine: asa007_cholesky
!> Compute cholesky decomposition
!----------------------------------------------------------------------
subroutine asa007_cholesky(mpl,n,nn,a,u)
 
implicit none
 
! Passed variables
type(mpl_type),intent(inout) :: mpl  !< MPI data
integer,intent(in) :: n              !< Matrix rank
integer,intent(in) :: nn             !< Half-matrix size (n*(n-1)/2)
real(kind_real),intent(in) :: a(nn)  !< Matrix
real(kind_real),intent(out) :: u(nn) !< Matrix square-root
 
! Local variables
integer :: i,icol,ii,irow,j,k,kk,l,m
real(kind_real) :: w,x
 
! Set name
 
 
! Probe in
 
 
! Initialization
ii = 0
j = 1
k = 0
if (nn/=(n*(n+1))/2) call mpl%abort('asa007_cholesky','wrong size in Cholesky decomposition')
w = zero
 
! Factorize column by column, ICOL = column number
do icol=1,n
   ii = ii+icol
   x = eta**2*a(ii)
   l = 0
   kk = 0
 
   ! IROW = row number within column ICOL
   do irow=1,icol
      kk = kk+irow
      k = k+1
      w = a(k)
      m = j
      do i=1,irow-1
        l = l+1
        w = w-u(l)*u(m)
        m = m+1
      end do
      l = l+1
      if (irow==icol) exit
      if (abs(u(l))>zero) then
         u(k) = w/u(l)
      else
         u(k) = zero
         if (inf(abs(x*a(k)),w**2)) call mpl%abort('asa007_cholesky','matrix is not positive semi-definite')
      end if
   end do
 
   ! End of row, estimate relative accuracy of diagonal element
   if (infeq(abs(w),abs(eta*a(k)))) then
      u(k) = zero
   else
      if (w<zero) call mpl%abort('asa007_cholesky','matrix is not positive semi-definite')
      u(k) = sqrt(w)
   end if
   j = j+icol
end do
 
! Probe out
 
 
end subroutine asa007_cholesky
 
!----------------------------------------------------------------------
! Subroutine: asa007_syminv
!> Compute inverse of a symmetric matrix
!----------------------------------------------------------------------
subroutine asa007_syminv(mpl,n,nn,a,c)
 
implicit none
 
! Passed variables
type(mpl_type),intent(inout) :: mpl  !< MPI data
integer,intent(in) :: n              !< Matrix rank
integer,intent(in) :: nn             !< Half-matrix size (n*(n-1)/2)
real(kind_real),intent(in) :: a(nn)  !< Matrix
real(kind_real),intent(out) :: c(nn) !< Matrix inverse
 
! Local variables
integer :: i,icol,irow,j,jcol,k,l,mdiag,ndiag,nrow
real(kind_real) :: w(n),x
 
! Set name
 
 
! Probe in
 
 
! Initialization
nrow = n
if (nn/=(n*(n+1))/2) then
   call mpl%abort('asa007_syminv','wrong size in matrix inversion')
end if
w = zero
 
! Compute the Cholesky factorization of A
call cholesky(mpl,n,nn,a,c)
 
! Invert C and form the product (Cinv)' * Cinv, where Cinv is the inverse of C, row by row starting with the last row
irow = nrow
ndiag = nn
 
do
   if (abs(c(ndiag))>zero) then
      ! General case
      l = ndiag
      do i=irow,nrow
         w(i) = c(l)
         l = l+i
      end do
      icol = nrow
      jcol = nn
      mdiag = nn
      do
         l = jcol
         if (icol==irow) then
            x = one/w(irow)
         else
            x = zero
         end if
         k = nrow
         do while (irow<k)
            x = x-w(k)*c(l)
            k = k-1
            l = l-1
            if (mdiag<l) l = l-k+1
         end do
         c(l) = x/w(irow)
         if (icol<=irow) exit
         mdiag = mdiag-icol
         icol = icol-1
         jcol = jcol-1
      end do
   else
      ! Special case, zero diagonal element
      l = ndiag
      do j=irow,nrow
         c(l) = zero
         l = l+j
      end do
   end if
   ndiag = ndiag-irow
   irow = irow-1
   if (irow<=0) exit
end do
 
! Probe out
 
 
end subroutine asa007_syminv
 
end module tools_asa007