blockpostcov()

subroutine, public kronlininv::blockpostcov	(	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)	iUCm1,
		real(dp), dimension(:,:), intent(in)	iUCm2,
		real(dp), dimension(:,:), intent(in)	iUCm3,
		integer, intent(in)	astart,
		integer, intent(in)	aend,
		integer, intent(in)	bstart,
		integer, intent(in)	bend,
		real(dp), dimension(:,:), intent(out)	postC
	)

Computes a block of the posterior covariance.

Parameters

[in]	U1,U2,U3	\( \mathbf{U}_1 \), \( \mathbf{U}_2 \), \( \mathbf{U}_3 \) of \( F_{\sf{A}} \)
[in]	diaginvlambda	\( F_{\sf{B}} \)
[in]	iUCm1,iUCm2,iUCm3	\(\mathbf{U}_1^{-1} \mathbf{C}_{\rm{M}}^{\rm{x}}\), \(\mathbf{U}_2^{-1} \mathbf{C}_{\rm{M}}^{\rm{y}}\), \(\mathbf{U}_2^{-1} \mathbf{C}_{\rm{M}}^{\rm{z}}\) of \( F_{\sf{C}} \)
[in]	astart,aend	index of first/last row of the requested block
[in]	bstart,bend	index of first/last column of the requested block
[out]	postC	block of the posterior covariance

Date

5/8/2016 - Initial Version

4/11/2016 - Removed indices lv,mv,kv

27/13/2016 - New ETA calculations

Definition at line 314 of file kronlininv.f08.

    !!
    !! Calculates a block of the posterior covariance
    !! 
    real(dp),intent(in) :: u1(:,:),u2(:,:),u3(:,:),iucm1(:,:), &
         iucm2(:,:),iucm3(:,:)
    !! diaginvlambda = (I + lam1 x lam2 x lam3 )^-1
    real(dp),intent(in) :: diaginvlambda(:) !! diagonal/vector central factor
    real(dp),intent(out) :: postc(:,:)
    integer,intent(in) :: astart,aend,bstart,bend
    
    integer :: a,b
    real(dp),allocatable :: row1(:),row2(:),col1(:)
    integer :: nr12,nc1
    integer :: ni,nj,nk,nl,nm,nn,na,nb    

    integer,allocatable :: av(:)!,bv(:)
    integer,allocatable :: iv(:),jv(:),kv(:) !,lv(:),mv(:),nv(:) 
    integer :: p,csha1,csha2,everynit,nthr,privcount,totcount 
    real(dp) :: eta,frac,startt,endt
    
    ni = size(u1,1)
    nl = size(u1,2)
    nj = size(u2,1)
    nm = size(u2,2)
    nk = size(u3,1)
    nn = size(u3,2)
    na = ni*nj*nk
    nb = nl*nm*nn
    
    !! check shape of output array
    if (na /= nb) then
       write(*,*) '(Na /= Nb)', na,nb
       stop
    end if
    csha1 = aend-astart+1
    csha2 = bend-bstart+1
    if ( (size(postc,1)/=csha1) .or. (size(postc,2)/=csha2)) then
       write(*,*) "Wrong size of the posterior covariance block array."
       stop
    end if
    !! check limits of requested block
    if ( (astart<1) .or. (aend>na) .or. (astart>aend) .or. &
         (bstart<1) .or. (bend>na) .or. (bstart>bend) ) then
       write(*,*) "Wrong size of the requested block array."
       stop
    end if
        
    !! vectorize row and col calculations for Kron prod AxBxC
    allocate(av(na),iv(na),jv(na),kv(na))
    !allocate(bv(Nb),lv(Nb),mv(Nb),nv(Nb))
    forall(p = 1:na) av(p) = p
    !forall(p = 1:Nb) bv(p) = p

    !! vectors containing all possible indices for 
    !!    row calculations of Kron prod AxBxC
    iv =  (av-1)/(nk*nj)+1 
    jv =  (av-1-(iv-1)*nk*nj)/nk + 1 
    kv =  av-(jv-1)*nk-(iv-1)*nk*nj 
     
    !! allocate stuff
    nr12 = size(u1,2)* size(u2,2)* size(u3,2)
    nc1  = size(iucm1,1)* size(iucm2,1)* size(iucm3,1)
    allocate(row1(nr12),row2(nr12),col1(nc1))

    !!------------------------
    if (aend-astart+1<20) then
       everynit = 1
    else if (aend-astart+1<100000) then
       everynit = (aend-astart+1)/100
    else
       everynit = (aend-astart+1)/1000
    end if

    
    !!------------------------
    !$OMP PARALLEL PRIVATE(privcount)
    startt = omp_get_wtime()
    nthr = omp_get_num_threads()

    totcount = aend-astart+1
    privcount = -1
     if ( omp_get_thread_num()==0 ) write(output_unit,*)
    !$OMP DO PRIVATE(row2,col1,a,b)
    do a=astart,aend
       privcount = privcount + 1 
       if ( (omp_get_thread_num()==0 ) .and. (mod(privcount,everynit/nthr)==0) ) then          
          frac = real(privcount,dp)/(real(totcount,dp)/real(nthr,dp))
          eta = ( (omp_get_wtime()-startt) / real(privcount,dp) ) * &
               (real(totcount-nthr*privcount,dp)/real(nthr,dp))
          write(output_unit,fmt='(a19,f7.3,6x,a4,f12.2,1x,a3,a5)') 'blockpostcov():  %',&
               frac*100_dp,"ETA:",eta/60.0,"min",char(27)//'[1A'//achar(13)
          flush(output_unit) !! to make sure it prints immediately          
       end if

       !! row first two factors
       !! a row x diag matrix 
       !!row2 =  U1(iv(a),lv) * U2(jv(a),mv) * U3(kv(a),nv) * diaginvlambda
       row2 =  u1(iv(a),iv) * u2(jv(a),jv) * u3(kv(a),kv) * diaginvlambda
       
       do b=bstart,bend
          
          !! calculate one row of first TWO factors
          !!call columnAxBxC(Ni,Nj,Nk,Nl,Nm,Nn, iUCm1,iUCm2,iUCm3,b,col1)
          col1 = iucm1(iv,iv(b)) * iucm2(jv,jv(b)) * iucm3(kv,kv(b))

          !! calculate one element of the posterior covariance
          postc(a,b) = sum(row2*col1)

       end do       
    end do
    !$OMP END DO
    endt = omp_get_wtime()
    !$OMP END PARALLEL
    write(output_unit,*)
    print*,"blockpostcov():",endt-startt," OMP wall clock time"