KronLinInv  0.3
Kronecker-product-based linear inversion

◆ calcfactors()

subroutine, public 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 187 of file kronlininv.f08.

187 
188  real(dp),intent(in) :: g1(:,:),g2(:,:),g3(:,:),&
189  cm1(:,:),cm2(:,:),cm3(:,:),&
190  cd1(:,:),cd2(:,:),cd3(:,:)
191 
192  real(dp),intent(out) :: u1(:,:),u2(:,:),u3(:,:),&
193  iucm1(:,:),iucm2(:,:),iucm3(:,:),&
194  iucmgticd1(:,:),iucmgticd2(:,:),iucmgticd3(:,:),&
195  diaginvlambda(:)
196 
197  real(dp),allocatable :: lambda1(:),lambda2(:),lambda3(:)
198  real(dp),allocatable :: icdg1(:,:),icdg2(:,:),icdg3(:,:),&
199  gticdg1(:,:),gticdg2(:,:),gticdg3(:,:),&
200  ide1(:,:),ide2(:,:),ide3(:,:),&
201  icd1(:,:),icd2(:,:),icd3(:,:)
202  integer :: i,j,k,p,nm,nm1,nm2,nm3,nd1,nd2,nd3
203  character :: uplo
204 
205  !! Lapack routines
206  external dsygv
207 
208  nm1 = size(cm1,1)
209  nm2 = size(cm2,1)
210  nm3 = size(cm3,1)
211  nd1 = size(cd1,1)
212  nd2 = size(cd2,1)
213  nd3 = size(cd3,1)
214  nm = nm1*nm2*nm3
215 
216  !! compute preliminary stuff
217  print*,'calcfactors(): compute preliminary stuff'
218  allocate(icdg1(nd1,nm1),icdg2(nd2,nm2),icdg3(nd3,nm3))
219  call symsolvels( cd1, g1, icdg1 )
220  call symsolvels( cd2, g2, icdg2 )
221  call symsolvels( cd3, g3, icdg3 )
222 
223  !allocate(GtiCdG1(nm1,nm1),GtiCdG2(nm2,nm2),GtiCdG3(nm3,nm3))
224  gticdg1 = matmul( transpose(g1), icdg1 )
225  gticdg2 = matmul( transpose(g2), icdg2 )
226  gticdg3 = matmul( transpose(g3), icdg3 )
227 
228  !!--------------------------
229  !! fa
230  !!--------------------------
231  !! compute eigendecomposition
232  print*,'calcfactors(): compute fa'
233  allocate(lambda1(nm1),lambda2(nm2),lambda3(nm3))
234  uplo = 'L'
235  call symgeneigvv(gticdg1,uplo,cm1,lambda1,u1)
236  call symgeneigvv(gticdg2,uplo,cm2,lambda2,u2)
237  call symgeneigvv(gticdg3,uplo,cm3,lambda3,u3)
238 
239  !!--------------------------
240  !! fb
241  !!--------------------------
242  !! calculates diagonal central factor
243  print*,'calcfactors():compute fb'
244  p=1
245  do i=1,nm1
246  do j=1,nm2
247  do k=1,nm3
248  diaginvlambda(p) = 1.0_dp/(1.0_dp+lambda1(i)*lambda2(j)*lambda3(k))
249  p = p+1
250  end do
251  end do
252  end do
253 
254  !!--------------------------
255  !! fc
256  !!--------------------------
257  !! computes U1^-1 Cmx, U2^-1 Cmy, U3^-1 Cmz
258  print*,'calcfactors(): compute fc'
259  call solvels( u1, cm1, iucm1 )
260  call solvels( u2, cm2, iucm2 )
261  call solvels( u3, cm3, iucm3 )
262 
263  !!--------------------------
264  !! fd
265  !!--------------------------
266  print*,'calcfactors(): compute fd'
267  allocate(ide1(nd1,nd1),ide2(nd2,nd2),ide3(nd3,nd3))
268  ide1 = 0.0_dp ! Initialize the array.
269  ide2 = 0.0_dp
270  ide3 = 0.0_dp
271  forall(i = 1:nd1) ide1(i,i) = 1.0_dp ! Set the diagonal.
272  forall(i = 1:nd2) ide2(i,i) = 1.0_dp
273  forall(i = 1:nd3) ide3(i,i) = 1.0_dp
274  allocate(icd1(nd1,nd1),icd2(nd2,nd2),icd3(nd3,nd3))
275  call symsolvels( cd1, ide1, icd1 )
276  call symsolvels( cd2, ide2, icd2 )
277  call symsolvels( cd3, ide3, icd3 )
278 
279  iucmgticd1 = matmul( iucm1, matmul( transpose(g1), icd1 ))
280  iucmgticd2 = matmul( iucm2, matmul( transpose(g2), icd2 ))
281  iucmgticd3 = matmul( iucm3, matmul( transpose(g3), icd3 ))
282 
283  ! print*,
284  ! print*,shape(U1),shape(U2),shape(U3)
285  ! print*,shape(iUCm1),shape(iUCm2),shape(iUCm3)
286  ! print*,shape(iUCmGtiCd1),shape(iUCmGtiCd2),shape(iUCmGtiCd3)
287  ! print*,
288  print*,'calcfactors(): end '
Here is the call graph for this function:
Here is the caller graph for this function: