algebra_lin.F90

Go to the documentation of this file.

module algebra_lin

  use real_type, only: rp
  use basic_tools

  implicit none  
  private

  !       %----------------------------------------%
  !       |                                        |
  !       |          PUBLIC DATA                   |
  !       |                                        |
  !       %----------------------------------------%
   
  !!
  !! LINEAR ALGEBRA
  !!
  !! GENERIC
  public :: transpose          ! tested
  public :: matvecprod         ! tested         
  !!
  !! NON GENERIC
  public :: matmatprod         ! tested  
  public :: invtrisup          ! tested  

  !!
  !! LINEAR ALGEBRA in R**3
  !!
  !!
  public :: nrm2_3d           ! tested
  public :: matvecprod_3x3    ! tested
  public :: matvecprod_3x2    ! tested
  public :: matvecprod_3x1    ! tested
  public :: solve_3x3         ! tested
  public :: det_3x3           ! tested  
  public :: crossprod_3d      ! tested
  public :: cplmtmat_3x2      ! tested
  public :: cplmtmat_3x3      ! tested
  public :: det_3x2           ! tested


  !       %----------------------------------------%
  !       |                                        |
  !       |       GENERIc SUBROUTINES              |
  !       |                                        |
  !       %----------------------------------------%

  interface matvecprod
     module procedure algebra_lin_matvecprod   !! TESTED
  end interface matvecprod

  interface transpose
     module procedure array_transpose      !! TESTED
  end interface transpose

contains

  
  subroutine algebra_lin_matvecprod(y,A,X)
    real(RP), dimension(:)  , intent(in) :: X
    real(RP), dimension(:)  , intent(out):: y
    real(RP), dimension(size(y,1), size(x,1)), intent(in) :: A

    integer :: ii

#if DBG>0
    if ( size(a,2)/= size(x,1) )       &
         &  call quit( &
            &  'algebra_lin: algebra_lin_matVecProd: 1' )
    if ( size(a,1)/= size(y,1) )       &
         &  call quit( &
            &  'algebra_lin: algebra_lin_matVecProd: 2' )
#endif
    
    y = x(1)*a(:,1)
    do ii=2, size(x,1)
       y = y + x(ii)*a(:,ii)
    end do

  end subroutine algebra_lin_matvecprod

  subroutine  matmatprod(M, A, B)
    real(RP), dimension(:,:), intent(out) :: M
    real(RP), dimension(:,:), intent(in)  :: A, B

    integer :: ii, jj

#if DBG>0
    if ( size(a,2)/= size(b,1) )       &
         &  call quit( &
            &  'algebra_lin: matMatProd: 1' )
    if ( size(a,1)/= size(m,1) )       &
         &  call quit( &
            &  'algebra_lin: matMatProd: 2' )
    if ( size(b,2)/= size(m,2) )       &
         &  call quit( &
            &  'algebra_lin: matMatProd: 3' )
#endif

    do ii=1, size(a,1)
       do jj=1, size(b,2)
          m(ii,jj) = sum(a(ii,:)*b(:,jj))
       end do
    end do

  end subroutine matmatprod

  subroutine  array_transpose(M, A)
    real(RP), dimension(:,:), intent(out) :: M
    real(RP), dimension(:,:), intent(in)  :: A

    integer :: ii, jj

#if DBG>0
    if ( size(a,1)/= size(m,2) ) call quit( &
            &  'algebra_lin: array_transpose: 1' )
    if ( size(a,2)/= size(m,1) ) call quit( &
            &  'algebra_lin: array_transpose: 2' )
#endif

    do ii=1, size(a,1)
       do jj=1, size(a,2)
          m(jj,ii) = a(ii,jj)
       end do
    end do

  end subroutine array_transpose


  subroutine invtrisup(res,M,vec)

    real(RP), dimension(:,:), intent(in)  :: M
    real(RP), dimension(:)  , intent(in)  :: vec
    real(RP), dimension(:)  , intent(out) :: res

    integer :: ii,ji,ni

#if DBG>0
    if ( size(m,1)/= size(m,2) )   call quit( &
            &  'algebra_lin: invtrisup: 1' )
    if ( size(m,1)/= size(vec,1) ) call quit( &
            &  'algebra_lin: invtrisup: 2' )
    if ( size(m,1)/= size(res,1) ) call quit( &
            &  'algebra_lin: invtrisup: 3' )
