blockpostcov()

subroutine, public ompi_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

27/1/2017 - OpenMPI version

Definition at line 585 of file ompi_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,irow,i,j,sender
    real(dp),allocatable :: row2(:),col1(:)
    integer :: nr12,nc1
    integer :: ni,nj,nk,nl,nm,nn,na,nb,ntota,ntotb   

    integer,allocatable :: av(:)!,bv(:)
    integer,allocatable :: iv(:),jv(:),kv(:) !,lv(:),mv(:),nv(:) 
    integer :: p,everynit,impicount,mytotit,rirow
    real(dp) :: startt,firststartt,endt
    integer,allocatable :: scheduling(:),looping(:,:)
    real(dp),allocatable :: postcmpi(:,:),recvrow(:),sendrow(:)
    character(len=30) :: loopinfo

    if (myrank==masterrank) then
       write(output_unit,*) "posteriormean(): calculating posterior mean... [using OpenMPI]"
    end if
    
    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
    ntota = aend-astart+1
    ntotb = bend-bstart+1
    if ( (size(postc,1)/=ntota) .or. (size(postc,2)/=ntotb)) 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(row2(nr12),col1(nc1))

    !!------------------------
    if (ntota<1000) then
       everynit = 20
    else 
       everynit = ntota/1000
    end if
        
    !!---------------------------------
    firststartt = mpi_wtime()
    loopinfo = "blockpostcov()  "
    if (myrank==masterrank) write(output_unit,*)
    !! calculate scheduling for a
    call spreadwork(ntota,numcpus,scheduling,looping,1)
    allocate(postcmpi(scheduling(myrank+1),ntotb),recvrow(ntotb),sendrow(ntotb))
    mytotit = scheduling(myrank+1)

    !!print*,myrank,":",looping(myrank+1,1),"|",looping(myrank+1,2),"|",scheduling
    
    !!-----------------------
    startt = mpi_wtime()
    impicount = 0
    do irow=looping(myrank+1,1),looping(myrank+1,2)
       !! set a
       a=irow+astart-1 
       !! info
       if ( (myrank==masterrank) .and. (mod(irow,everynit)==0) ) then
          call timeinfo(mytotit,irow,startt,loopinfo)
       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

       j=0
       do b=bstart,bend
          j=j+1
          !! 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
          i = irow-looping(myrank+1,1)+1
          postcmpi(i,j) = sum(row2*col1)

       end do     
     
    end do

    !! collect all results in masterrank
    call para_barrier()
    do sender=0,numcpus-1   
       if ( (myrank==sender) .or. (myrank==masterrank) ) then
             do i=looping(sender+1,1),looping(sender+1,2)
             !!----------------------
             !! send/rec each row
             !print*,"==>",myrank,": sender=",sender,"i=",i
             if (myrank==sender) sendrow = postcmpi(i,:)
             call para_srarr1ddp(sender,masterrank,i,sendrow,recvrow)
             !call para_srint(myrank,masterrank,impicount,irow,rirow)
             if (myrank==masterrank)  postc(i,:) = recvrow
          end do
       end if
    end do
    
    call para_barrier()
    if (myrank==masterrank) then
       write(output_unit,*)
       endt = MPI_Wtime()
       print*,"blockpostcov():",endt-firststartt," MPI wall clock time"
    end if    

Here is the call graph for this function:

Here is the caller graph for this function: