test_ompi_kronlininv.f08

Go to the documentation of this file.

!------------------------------------------------------------------------
!
!    Copyright 2018, Andrea Zunino 
!
!    This file is part of KronLinInv.
!
!    KronLinInv 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.
!
!    KronLinInv 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 KronLinInv.  If not, see <http://www.gnu.org/licenses/>.
!
!------------------------------------------------------------------------

!-----------------------------------------------------
!
!
!
!-----------------------------------------------------
program test

  use realprec
  use parautil
  use ompi_kronlininv
  use readwriteh5

  implicit none

  real(dp),allocatable :: G1(:,:),G2(:,:),G3(:,:),&
         Cm1(:,:),Cm2(:,:),Cm3(:,:),&
         Cd1(:,:),Cd2(:,:),Cd3(:,:),&
         mprior(:),dobs(:),mtrue(:)
  real(dp),allocatable :: U1(:,:),U2(:,:),U3(:,:),&
         iUCm1(:,:),iUCm2(:,:),iUCm3(:,:),&
         iUCmGtiCd1(:,:),iUCmGtiCd2(:,:),iUCmGtiCd3(:,:),&
         invlambda(:)  
  character(len=1024) :: inpfile,outfile
  integer :: nm,nd,nm1,nm2,nm3,nd1,nd2,nd3
  real(dp),allocatable :: postm(:),postmFAST(:),postcov(:,:)
  real(dp),allocatable :: bandpostC(:,:)
  integer :: i1,i2,j1,j2, lowdiag, updiag
  real(dp) :: start,finish

  integer :: iar,mpierr
  character(len=256) :: simname
  !----------------------------------


  !! Initialize OpenMPI 
  call para_init()   

  !! simname -> name of the input file (without the .hdf5)
  
  iar = 1 
  call get_command_argument(iar, simname)

  inpfile = trim(simname)//'.h5'
  outfile = trim(simname)//'_output.h5'

  print*,"inpfile:",trim(inpfile)
  print*,"outfile:",trim(outfile)

  !----------------------------------

  !! read all input arrays
  print*,'reading input arrays'
  call readreal2darrh5(inpfile,'G1',g1)
  call readreal2darrh5(inpfile,'G2',g2)
  call readreal2darrh5(inpfile,'G3',g3)
  call readreal2darrh5(inpfile,'Cm1',cm1)
  call readreal2darrh5(inpfile,'Cm2',cm2)
  call readreal2darrh5(inpfile,'Cm3',cm3)
  call readreal2darrh5(inpfile,'Cd1',cd1)
  call readreal2darrh5(inpfile,'Cd2',cd2)
  call readreal2darrh5(inpfile,'Cd3',cd3)
  call readreal1darrh5(inpfile,'mtrue',mtrue)
  call readreal1darrh5(inpfile,'dobs',dobs)
  call readreal1darrh5(inpfile,'mprior',mprior)

  call para_barrier()
  
  if (myrank==masterrank) then
     !! write arrays to output file
     print*,"WRITE hello from ",myrank,masterrank 
     call writereal2darrh5(outfile,'G1',g1)
     call writereal2darrh5(outfile,'G2',g2)
     call writereal2darrh5(outfile,'G3',g3)
     call writereal2darrh5(outfile,'Cm1',cm1)
     call writereal2darrh5(outfile,'Cm2',cm2)
     call writereal2darrh5(outfile,'Cm3',cm3)
     call writereal2darrh5(outfile,'Cd1',cd1)
     call writereal2darrh5(outfile,'Cd2',cd2)
     call writereal2darrh5(outfile,'Cd3',cd3)
     call writereal1darrh5(outfile,'mtrue',mtrue)
     call writereal1darrh5(outfile,'dobs',dobs)
     call writereal1darrh5(outfile,'mprior',mprior)
  end if


  !! sizes
  nm1 = size(cm1,1)
  nm2 = size(cm2,1)
  nm3 = size(cm3,1)
  nd1 = size(cd1,1)
  nd2 = size(cd2,1)
  nd3 = size(cd3,1)
  nm = nm1*nm2*nm3
  nd = nd1*nd2*nd3

  print*,"Sizes:"
  print*,"M:",nm1,nm2,nm3,nm
  print*,"D:",nd1,nd2,nd3,nd
  print*

  !! allocate working arrays
  print*,'Allocating working arrays'
  allocate(u1(nm1,nm1),u2(nm2,nm2),u3(nm3,nm3))
  allocate(invlambda(nm))
  allocate(iucm1(nm1,nm1),iucm2(nm2,nm2),iucm3(nm3,nm3))
  allocate(iucmgticd1(nm1,nd1),iucmgticd2(nm2,nd2),iucmgticd3(nm3,nd3))


  !! compute the factors for post mean and covariance
  print*,'Computing the factors for post mean and covariance'
  call calcfactors(g1,g2,g3,cm1,cm2,cm3,cd1,cd2,cd3,&
       u1,u2,u3,invlambda,iucm1,iucm2,iucm3,&
       iucmgticd1,iucmgticd2,iucmgticd3 )
  print*,'END computing the factors for post mean and covariance'
  
    
  call para_barrier()
  if (myrank==masterrank) then 
     call writereal2darrh5(outfile,'U1',u1)
     call writereal2darrh5(outfile,'U2',u2)
     call writereal2darrh5(outfile,'U3',u3)
     call writereal1darrh5(outfile,'invlambda',invlambda)
     call writereal2darrh5(outfile,'iUCm1',iucm1)
     call writereal2darrh5(outfile,'iUCm2',iucm2)
     call writereal2darrh5(outfile,'iUCm3',iucm3)
     call writereal2darrh5(outfile,'iUCmGtiCd1',iucmgticd1)
     call writereal2darrh5(outfile,'iUCmGtiCd2',iucmgticd2)
     call writereal2darrh5(outfile,'iUCmGtiCd3',iucmgticd3)
  end if



    
  !!=====================================================================
  call para_barrier()
  allocate(postm(nm)) 
  postm=0.0_dp

  if (myrank==masterrank) then
     call cpu_time(start)
     print*,'Computing post mean'
  end if
  
  call posteriormean(u1,u2,u3, invlambda, iucmgticd1,iucmgticd2,iucmgticd3,&
       g1,g2,g3,mprior,dobs, postm ) 
  
  if (myrank==masterrank) then
     call cpu_time(finish)
     print*,'postm time:',finish-start
     !! write posterior mean
     call writereal1darrh5(outfile,'postm',postm)
  end if
  !!=====================================================================
  



 
  !!=====================================================================
  call para_barrier()
  allocate(postcov(nm,nm))
  postcov = 0.0_dp
  i1 = 10
  i2 = 1000
  j1 = 15
  j2 = 1000
  if (myrank==masterrank) then
     call cpu_time(start)
     print*,'Computing block postcov'
  end if
     
  call blockpostcov(u1,u2,u3, invlambda,iucm1,iucm2,iucm3, &
                 i1,i2,j1,j2, postcov(i1:i2,j1:j2))     

  if (myrank==masterrank) then
     call cpu_time(finish)
     print*,'postcov block time:',finish-start
     !! write posterior covariance
     call writereal2darrh5(outfile,'postcov',postcov)
  end if  
  !!=====================================================================

  



  !!=====================================================================
  call para_barrier()
  if (myrank==masterrank) print*,'Computing band of post cov'
  lowdiag = 10
  updiag  = 10
  !! allocate the array for band storage
  allocate(bandpostc(lowdiag+updiag+1,nm))
  call bandpostcov(u1,u2,u3, invlambda,iucm1,iucm2,iucm3, &
                  lowdiag, updiag, bandpostc)     

  !! write posterior covariance
  if (myrank==masterrank) then
     call writereal2darrh5(outfile,'bandpostcov',bandpostc)
  end if
  !!=====================================================================


  
  call para_finish()
  !stop
  !!-------------------------------------------------------------

end program