KronLinInv  0.3
Kronecker-product-based linear inversion

◆ posteriormean()

subroutine, public kronlininv::posteriormean ( real(dp), dimension(:,:), intent(in)  U1,
real(dp), dimension(:,:), intent(in)  U2,
real(dp), dimension(:,:), intent(in)  U3,
real(dp), dimension(:), intent(in)  diaginvlambda,
real(dp), dimension(:,:), intent(in)  Z1,
real(dp), dimension(:,:), intent(in)  Z2,
real(dp), dimension(:,:), intent(in)  Z3,
real(dp), dimension(:,:), intent(in)  G1,
real(dp), dimension(:,:), intent(in)  G2,
real(dp), dimension(:,:), intent(in)  G3,
real(dp), dimension(:), intent(in)  mprior,
real(dp), dimension(:), intent(in)  dobs,
real(dp), dimension(:), intent(out)  postm 
)

Computes the posterior mean

Parameters
[in]U1,U2,U3\( \mathbf{U}_1 \), \( \mathbf{U}_2 \), \( \mathbf{U}_3 \) of \( F_{\sf{A}} \)
[in]diaginvlambda\( F_{\sf{B}} \)
[in]Z1,Z2,Z3\( \mathbf{U}_1^{-1} \mathbf{C}_{\rm{M}}^{\rm{x}} (\mathbf{G}^{\rm{x}})^{\sf{T}}(\mathbf{C}_{\rm{D}}^{\rm{x}})^{-1} \), \( \mathbf{U}_2^{-1} \mathbf{C}_{\rm{M}}^{\rm{y}} (\mathbf{G}^{\rm{y}})^{\sf{T}} (\mathbf{C}_{\rm{D}}^{\rm{y}})^{-1}\), \( \mathbf{U}_3^{-1} \mathbf{C}_{\rm{M}}^{\rm{z}} (\mathbf{G}^{\rm{z}})^{\sf{T}} (\mathbf{C}_{\rm{D}}^{\rm{z}})^{-1} \) of \( F_{\sf{D}} \)
[in]G1,G2,G3The 3 forward modeling matrices \( \mathbf{G} = \mathbf{G_1} \otimes \mathbf{G_2} \otimes \mathbf{G_3} \)
[in]mpriorPrior model (vector)
[in]dobsObserved data (vector)
[out]postmCalculated posterior mean model (vector)
Date
5/8/2016 - Initial Version
27/12/2016 - New ETA calculations

Definition at line 627 of file kronlininv.f08.