#endif

    ni=size(m,1)
    res(ni)=vec(ni)/m(ni,ni)

    do ii=ni-1,1,-1
       res(ii)=vec(ii)
       do ji=ii+1,ni
          res(ii)=res(ii)-m(ii,ji)*res(ji)
       end do
       res(ii)=res(ii)/m(ii,ii)
    end do

  end subroutine invtrisup

  function crossprod_3d(u,v) result(res)
    real(RP), dimension(3), intent(in) :: u
    real(RP), dimension(3), intent(in) :: v
    real(RP), dimension(3)             :: res

    res(1) =  u(2)*v(3)-u(3)*v(2)
    res(2) = -u(1)*v(3)+u(3)*v(1)
    res(3) =  u(1)*v(2)-u(2)*v(1)

  end function crossprod_3d

  function det_3x3(M) result(d)
    real(RP), dimension(3,3), intent(in) :: M
    real(RP)                             :: d

    real(RP) :: u, v, w

    u = m(2,2)*m(3,3)-m(3,2)*m(2,3)
    v = m(1,2)*m(3,3)-m(1,3)*m(3,2)
    w = m(1,2)*m(2,3)-m(1,3)*m(2,2)

    d = m(1,1)*u - m(2,1)*v + m(3,1)*w

  end function det_3x3


  function det_3x2(M) result(d)
    real(RP), dimension(3,2), intent(in) :: M
    real(RP)                             :: d

    real(RP) :: u

    u = m(2,1)*m(3,2)-m(3,1)*m(2,2)
    d = u**2
    u = -m(1,1)*m(3,2)+m(3,1)*m(1,2)
    d = d + u**2
    u =  m(1,1)*m(2,2)-m(2,1)*m(1,2)
    d = d + u**2
    d = sqrt(d)

  end function det_3x2



  function cplmtmat_3x3(M) result(N)
    real(RP), dimension(3,3)             :: N
    real(RP), dimension(3,3), intent(in) :: M

    n(:,1) = crossprod_3d(m(:,2), m(:,3))
    n(:,2) = crossprod_3d(m(:,3), m(:,1))
    n(:,3) = crossprod_3d(m(:,1), m(:,2))

  end function cplmtmat_3x3

  function cplmtmat_3x2(M) result(N)
    real(RP), dimension(3,2)             :: N
    real(RP), dimension(3,2), intent(in) :: M

    real(RP), dimension(3) :: e1, e2, e3
    real(RP):: xB, xC, yC

    e1 = m(:,1) 
    e2 = m(:,2)
    e3 = crossprod_3d( e1, e2 )
    e2 = crossprod_3d(e3, e1)
    e1 = e1 / nrm2_3d(e1)
    e2 = e2 / nrm2_3d(e2)
    xb = dot_product(m(:,1), e1)
    xc = dot_product(m(:,2), e1)
    yc = dot_product(m(:,2), e2)

    n(:,1) = yc*e1 - xc*e2
    n(:,2) = xb*e2

  end function cplmtmat_3x2

  function nrm2_3d(v) result(x)
    real(RP), dimension(3) :: v
    real(RP) :: x

    x = sum(v*v)
    x = sqrt(x)

  end function nrm2_3d

  subroutine matvecprod_3x3(y,A,X)
    real(RP), dimension(3)  , intent(out):: y
    real(RP), dimension(3)  , intent(in) :: X
    real(RP), dimension(3,3), intent(in) :: A

    y(1:3) =          a(1:3, 1) * x(1)
    y(1:3) = y(1:3) + a(1:3, 2) * x(2)
    y(1:3) = y(1:3) + a(1:3, 3) * x(3)

  end subroutine matvecprod_3x3


  
  subroutine matvecprod_3x2(y,A,X)
    real(RP), dimension(3)  , intent(out):: y
    real(RP), dimension(2)  , intent(in) :: X
    real(RP), dimension(3,2), intent(in) :: A

    y(1:3) =          a(1:3, 1) * x(1)
    y(1:3) = y(1:3) + a(1:3, 2) * x(2)
    
  end subroutine matvecprod_3x2

  
  subroutine matvecprod_3x1(y,A,X)
    real(RP), dimension(3)  , intent(out):: y
    real(RP), dimension(1)  , intent(in) :: X
    real(RP), dimension(3,1), intent(in) :: A

    y(1:3) = a(1:3, 1) * x(1)

  end subroutine matvecprod_3x1


  subroutine solve_3x3(x, A, b)
    real(RP), dimension(3)  , intent(out) :: x
    real(RP), dimension(3,3), intent(in)  :: A
    real(RP), dimension(3)  , intent(in)  :: b

    real(RP), dimension(3) :: v

    v    = crossprod_3d(a(:,2),a(:,3))
    x(1) = dot_product(b,v) / dot_product(a(:,1),v)

    v    = crossprod_3d(a(:,1),a(:,3))
    x(2) = dot_product(b,v) / dot_product(a(:,2),v)

    v    = crossprod_3d(a(:,1),a(:,2))
    x(3) = dot_product(b,v) / dot_product(a(:,3),v)

  end subroutine solve_3x3

end module algebra_lin