test_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 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(:),postcov(:,:),postcov_diag(:,:)
  real(dp),allocatable :: bandpostC(:,:)
  integer :: i1,i2,j1,j2, lowdiag, updiag
  real(dp) :: start,finish

  integer :: iar
  character(len=256) :: simname
  !----------------------------------
  
  !! 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)
  

  ! write arrays to output file
  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)
  
  !! 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 )


  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)
  ! call writereal1Darrh5(outfile,'nm',[nm1,nm2,nm3])
  ! call writereal1Darrh5(outfile,'nd',[nd1,nd2,nd3])

  

  allocate(postm(nm))
  postm=0.0_dp

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


  
  ! print*,'Computing postcov'
  ! i1 = 1
  ! i2 = 1000
  ! j1 = 1
  ! j2 = 1000
  ! allocate(postcov(i2-i1+1,j2-j1+1))
  ! postcov = 0.0_dp
  ! call cpu_time(start)
  ! call blockpostcov(U1,U2,U3, invlambda,iUCm1,iUCm2,iUCm3, &
  !                i1,i2,j1,j2, postcov(i1:i2,j1:j2))     
  ! call cpu_time(finish)
  ! print*,'postCov block time:',finish-start
  ! !! write posterior covariance
  ! call writereal2Darrh5(outfile,'postcov',postcov)

  
  ! print*,
  ! print*,'computing postcov_diag'
  ! !!http://www.netlib.org/lapack/lug/node124.html
  ! lowdiag = 0
  ! updiag  = 0
  ! allocate(postcov_diag(1,nm))
  ! postcov_diag = 0.0_dp
  ! call cpu_time(start)
  ! call bandpostcov(U1,U2,U3, invlambda,iUCm1,iUCm2,iUCm3, &
  !                lowdiag,updiag, postcov_diag)     
  ! call cpu_time(finish)
  ! print*,'postCov band time:',finish-start
  ! !! write posterior covariance
  ! call writereal2Darrh5(outfile,'postcov_diag',postcov_diag)


  ! 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
  ! call writereal2Darrh5(outfile,'bandpostcov',bandpostC)

  

end program test