bandpostcov()

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

30/1/2017 - OpenMPI version

Definition at line 753 of file ompi_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 :: row2(:),col1(:),recvrow(:),bandpostcmpi(:)
    integer :: nr12,nc1
    integer :: ni,nj,nk,nl,nm,nn,na,nb    

    real(dp) :: firststartt,startt
    integer,allocatable :: av(:)!,bv(:)
    integer,allocatable :: iv(:),jv(:),kv(:) !,lv(:),mv(:),nv(:) 
    integer :: p,aband,aend,bband,astart,d,bband1,bband2,ntota,everynit,mytotit
    integer,allocatable :: scheduling(:),looping(:,:)
    character(len=30) :: loopinfo
    character(len=12) :: dnum
    
    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(row2(nr12),col1(nc1))

    ! 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

    !!---------------------------------
    firststartt = mpi_wtime()
    ! !! calculate scheduling for a
    ! call spreadwork(ntota,numcpus,scheduling,looping,1)
    ! allocate(postcmpi(scheduling(myrank+1),ntotb),recvrow(ntotb),sendrow(ntotb))
    ! mytotit = looping(myrank+1,2)-looping(myrank+1,1)+1

    allocate(bandpostcmpi(size(bandpostc,2)), recvrow(size(bandpostc,2)))
    ! initialize postC
    
    bandpostc = 0.0_dp       

    do d=-lowdiag,updiag
       call para_barrier()
       
       !!--------------------------
       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

       !! calculate scheduling for a
       ntota = aend-astart+1
       call spreadwork(ntota,numcpus,scheduling,looping,astart)
       if (ntota<1000) then
          everynit = ntota/100
       else
          everynit = ntota/1000
       end if
       
       write(dnum,"(i9)") d
       loopinfo = "bandpostcov()  diag "//trim(dnum)//" "
       startt = mpi_wtime()
       if (myrank==masterrank) write(output_unit,*)
       !! indices of normal matrix
       bandpostcmpi(:) = 0.0_dp
       aband = updiag+1-d !! aband is constant within this loop
       do a=looping(myrank+1,1),looping(myrank+1,2) !!astart,aend
          !! info
          if ( (myrank==masterrank) .and. (mod(a,everynit)==0) ) then
             call timeinfo(scheduling(myrank+1),a-astart+1,startt,loopinfo)
          end if

          !if ( (myrank==masterrank) .and. (mod(a,250)==0) ) print*,"d",d,"a",a
          !!--------------------------------------------------------
          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)
          bandpostcmpi(b) = sum(row2*col1)
          !print*, a,b,aband,bband,bandpostC(aband,bband)
          !!----------------------------------------------------------------
          
       end do ! a=astart,aend
       
       !! send/recv rows                                                      
       !! collect all results in masterrank
       

       bband1 = looping(myrank+1,1)+d !!b = a+d
       bband2 = looping(myrank+1,2)+d
       
       call para_gathv1ddp(scheduling, looping(:,1)-1, masterrank, &
            bandpostcmpi(bband1:bband2), recvrow )
       if (myrank==masterrank) bandpostc(aband,:) = recvrow
      
       
    end do

Here is the caller graph for this function: