generate_pf.f90

Go to the documentation of this file.

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!! Time-stamp: <2016-10-18 15:37:32 pbrowne>
!!!
!!!    Subroutine to generate Pf matrix given ensemble members on a communicator
!!!    Copyright (C) 2014  Philip A. Browne
!!!
!!!    This program is free software: you can redistribute it and/or modify
!!!    it under the terms of the GNU General Public License as published by
!!!    the Free Software Foundation, either version 3 of the License, or
!!!    (at your option) any later version.
!!!
!!!    This program is distributed in the hope that it will be useful,
!!!    but WITHOUT ANY WARRANTY; without even the implied warranty of
!!!    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!!!    GNU General Public License for more details.
!!!
!!!    You should have received a copy of the GNU General Public License
!!!    along with this program.  If not, see <http://www.gnu.org/licenses/>.
!!!
!!!    Email: p.browne @ reading.ac.uk
!!!    Mail:  School of Mathematical and Physical Sciences,
!!!           University of Reading,
!!!       Reading, UK
!!!       RG6 6BB
!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


subroutine generate_pf(stateDim,cnt,comm,x,pf)
  use output_empire, only : emp_e
  implicit none
  include 'mpif.h'
  integer, parameter :: rk = kind(1.0d0)
  integer, intent(in) :: stateDim
  integer, intent(in) :: cnt
  integer, intent(in) :: comm
  real(kind=rk), dimension(stateDim,cnt), intent(in) :: x
  real(kind=rk), dimension(stateDim*(stateDim+1)/2), intent(out) :: Pf

  
  
  real(kind=rk), dimension(stateDim) :: mean
  real(kind=rk), dimension(stateDim,cnt) :: xp

  integer :: m
  integer :: i ! counter
  integer :: mpi_err !error flag for mpi

  ! get the total number of ensemble members used:
  call mpi_allreduce(cnt,m,1,mpi_integer,mpi_sum,comm,mpi_err)
  if(mpi_err .ne. mpi_success) then
     write(emp_e,*) 'ERROR in generate_pf: mpi_allreduce 1 failed with flag ',mpi_err
     stop
  end if
  
  !get the ensemble mean
  xp(:,1) = sum(x,dim=2)
  call mpi_allreduce(xp(:,1),mean,statedim,mpi_double_precision&
       &,mpi_sum,comm,mpi_err)
  if(mpi_err .ne. mpi_success) then
     write(emp_e,*) 'ERROR in generate_pf: mpi_allreduce 2 failed with flag ',mpi_err
     stop
  end if
  ! use BLAS to perform mean = mean/real(m,rk)
  call dscal(statedim,1.0d0/real(m,rk),mean,1)


  !compute the unscaled ensemble perturbation matrix on this process
  do i = 1,cnt
     xp(:,i) = x(:,i)-mean
  end do

  !form pf with only ensemble members on this process
  call dsfrk('N',&      !TRANSR The normal form of RFP A is stored
             'U',&      !UPLO   upper part of matrix formed
             'N',&      !TRANS  no transposed done
             statedim,& !N      dim of output matrix [square]
             cnt,&      !K      number of columns in input
             1.0d0,&    !ALPHA  alpha scalar
             xp,&       !A      input matrix
             statedim,& !LDA    leading dimension of A
             0.0d0,&    !BETA   add to empty array
             pf,&       !C      output matrix
             statedim&  !LDC    dimension of output matrix
             )

  ! get pf on all processors
  call mpi_allreduce(mpi_in_place,pf,statedim*(statedim+1)/2,mpi_double_precision&
       &,mpi_sum,comm,mpi_err)
  if(mpi_err .ne. mpi_success) then
     write(emp_e,*) 'ERROR in generate_pf: mpi_allreduce 3 failed with flag ',mpi_err
     stop
  end if
  
  !scale pf appropriately
  ! use BLAS to perform pf = pf/(m-1)
  call dscal(statedim*(statedim+1)/2,1.0d0/real(m-1,rk),pf,1)
  
end subroutine generate_pf