kronlininv.f08

Go to the documentation of this file.

!------------------------------------------------------------------------
!
!    Copyright 2018, Andrea Zunino 
!
!    This file is part of KronLinInv.
!
!    KronLinInv is free software: you can redistribute it and/or modify
!    it under the terms of the GNU General Public License as published by
!    the Free Software Foundation, either version 3 of the License, or
!    (at your option) any later version.
!
!    KronLinInv is distributed in the hope that it will be useful,
!    but WITHOUT ANY WARRANTY; without even the implied warranty of
!    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!    GNU General Public License for more details.
!
!    You should have received a copy of the GNU General Public License
!    along with KronLinInv.  If not, see <http://www.gnu.org/licenses/>.
!
!------------------------------------------------------------------------





!-------------------------------------------------------------
!
!-------------------------------------------------------------
module realprec

  ! http://fortranwiki.org/fortran/show/Real+precision
  integer,parameter  :: pdigits=15  !15 !! used also by MPI
  integer,parameter  :: rexprange=307  !307 !! used also by MPI
    
  ! RESULT = SELECTED_REAL_KIND([P, R, RADIX]) 
  integer,parameter  :: dp = selected_real_kind(pdigits, rexprange)

  ! integer, parameter :: sp = selected_real_kind(6, 37)
  ! integer, parameter :: dp = selected_real_kind(15, 307)
  ! integer, parameter :: qp = selected_real_kind(33, 4931)

  private :: pdigits,rexprange
  public  :: dp

end module realprec


!################################################################
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!################################################################

!-------------------------------------------------------------------------
!
!
!
!-------------------------------------------------------------------------
module kronlininv
  
  use realprec
  use, intrinsic :: iso_fortran_env, only: output_unit

  implicit none
!#ifdef omppara  
  include 'omp_lib.h'
!#endif

  
  private :: solvels,symgeneigvv,symsolvels
  public :: calcfactors,blockpostcov,posteriormean,bandpostcov
  
contains

  !!!=====================================================
  !--------------------------------------------------
  ! DESCRIPTION: 
  !  REVISION HISTORY:
  !
  !--------------------------------------------------
  subroutine calcfactors(G1,G2,G3,Cm1,Cm2,Cm3,Cd1,Cd2,Cd3,&
       U1,U2,U3,diaginvlambda,iUCm1,iUCm2,iUCm3,&
       iUCmGtiCd1,iUCmGtiCd2,iUCmGtiCd3 ) 
 
    real(dp),intent(in) :: G1(:,:),G2(:,:),G3(:,:),&
         Cm1(:,:),Cm2(:,:),Cm3(:,:),&
         Cd1(:,:),Cd2(:,:),Cd3(:,:)
    
    real(dp),intent(out) :: U1(:,:),U2(:,:),U3(:,:),&
         iUCm1(:,:),iUCm2(:,:),iUCm3(:,:),&
         iUCmGtiCd1(:,:),iUCmGtiCd2(:,:),iUCmGtiCd3(:,:),&
         diaginvlambda(:)
    
    real(dp),allocatable :: lambda1(:),lambda2(:),lambda3(:)
    real(dp),allocatable :: iCdG1(:,:),iCdG2(:,:),iCdG3(:,:),&
         GtiCdG1(:,:),GtiCdG2(:,:),GtiCdG3(:,:),&
         ide1(:,:),ide2(:,:),ide3(:,:),&
         iCd1(:,:),iCd2(:,:),iCd3(:,:)
    integer :: i,j,k,p,nm,nm1,nm2,nm3,nd1,nd2,nd3
    character :: uplo
   
    !! Lapack routines
    external dsygv

    nm1 = size(cm1,1)
    nm2 = size(cm2,1)
    nm3 = size(cm3,1)
    nd1 = size(cd1,1)
    nd2 = size(cd2,1)
    nd3 = size(cd3,1)
    nm = nm1*nm2*nm3

    !! compute preliminary stuff
    print*,'calcfactors(): compute preliminary stuff'
    allocate(icdg1(nd1,nm1),icdg2(nd2,nm2),icdg3(nd3,nm3))
    call symsolvels( cd1, g1, icdg1 )
    call symsolvels( cd2, g2, icdg2 )
    call symsolvels( cd3, g3, icdg3 )

    !allocate(GtiCdG1(nm1,nm1),GtiCdG2(nm2,nm2),GtiCdG3(nm3,nm3))
    gticdg1  = matmul( transpose(g1), icdg1 )
    gticdg2  = matmul( transpose(g2), icdg2 ) 
    gticdg3  = matmul( transpose(g3), icdg3 )

    !!--------------------------
    !!          fa
    !!--------------------------
    !! compute eigendecomposition
    print*,'calcfactors(): compute fa'
    allocate(lambda1(nm1),lambda2(nm2),lambda3(nm3))
    uplo = 'L'
    call symgeneigvv(gticdg1,uplo,cm1,lambda1,u1)
    call symgeneigvv(gticdg2,uplo,cm2,lambda2,u2)
    call symgeneigvv(gticdg3,uplo,cm3,lambda3,u3)

    !!--------------------------
    !!          fb
    !!--------------------------
    !! calculates diagonal central factor
    print*,'calcfactors():compute fb'
    p=1
    do i=1,nm1
       do j=1,nm2
          do k=1,nm3
             diaginvlambda(p) = 1.0_dp/(1.0_dp+lambda1(i)*lambda2(j)*lambda3(k))
             p = p+1
          end do 
       end do
    end do

    !!--------------------------
    !!          fc
    !!--------------------------
    !! computes U1^-1 Cmx, U2^-1 Cmy, U3^-1 Cmz
    print*,'calcfactors(): compute fc'
    call solvels( u1, cm1, iucm1 )
    call solvels( u2, cm2, iucm2 )
    call solvels( u3, cm3, iucm3 )
    
    !!--------------------------
    !!          fd
    !!--------------------------
    print*,'calcfactors(): compute fd'
    allocate(ide1(nd1,nd1),ide2(nd2,nd2),ide3(nd3,nd3))
    ide1 = 0.0_dp  ! Initialize the array.
    ide2 = 0.0_dp
    ide3 = 0.0_dp
    forall(i = 1:nd1) ide1(i,i) = 1.0_dp  ! Set the diagonal.
    forall(i = 1:nd2) ide2(i,i) = 1.0_dp 
    forall(i = 1:nd3) ide3(i,i) = 1.0_dp
    allocate(icd1(nd1,nd1),icd2(nd2,nd2),icd3(nd3,nd3))
    call symsolvels( cd1, ide1, icd1 )
    call symsolvels( cd2, ide2, icd2 )
    call symsolvels( cd3, ide3, icd3 )

    iucmgticd1 = matmul( iucm1, matmul( transpose(g1), icd1 ))     
    iucmgticd2 = matmul( iucm2, matmul( transpose(g2), icd2 )) 
    iucmgticd3 = matmul( iucm3, matmul( transpose(g3), icd3 ))

    ! print*,
    ! print*,shape(U1),shape(U2),shape(U3)
    ! print*,shape(iUCm1),shape(iUCm2),shape(iUCm3)
    ! print*,shape(iUCmGtiCd1),shape(iUCmGtiCd2),shape(iUCmGtiCd3)
    ! print*,
    print*,'calcfactors(): end '
  end subroutine calcfactors

  !!!==================================================
  !-------------------------------------------------
  !
  !
  !
  ! REVISION HISTORY:
  !
  !-------------------------------------------------  
  subroutine blockpostcov(U1,U2,U3, diaginvlambda, &
       iUCm1,iUCm2,iUCm3, astart,aend,bstart,bend, postC) 
    !!
    !! 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"

  end subroutine blockpostcov

  !!!==================================================
  !-------------------------------------------------
  !
  !
  ! REVISION HISTORY:
  !
  !-------------------------------------------------  
  subroutine bandpostcov(U1,U2,U3, diaginvlambda, &
       iUCm1,iUCm2,iUCm3, lowdiag, updiag, bandpostC) 
    !!
    !! 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
 
    
  end subroutine bandpostcov

  !!!===============================================================
  !-------------------------------------------------
  !
  !
  !
  ! REVISION HISTORY:
  ! 
  !-------------------------------------------------
  subroutine posteriormean(U1,U2,U3, diaginvlambda, Z1,Z2,Z3,&
       G1,G2,G3, mprior, dobs, postm) 
    !!
    !! Calculate the posterior mean model
    !! 
    real(dp),intent(in)  :: U1(:,:),U2(:,:),U3(:,:), G1(:,:),G2(:,:),G3(:,:)
    real(dp),intent(in)  :: Z1(:,:),Z2(:,:),Z3(:,:)
    real(dp),intent(in)  :: mprior(:),dobs(:),diaginvlambda(:)
    real(dp),intent(out) :: postm(:)

    real(dp),allocatable :: row2(:),col1(:),bigmatrow(:),datrow(:)
    real(dp),allocatable :: ddiff(:)
    integer :: Nr12,Nc1,a,b,Na,Nb
    integer :: Ni,Nj,Nk,Nl,Nm,Nn

    integer,allocatable :: av(:),bv(:)
    integer,allocatable :: iv(:),jv(:),kv(:),lv(:),mv(:),nv(:)
        !! ivo(:),jvo(:),kvo(:)
    integer :: p,j,i,everynit

    real(dp),allocatable :: Zh(:),elUDZh(:)
    real(dp) :: datp,elg,elrowud,tZZ

    integer :: nsteps,hr,min,nthread,privcount,nthr
    real(dp) :: eta,et1,et2,sec,startt,firststartt,endt,finish,frac
    character(8)  :: curdate
    character(10) :: curtime

    write(output_unit,*) "posteriormean(): calculating posterior mean... [using OpenMP]"
    
    !! sizes
    ni = size(z1,1)
    nl = size(z1,2)
    nj = size(z2,1)
    nm = size(z2,2)
    nk = size(z3,1)
    nn = size(z3,2)

    !! sizes
    ! Nr12 = size(U1,2)*size(U2,2)*size(U3,2)
    ! Nc1  = size(Z1,1)*size(Z2,1)*size(Z3,1)
    na   = size(mprior,1)
    nb   = size(dobs,1)
    
    !! allocate stuff
    allocate(ddiff(nb))
    allocate(zh(na),eludzh(na))
 
    !! vectorize row and col calculations for Kron prod AxBxC
    !!allocate(ivo(Nb),jvo(Nb),kvo(Nb))
    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 
    !!  Gs have different shape than Us !!

    !!#######################
    if (na<1000)then
       everynit = na/50
    else if (na<100000) then
       everynit = na/100
    else
       everynit = na/1000
    end if

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

    !!#######################
    !!#    dobs - dcalc     #
    !!#######################
    startt= omp_get_wtime()
    privcount = -1
    !$OMP DO PRIVATE(b)
    !! difference obs-calc data
    do b=1,nb
       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(nb-nthr*privcount,dp)/real(nthr,dp))
          write(output_unit,fmt='(a27,f7.3,6x,a4,f12.2,1x,a3,a5)') 'posteriormean() loop 1/3: %',&
               frac*100_dp,"ETA:",eta/60.0,"min",char(27)//'[1A'//achar(13)
          flush(output_unit) !! to make sure it prints immediately          
       end if
       
       !!---------------------------------------------------------------------------
       ddiff(b) = dobs(b) - sum(mprior * g1(lv(b),iv) * g2(mv(b),jv) * g3(nv(b),kv))
       !!---------------------------------------------------------------------------

    end do
    !$OMP END DO
       
    !!#######################
    !!#   === U d Z h ===   #
    !!#######################
    privcount = -1
    startt= omp_get_wtime()
    !$OMP DO PRIVATE(a)
    do a=1,na
       privcount = privcount + 1 
       if ( (omp_get_thread_num()==0 ) .and. (mod(privcount,everynit/nthr)==0) ) then          
          frac = real(privcount,dp)/(real(na,dp)/real(nthr,dp))
          eta = ( (omp_get_wtime()-startt) / real(privcount,dp) ) * &
               (real(na-nthr*privcount,dp)/real(nthr,dp))
          write(output_unit,fmt='(a27,f7.3,6x,a4,f12.2,1x,a3,a5)') 'posteriormean() loop 2/3: %',&
               frac*100_dp,"ETA:",eta/60.0,"min",char(27)//'[1A'//achar(13)
          flush(output_unit) !! to make sure it prints immediately          
       end if

       !!---------------------------------------------------------------------------
       zh(a) = sum( z1(iv(a),lv) * z2(jv(a),mv) * z3(kv(a),nv) * ddiff )
       !!---------------------------------------------------------------------------
       
    end do
    !$OMP END DO
        
    !!-----------------------------------------------
    !! Second loop
    !!### need to re-loop because full Zh is needed
    !!#######################
    !!#       post(a)       #
    !!#######################
    privcount = -1
    startt= omp_get_wtime()
    !$OMP DO PRIVATE(a)
    do a=1,na
       privcount = privcount + 1
       if ( (omp_get_thread_num()==0 ) .and. (mod(privcount,everynit/nthr)==0) ) then          
          frac = real(privcount,dp)/(real(na,dp)/real(nthr,dp))
          eta = ( (omp_get_wtime()-startt) / real(privcount,dp) ) * &
               (real(na-nthr*privcount,dp)/real(nthr,dp))
          write(output_unit,fmt='(a27,f7.3,6x,a4,f12.2,1x,a3,a5)') 'posteriormean() loop 3/3: %',&
               frac*100_dp,"ETA:",eta/60.0,"min",char(27)//'[1A'//achar(13)
          flush(output_unit) !! to make sure it prints immediately          
       end if

       !!--------------------------------------------------------------------------------
       eludzh(a) = sum( u1(iv(a),iv) * u2(jv(a),jv) * u3(kv(a),kv) * diaginvlambda * zh )
       !! element of the posterior mean
       postm(a) = mprior(a) + eludzh(a)
       !!--------------------------------------------------------------------------------
       
    end do
    !$OMP END DO
    endt = omp_get_wtime()
    !$OMP END PARALLEL
    write(output_unit,*)
    print*,"posteriormean():",endt-firststartt," OMP wall clock time"

  end subroutine posteriormean
  
  !!!==================================================
  !-------------------------------------------------
  !
  !
  ! REVISION HISTORY:
  ! 
  !-------------------------------------------------  
  subroutine symgeneigvv(A,uplo,Bpd,lambda,U)
    real(dp),intent(in) :: A(:,:),Bpd(:,:)
    character,intent(in) :: uplo
    real(dp),intent(out) :: lambda(:),U(:,:)
    integer :: n,lwork,itype,info
    character :: jobz
    real(dp),allocatable :: work(:),tmpB(:,:)
    ! DSYGV computes all the eigenvalues, and optionally, the eigenvectors
    ! of a real generalized symmetric-definite eigenproblem, of the form
    ! A*x=(lambda)*B*x,  A*Bx=(lambda)*x,  or B*A*x=(lambda)*x.
    ! Here A and B are assumed to be symmetric and B is also
    ! positive definite.
    itype = 3
    jobz = 'V'
    n=size(a,1)
    if ( (size(a,2)/=n) .or. (size(bpd,1)/=n) .or. (size(bpd,2)/=n) ) then
       write(*,*) "symgeneigvv(): Matrices have wrong sizes"
       stop 
    end if
    allocate(tmpb(n,n))
    u = a
    tmpb = bpd
    lwork=3*n-1
    ! lwork=-1 
    ! allocate(work(1))
    ! call DSYGV(itype,jobz,uplo,nm1,U1,nm1,choCm1,nm1,lambda1,work,lwork,info)
    ! deallocate(work)
    allocate(work(lwork))
    call dsygv(itype,jobz,uplo,n,u,n,tmpb,n,lambda,work,lwork,info)
    if (info /= 0) stop 'symgeneigvv(): Matrix eigendecomposition failed!'
  end subroutine symgeneigvv

  !!!==================================================
  !-------------------------------------------------
  !
  !
  !
  ! REVISION HISTORY:
  ! 
  !-------------------------------------------------
  subroutine solvels( A, B, sol )
    real(dp),intent(in) :: A(:,:),B(:,:)
    real(dp),intent(out) :: sol(:,:)
    integer,allocatable :: ipiv(:)
    integer :: n,nrhs,lda,ldb
    integer :: info
    real(dp),allocatable :: A2(:,:)
    ! External procedures defined in LAPACK
    external dgesv
    n = size(a,1)
    nrhs = size(b,2)
    allocate(a2(n,n))
    a2 = a
    lda = size(a,1)
    allocate(ipiv(n))
    ldb = size(b,1)
    sol = b
    call dgesv(n,nrhs,a2,lda,ipiv,sol,ldb,info)
    if (info/=0) then
       write(*,*) "linear system solver failed...'"
       print*,'info: ',info
       stop
    end if
  end subroutine solvels

  !!!============================================================
  !-------------------------------------------------
  !
  !
  ! REVISION HISTORY:
  ! 
  !-------------------------------------------------
  subroutine symsolvels( A, B, sol )
    real(dp),intent(in) :: A(:,:),B(:,:)
    real(dp),intent(out) :: sol(:,:)
    integer,parameter :: wmax=3000 ! max size allowed for work
    integer,allocatable :: ipiv(:)
    integer :: n,nrhs,lda,ldb,lwork
    character :: uplo
    real(dp),allocatable :: work(:),A2(:,:)
    integer :: info
    ! External procedures defined in LAPACK
    external dsysv
    uplo = 'L'
    n = size(a,1)
    nrhs = size(b,2)
    allocate(a2(n,n))
    a2 = a
    lda = size(a,1)
    allocate(ipiv(n))
    ldb = size(b,1)
    ! copy B to sol to avoid it to be overwritten
    sol = b
    ! apparently work must be allocated for the query to dsysv to work...
    allocate(work(1)) 
    lwork = -1
    call dsysv(uplo,n,nrhs,a2,lda,ipiv,sol,ldb,work,lwork,info)
    lwork = min(wmax,int(work(1)))
    !! now lwork is assigned, so actually calculate eigvec/val
    deallocate(work)
    allocate(work(lwork))
    call dsysv(uplo,n,nrhs,a2,lda,ipiv,sol,ldb,work,lwork,info)
    if (info/=0) then
       write(*,*) "linear system solver failed...'"
       print*,'info: ',info
       stop
    end if
  end subroutine symsolvels

  !!====================================================
 
end module kronlininv

!!====================================================
!!====================================================