compute_coeff Subroutine

private impure subroutine compute_coeff(this)

Uses

    • formatmul
    • forsvd
  • proc~~compute_coeff~~UsesGraph proc~compute_coeff forpca::tpca%compute_coeff formatmul formatmul proc~compute_coeff->formatmul forsvd forsvd proc~compute_coeff->forsvd

Type Bound

tpca

Arguments

Type IntentOptional Attributes Name
class(tpca), intent(inout) :: this

Calls

proc~~compute_coeff~~CallsGraph proc~compute_coeff forpca::tpca%compute_coeff eig eig proc~compute_coeff->eig mask mask proc~compute_coeff->mask order order proc~compute_coeff->order svd svd proc~compute_coeff->svd

Called by

proc~~compute_coeff~~CalledByGraph proc~compute_coeff forpca::tpca%compute_coeff proc~pca forpca::tpca%pca proc~pca->proc~compute_coeff program~benchmark benchmark program~benchmark->proc~pca program~test1 test1 program~test1->proc~pca

Source Code

   impure subroutine compute_coeff(this)
#else
   pure subroutine compute_coeff(this)
#endif
      use forsvd
      use formatmul, only: matmul

      class(tpca), intent(inout)               :: this
      real(rk), dimension(this%ncol)           :: mean
      real(rk), dimension(this%ncol,this%ncol) :: cov
      integer                                  :: i
      real(rk), dimension(:,:), allocatable    :: mean_data_T

      allocate(this%mean_data(this%nrow,this%ncol))
      allocate(mean_data_T(this%ncol,this%nrow))

#if defined(USE_COARRAY)
      if (this_image() == 1) then
#endif
         mean = sum(this%matrix, dim=1) / real(this%nrow, kind=rk)

         do i = 1, this%nrow
            this%mean_data(i, :) = this%matrix(i, :) - mean
         end do
         mean_data_T = transpose(this%mean_data)
#if defined(USE_COARRAY)
      end if
#endif

#if defined(USE_COARRAY)
      call co_broadcast(mean_data_T, source_image=1)
      call co_broadcast(this%mean_data, source_image=1)
      cov = matmul(mean_data_T, this%mean_data, method='coarray', option='m2')/real(this%nrow - 1, kind=rk)
#else
      cov = matmul(mean_data_T, this%mean_data, method='default', option='m2')/real(this%nrow - 1, kind=rk)
#endif

      allocate(this%latent(this%ncol))
      allocate(this%coeff(this%ncol,this%ncol))
#if defined(USE_COARRAY)
      if (this_image() == 1) then
#endif
         select case(this%method)
          case('svd')
            block
               real(rk), dimension(this%ncol,this%ncol) :: U
               real(rk), dimension(this%ncol,this%ncol) :: VT
               real(rk), dimension(this%ncol)           :: S

               call svd(cov, U,S,VT)
               this%latent = S**2/(this%ncol-1)
               this%coeff  = transpose(VT)
            end block
          case('eig')
            block
               logical,  dimension(this%ncol) :: mask
               integer,  dimension(this%ncol) :: order

               call eig(cov, this%coeff, this%latent)

               ! Sort
               mask = .true.
               do i = lbound(this%latent,1), ubound(this%latent,1)
                  order(i) = maxloc(this%latent,1,mask)
                  mask(order(i)) = .false.
               end do

               this%latent = this%latent(order)
               this%coeff  = this%coeff(:,order)
            end block
         end select

#if defined(USE_COARRAY)
      end if
#endif

   end subroutine compute_coeff