627  !!
628  !! Calculate the posterior mean model
629  !!
630  real(dp),intent(in) :: u1(:,:),u2(:,:),u3(:,:), g1(:,:),g2(:,:),g3(:,:)
631  real(dp),intent(in) :: z1(:,:),z2(:,:),z3(:,:)
632  real(dp),intent(in) :: mprior(:),dobs(:),diaginvlambda(:)
633  real(dp),intent(out) :: postm(:)
634 
635  real(dp),allocatable :: row2(:),col1(:),bigmatrow(:),datrow(:)
636  real(dp),allocatable :: ddiff(:)
637  integer :: nr12,nc1,a,b,na,nb
638  integer :: ni,nj,nk,nl,nm,nn
639 
640  integer,allocatable :: av(:),bv(:)
641  integer,allocatable :: iv(:),jv(:),kv(:),lv(:),mv(:),nv(:)
642  !! ivo(:),jvo(:),kvo(:)
643  integer :: p,j,i,everynit
644 
645  real(dp),allocatable :: zh(:),eludzh(:)
646  real(dp) :: datp,elg,elrowud,tzz
647 
648  integer :: nsteps,hr,min,nthread,privcount,nthr
649  real(dp) :: eta,et1,et2,sec,startt,firststartt,endt,finish,frac
650  character(8) :: curdate
651  character(10) :: curtime
652 
653  write(output_unit,*) "posteriormean(): calculating posterior mean... [using OpenMP]"
654 
655  !! sizes
656  ni = size(z1,1)
657  nl = size(z1,2)
658  nj = size(z2,1)
659  nm = size(z2,2)
660  nk = size(z3,1)
661  nn = size(z3,2)
662 
663  !! sizes
664  ! Nr12 = size(U1,2)*size(U2,2)*size(U3,2)
665  ! Nc1 = size(Z1,1)*size(Z2,1)*size(Z3,1)
666  na = size(mprior,1)
667  nb = size(dobs,1)
668 
669  !! allocate stuff
670  allocate(ddiff(nb))
671  allocate(zh(na),eludzh(na))
672 
673  !! vectorize row and col calculations for Kron prod AxBxC
674  !!allocate(ivo(Nb),jvo(Nb),kvo(Nb))
675  allocate(av(na),iv(na),jv(na),kv(na))
676  allocate(bv(nb),lv(nb),mv(nb),nv(nb))
677  forall(p = 1:na) av(p) = p
678  forall(p = 1:nb) bv(p) = p
679 
680  !! vectors containing all possible indices for
681  !! row calculations of Kron prod AxBxC
682  iv = (av-1)/(nk*nj)+1
683  jv = (av-1-(iv-1)*nk*nj)/nk+1
684  kv = av-(jv-1)*nk-(iv-1)*nk*nj
685  !! vectors containing all possible indices for
686  !! column calculations of Kron prod AxBxC
687  lv = (bv-1)/(nn*nm) + 1
688  mv = (bv-1-(lv-1)*nn*nm)/nn + 1
689  nv = bv-(mv-1)*nn-(lv-1)*nn*nm
690  !! Gs have different shape than Us !!
691 
692  !!#######################
693  if (na<1000)then
694  everynit = na/50
695  else if (na<100000) then
696  everynit = na/100
697  else
698  everynit = na/1000
699  end if
700 
701 
702  !!---------------------------------
703  !$OMP PARALLEL PRIVATE(privcount)
704  firststartt = omp_get_wtime()
705  nthr = omp_get_num_threads()
706 
707  !!#######################
708  !!# dobs - dcalc #
709  !!#######################
710  startt= omp_get_wtime()
711  privcount = -1
712  !$OMP DO PRIVATE(b)
713  !! difference obs-calc data
714  do b=1,nb
715  privcount = privcount + 1
716  if ( (omp_get_thread_num()==0 ) .and. (mod(privcount,everynit/nthr)==0) ) then
717  frac = real(privcount,dp)/(real(nb,dp)/real(nthr,dp))
718  eta = ( (omp_get_wtime()-startt) / real(privcount,dp) ) * &
719  (real(nb-nthr*privcount,dp)/real(nthr,dp))
720  write(output_unit,fmt='(a27,f7.3,6x,a4,f12.2,1x,a3,a5)') 'posteriormean() loop 1/3: %',&
721  frac*100_dp,"ETA:",eta/60.0,"min",char(27)//'[1A'//achar(13)
722  flush(output_unit) !! to make sure it prints immediately
723  end if
724 
725  !!---------------------------------------------------------------------------
726  ddiff(b) = dobs(b) - sum(mprior * g1(lv(b),iv) * g2(mv(b),jv) * g3(nv(b),kv))
727  !!---------------------------------------------------------------------------
728 
729  end do
730  !$OMP END DO
731 
732  !!#######################
733  !!# === U d Z h === #
734  !!#######################
735  privcount = -1
736  startt= omp_get_wtime()
737  !$OMP DO PRIVATE(a)
738  do a=1,na
739  privcount = privcount + 1
740  if ( (omp_get_thread_num()==0 ) .and. (mod(privcount,everynit/nthr)==0) ) then
741  frac = real(privcount,dp)/(real(na,dp)/real(nthr,dp))
742  eta = ( (omp_get_wtime()-startt) / real(privcount,dp) ) * &
743  (real(na-nthr*privcount,dp)/real(nthr,dp))
744  write(output_unit,fmt='(a27,f7.3,6x,a4,f12.2,1x,a3,a5)') 'posteriormean() loop 2/3: %',&
745  frac*100_dp,"ETA:",eta/60.0,"min",char(27)//'[1A'//achar(13)
746  flush(output_unit) !! to make sure it prints immediately
747  end if
748 
749  !!---------------------------------------------------------------------------
750  zh(a) = sum( z1(iv(a),lv) * z2(jv(a),mv) * z3(kv(a),nv) * ddiff )
751  !!---------------------------------------------------------------------------
752 
753  end do
754  !$OMP END DO
755 
756  !!-----------------------------------------------
757  !! Second loop
758  !!### need to re-loop because full Zh is needed
759  !!#######################
760  !!# post(a) #
761  !!#######################
762  privcount = -1
763  startt= omp_get_wtime()
764  !$OMP DO PRIVATE(a)
765  do a=1,na
766  privcount = privcount + 1
767  if ( (omp_get_thread_num()==0 ) .and. (mod(privcount,everynit/nthr)==0) ) then
768  frac = real(privcount,dp)/(real(na,dp)/real(nthr,dp))
769  eta = ( (omp_get_wtime()-startt) / real(privcount,dp) ) * &
770  (real(na-nthr*privcount,dp)/real(nthr,dp))
771  write(output_unit,fmt='(a27,f7.3,6x,a4,f12.2,1x,a3,a5)') 'posteriormean() loop 3/3: %',&
772  frac*100_dp,"ETA:",eta/60.0,"min",char(27)//'[1A'//achar(13)
773  flush(output_unit) !! to make sure it prints immediately
774  end if
775 
776  !!--------------------------------------------------------------------------------
777  eludzh(a) = sum( u1(iv(a),iv) * u2(jv(a),jv) * u3(kv(a),kv) * diaginvlambda * zh )
778  !! element of the posterior mean
779  postm(a) = mprior(a) + eludzh(a)
780  !!--------------------------------------------------------------------------------
781 
782  end do
783  !$OMP END DO
784  endt = omp_get_wtime()
785  !$OMP END PARALLEL
786  write(output_unit,*)
787  print*,"posteriormean():",endt-firststartt," OMP wall clock time"
788 
realprec::dp
integer, parameter, public dp
Definition: ompi_kronlininv.f08:45
Here is the caller graph for this function: