KronLinInv  0.3
Kronecker-product-based linear inversion

◆ calcfactors()

subroutine, public ompi_kronlininv::calcfactors ( real(dp), dimension(:,:), intent(in)  G1,
real(dp), dimension(:,:), intent(in)  G2,
real(dp), dimension(:,:), intent(in)  G3,
real(dp), dimension(:,:), intent(in)  Cm1,
real(dp), dimension(:,:), intent(in)  Cm2,
real(dp), dimension(:,:), intent(in)  Cm3,
real(dp), dimension(:,:), intent(in)  Cd1,
real(dp), dimension(:,:), intent(in)  Cd2,
real(dp), dimension(:,:), intent(in)  Cd3,
real(dp), dimension(:,:), intent(out)  U1,
real(dp), dimension(:,:), intent(out)  U2,
real(dp), dimension(:,:), intent(out)  U3,
real(dp), dimension(:), intent(out)  diaginvlambda,
real(dp), dimension(:,:), intent(out)  iUCm1,
real(dp), dimension(:,:), intent(out)  iUCm2,
real(dp), dimension(:,:), intent(out)  iUCm3,
real(dp), dimension(:,:), intent(out)  iUCmGtiCd1,
real(dp), dimension(:,:), intent(out)  iUCmGtiCd2,
real(dp), dimension(:,:), intent(out)  iUCmGtiCd3 
)

Computes the factors necessary to solve the inverse problem.

The factors are the ones to be stored to subsequently calculate posterior mean and covariance. First an eigen decomposition is performed, to get

\[ \mathbf{U}_1 \mathbf{\Lambda}_1 \mathbf{U}_1^{-1} = \mathbf{C}_{\rm{M}}^{\rm{x}} (\mathbf{G}^{\rm{x}})^{\sf{T}} (\mathbf{C}_{\rm{D}}^{\rm{x}})^{-1} \mathbf{G}^{\rm{x}} \]

\[ \mathbf{U}_2 \mathbf{\Lambda}_2 \mathbf{U}_2^{-1} = \mathbf{C}_{\rm{M}}^{\rm{y}} (\mathbf{G}^{\rm{y}})^{\sf{T}} (\mathbf{C}_{\rm{D}}^{\rm{y}})^{-1} \mathbf{G}^{\rm{y}} \]

\[ \mathbf{U}_3 \mathbf{\Lambda}_3 \mathbf{U}_3^{-1} = \mathbf{C}_{\rm{M}}^{\rm{z}} (\mathbf{G}^{\rm{z}})^{\sf{T}} (\mathbf{C}_{\rm{D}}^{\rm{z}})^{-1} \mathbf{G}^{\rm{z}} \]

The principal factors involved in the computation of the posterior covariance and mean are

\[ F_{\sf{A}} = \mathbf{U}_1 \otimes \mathbf{U}_2 \otimes \mathbf{U}_3 \]

\[ F_{\sf{B}} = \big( \mathbf{I} + \mathbf{\Lambda}_1 \! \otimes \! \mathbf{\Lambda}_2 \! \otimes \! \mathbf{\Lambda}_3 \big)^{-1} \]

\[ F_{\sf{C}} = \mathbf{U}_1^{-1} \mathbf{C}_{\rm{M}}^{\rm{x}} \otimes \mathbf{U}_2^{-1} \mathbf{C}_{\rm{M}}^{\rm{y}} \otimes \mathbf{U}_3^{-1} \mathbf{C}_{\rm{M}}^{\rm{z}} \]

\[ F_{\sf{D}} = \left( \mathbf{U}_1^{-1} \mathbf{C}_{\rm{M}}^{\rm{x}} (\mathbf{G}^{\rm{x}})^{\sf{T}} (\mathbf{C}_{\rm{D}}^{\rm{x}})^{-1} \right) \! \otimes \left( \mathbf{U}_2^{-1} \mathbf{C}_{\rm{M}}^{\rm{y}} (\mathbf{G}^{\rm{y}})^{\sf{T}} (\mathbf{C}_{\rm{D}}^{\rm{y}})^{-1} \right) \! \otimes \left( \mathbf{U}_3^{-1} \mathbf{C}_{\rm{M}}^{\rm{z}} (\mathbf{G}^{\rm{z}})^{\sf{T}} (\mathbf{C}_{\rm{D}}^{\rm{z}})^{-1} \right) \]

http://www.netlib.org/lapack/lug/node54.html ?sygst Reduces a real symmetric-definite generalized eigenvalue problem to the standard form.
\(B A z = \lambda z\) B = LLT C = LT A L z = L y

  • A is symmetric
  • B is symmetric, positive definite

This routines calculates

Parameters
[in]G1,G2,G3The 3 forward modeling matrices \( \mathbf{G} = \mathbf{G_1} \otimes \mathbf{G_2} \otimes \mathbf{G_3} \)
[in]Cm1,Cm2,Cm3,Cd1,Cd2,Cd3The 6 covariance matrices
[out]U1,U2,U3\( \mathbf{U}_1 \), \( \mathbf{U}_2 \), \( \mathbf{U}_3 \) of \( F_{\sf{A}} \)
[out]diaginvlambda\( F_{\sf{B}} \)
[out]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}} \)
[out]iUCmGtiCd1,iUCmGtiCd1,iUCmGtiCd1\( \mathbf{U}_1^{-1} \mathbf{C}_{\rm{M}}^{\rm{x}} (\mathbf{G}^{\rm{x}})^{\sf{T}}(\mathbf{C}_{\rm{D}}^{\rm{x}})^{-1} \), \( \mathbf{U}_2^{-1} \mathbf{C}_{\rm{M}}^{\rm{y}} (\mathbf{G}^{\rm{y}})^{\sf{T}} (\mathbf{C}_{\rm{D}}^{\rm{y}})^{-1}\), \( \mathbf{U}_3^{-1} \mathbf{C}_{\rm{M}}^{\rm{z}} (\mathbf{G}^{\rm{z}})^{\sf{T}} (\mathbf{C}_{\rm{D}}^{\rm{z}})^{-1} \) of \( F_{\sf{D}} \)
Date
10/8/2016 - Initial Version

