bandpostcov()

subroutine, public kronlininv::bandpostcov	(	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)	lowdiag,
		integer, intent(in)	updiag,
		real(dp), dimension(:,:), intent(inout)	bandpostC
	)

Computes a band of the posterior covariance. See http://www.netlib.org/lapack/lug/node124.html

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]	lowdiag,updiag	Lower and upper diagonal number of requested band
[out]	postC	band of the posterior covariance stored following Lapack convention

Date

11/8/2016 - Initial Version

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

Definition at line 453 of file kronlininv.f08.

    !!
    !! Calculate a band 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(inout) :: bandpostc(:,:)
    integer,intent(in) :: lowdiag, updiag
    
    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,aband,aend,bband,astart,d
    integer :: nthr,privcount,everynit
    real(dp) :: eta,frac,startt,firststartt,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    

    if (na /= nb) then
       write(*,*) '(Na /= Nb)', na,nb
       stop
    end if
    if ( (updiag>=na) .or. (lowdiag>=na) .or. (lowdiag<0) .or. (updiag<0) ) then
       write(*,*) .or."(updiag<Na)  (lowdiag<Na)"
       write(*,*) "updiag",updiag,"Na",na,"lowdiag",lowdiag,"Na",na
       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 
    !! vectors containing all possible indices for
    !!    column calculations of Kron prod AxBxC
    ! lv =  (bv-1)/(Nn*Nm) + 1 
    ! mv =  (bv-1-(lv-1)*Nn*Nm)/Nn + 1 
    ! nv =  bv-(mv-1)*Nn-(lv-1)*Nn*Nm
    
    !! 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))
   
    
    everynit = 250
    

    ! Lapack: http://www.netlib.org/lapack/lug/node124.html
    ! aij is stored in AB(ku+1+i-j,j) for max(1,j-ku) <= i <= \min(m,j+kl).
    !------------
    ! Diagonals of a matrix
    ! i + d = j
    ! main diag d = 0
    ! upper d > 0
    ! lower d < 0
    ! diagonals of a matrix and indices of related band matrix
    ! ONLY for square matrix
   
    ! initialize postC
    bandpostc = 0.0_dp

    !!=====================================
    firststartt = omp_get_wtime()
    
    do d=-lowdiag,updiag
       if (d<0) then
          astart = abs(d)+1
       else
          astart = 1
       end if
       if (d>0) then
          aend = na-abs(d)
       else
          aend = na
       end if
       !print*,'diagonal',d
       !! indices of normal matrix
       
       print*,"Diagonal ",d," from range [",-lowdiag,",",updiag,"]"
       !-------------------------
       !$OMP PARALLEL
       privcount = -1
       startt = omp_get_wtime()
       nthr = omp_get_num_threads()
       
       !$OMP DO PRIVATE(row2,col1,a,b,aband,bband,privcount)
       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(nb,dp)/real(nthr,dp))
             eta = ( (omp_get_wtime()-startt) / real(privcount,dp) ) * &
                  (real(aend-astart+1-nthr*privcount,dp)/real(nthr,dp))
             write(output_unit,fmt='(a27,f7.3,6x,a4,f12.2,1x,a3,a5)') 'bandpostcov():  %',&
                  frac*100_dp,"ETA:",eta/60.0,"min",char(27)//'[1A'//achar(13)
             flush(output_unit) !! to make sure it prints immediately          
          end if

          b = a+d
          !! indices of the band matrix
          aband = updiag+1+a-b
          bband = b

          !! row first two factors
          row2 = diaginvlambda * u1(iv(a),iv) * u2(jv(a),jv) * u3(kv(a),kv)

          !! 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
          !! store it in the band storage format
          bandpostc(aband,bband) = sum(row2*col1)
          !print*, a,b,aband,bband,bandpostC(aband,bband)

       end do ! a=astart,aend
       !$OMP END DO
       endt = omp_get_wtime()
       !$OMP END PARALLEL
       write(output_unit,*)
       print*,"bandpostcov():",endt-firststartt," OMP wall clock time"
    
    end do