inverse_matrix.F90

Go to the documentation of this file.

!SFX_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
!SFX_LIC This is part of the SURFEX software governed by the CeCILL-C licence
!SFX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt  
!SFX_LIC for details. version 1.
SUBROUTINE inverse_matrix(N,A,P)
!--------------------------------------------------------
!
! Explicit inversed matrix after Cholesky decomposition
!
!--------------------------------------------------------
 USE yomhook   ,ONLY : lhook,   dr_hook
 USE parkind1  ,ONLY : jprb
!
 IMPLICIT NONE
 INTEGER, INTENT(IN) :: n
 REAL, DIMENSION (N,N),  INTENT(INOUT) :: a
 REAL, DIMENSION (N),  INTENT(IN)      :: p
 REAL zsum
 INTEGER :: i, j, k
 REAL(KIND=JPRB) :: zhook_handle

 IF (lhook) CALL dr_hook('INVERSE_MATRIX',0,zhook_handle)
 DO i=1,n
   a(i,i)=1./p(i)
   DO j=i+1,n
     zsum = 0.
     DO k=i,j-1
       zsum = zsum - a(j,k)*a(k,i)
     ENDDO
     a(j,i) = zsum/p(j)
   ENDDO
 ENDDO  
 DO i=1,n
   DO j=i+1,n
      a(i,j) =0.0
   ENDDO
 ENDDO
 a = matmul(transpose(a),a)

 IF (lhook) CALL dr_hook('INVERSE_MATRIX',1,zhook_handle)

END SUBROUTINE inverse_matrix