Definition at line 461 of file ompi_kronlininv.f08.

461 
462  real(dp),intent(in) :: g1(:,:),g2(:,:),g3(:,:),&
463  cm1(:,:),cm2(:,:),cm3(:,:),&
464  cd1(:,:),cd2(:,:),cd3(:,:)
465 
466  real(dp),intent(out) :: u1(:,:),u2(:,:),u3(:,:),&
467  iucm1(:,:),iucm2(:,:),iucm3(:,:),&
468  iucmgticd1(:,:),iucmgticd2(:,:),iucmgticd3(:,:),&
469  diaginvlambda(:)
470 
471  real(dp),allocatable :: lambda1(:),lambda2(:),lambda3(:)
472  real(dp),allocatable :: icdg1(:,:),icdg2(:,:),icdg3(:,:),&
473  gticdg1(:,:),gticdg2(:,:),gticdg3(:,:),&
474  ide1(:,:),ide2(:,:),ide3(:,:),&
475  icd1(:,:),icd2(:,:),icd3(:,:)
476  integer :: i,j,k,p,nm,nm1,nm2,nm3,nd1,nd2,nd3
477  character :: uplo
478 
479  !! Lapack routines
480  external dsygv
481 
482  nm1 = size(cm1,1)
483  nm2 = size(cm2,1)
484  nm3 = size(cm3,1)
485  nd1 = size(cd1,1)
486  nd2 = size(cd2,1)
487  nd3 = size(cd3,1)
488  nm = nm1*nm2*nm3
489 
490  !! compute preliminary stuff
491  print*,'calcfactors(): compute preliminary stuff'
492  allocate(icdg1(nd1,nm1),icdg2(nd2,nm2),icdg3(nd3,nm3))
493  call symsolvels( cd1, g1, icdg1 )
494  call symsolvels( cd2, g2, icdg2 )
495  call symsolvels( cd3, g3, icdg3 )
496 
497  allocate(gticdg1(nm1,nm1),gticdg2(nm2,nm2),gticdg3(nm3,nm3))
498  gticdg1 = matmul( transpose(g1), icdg1 )
499  gticdg2 = matmul( transpose(g2), icdg2 )
500  gticdg3 = matmul( transpose(g3), icdg3 )
501 
502  !!--------------------------
503  !! fa
504  !!--------------------------
505  !! compute eigendecomposition
506  print*,'calcfactors(): compute fa'
507  allocate(lambda1(nm1),lambda2(nm2),lambda3(nm3))
508  uplo = 'L'
509  call symgeneigvv(gticdg1,uplo,cm1,lambda1,u1)
510  call symgeneigvv(gticdg2,uplo,cm2,lambda2,u2)
511  call symgeneigvv(gticdg3,uplo,cm3,lambda3,u3)
512 
513  !!--------------------------
514  !! fb
515  !!--------------------------
516  !! calculates diagonal central factor
517  print*,'calcfactors(): compute fb'
518  p=1
519  do i=1,nm1
520  do j=1,nm2
521  do k=1,nm3
522  diaginvlambda(p) = 1.0_dp/(1.0_dp+lambda1(i)*lambda2(j)*lambda3(k))
523  p = p+1
524  end do
525  end do
526  end do
527 
528  !!--------------------------
529  !! fc
530  !!--------------------------
531  !! computes U1^-1 Cmx, U2^-1 Cmy, U3^-1 Cmz
532  print*,'calcfactors(): compute fc'
533  call solvels( u1, cm1, iucm1 )
534  call solvels( u2, cm2, iucm2 )
535  call solvels( u3, cm3, iucm3 )
536 
537 
538  !!--------------------------
539  !! fd
540  !!--------------------------
541  print*,'calcfactors(): compute fd'
542  allocate(ide1(nd1,nd1),ide2(nd2,nd2),ide3(nd3,nd3))
543  ide1 = 0.0_dp ! Initialize the array.
544  ide2 = 0.0_dp
545  ide3 = 0.0_dp
546  forall(i = 1:nd1) ide1(i,i) = 1.0_dp ! Set the diagonal.
547  forall(i = 1:nd2) ide2(i,i) = 1.0_dp
548  forall(i = 1:nd3) ide3(i,i) = 1.0_dp
549  allocate(icd1(nd1,nd1),icd2(nd2,nd2),icd3(nd3,nd3))
550  call symsolvels( cd1, ide1, icd1 )
551  call symsolvels( cd2, ide2, icd2 )
552  call symsolvels( cd3, ide3, icd3 )
553 
554  iucmgticd1 = matmul( iucm1, matmul( transpose(g1), icd1 ))
555  iucmgticd2 = matmul( iucm2, matmul( transpose(g2), icd2 ))
556  iucmgticd3 = matmul( iucm3, matmul( transpose(g3), icd3 ))
557 
558  print*,"calcfactors(): end "
Here is the call graph for this function:
Here is the caller graph for this function: