getLU Subroutine

public subroutine getLU(dimA, A, L, U, P)

Return P L U matrices of the LU decomposition, cmoputed from LAPACK Routine (if P is not passed, L=P*L)

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: dimA
real(kind=wp), intent(in) :: A(1:dimA,1:dimA)

matrix
real(kind=wp), intent(out) :: L(1:dimA,1:dimA)

L or P*L (if P omitted)
real(kind=wp), intent(out) :: U(1:dimA,1:dimA)
real(kind=wp), intent(out), optional :: P(1:dimA,1:dimA)

Calls

proc~~getlu~~CallsGraph proc~getlu getLU dgetrf dgetrf proc~getlu->dgetrf

Called by

proc~~getlu~~CalledByGraph proc~getlu getLU proc~sbase_init t_sBase%sBase_init proc~sbase_init->proc~getlu proc~sbase_test sBase_test proc~sbase_init->proc~sbase_test proc~sbase_copy t_sBase%sBase_copy proc~sbase_copy->proc~sbase_init proc~sbase_new sBase_new proc~sbase_new->proc~sbase_init proc~base_copy t_base%base_copy proc~base_copy->proc~sbase_copy proc~base_new Base_new proc~base_new->proc~sbase_new proc~readstatefilefromascii ReadStateFileFromASCII proc~readstatefilefromascii->proc~sbase_new proc~readstatefilefromascii->proc~base_new proc~sbase_test->proc~sbase_new interface~readstate ReadState interface~readstate->proc~readstatefilefromascii proc~init_base Init_Base proc~init_base->proc~base_new proc~initmhd3d t_functional_mhd3d%InitMHD3D proc~initmhd3d->proc~base_new proc~transform_sfl_init t_transform_sfl%transform_SFL_init proc~transform_sfl_init->proc~base_new proc~init Init proc~init->proc~initmhd3d proc~init_gvec_to_jorek init_gvec_to_jorek proc~init_gvec_to_jorek->interface~readstate proc~init_gvec_to_jorek->proc~init_base proc~restartfromstate RestartFromState proc~restartfromstate->interface~readstate proc~rungvec rungvec proc~rungvec->proc~initmhd3d proc~transform_sfl_new transform_sfl_new proc~transform_sfl_new->proc~transform_sfl_init program~gvec_post GVEC_POST program~gvec_post->proc~initmhd3d proc~start_rungvec start_rungvec proc~start_rungvec->proc~rungvec program~gvec GVEC program~gvec->proc~rungvec

Source Code

SUBROUTINE getLU(dimA,A,L,U,P)
! MODULES
IMPLICIT NONE
!-----------------------------------------------------------------------------------------------------------------------------------
! INPUT VARIABLES
INTEGER, INTENT(IN) :: dimA
REAL(wp),INTENT(IN) :: A(1:dimA,1:dimA) !! matrix
!-----------------------------------------------------------------------------------------------------------------------------------
! OUTPUT VARIABLES
REAL(wp),INTENT(OUT) :: L(1:dimA,1:dimA)   !! L or P*L (if P omitted)
REAL(wp),INTENT(OUT) :: U(1:dimA,1:dimA)
REAL(wp),INTENT(OUT),OPTIONAL :: P(1:dimA,1:dimA)
!-----------------------------------------------------------------------------------------------------------------------------------
! LOCAL VARIABLES
! External procedures defined in LAPACK
EXTERNAL DGETRF
! LOCAL VARIABLES
REAL(wp)    :: tmp,Atmp(dimA,dimA)
INTEGER     :: ipiv(dimA)  ! pivot indices
INTEGER     :: i,j,info
!===================================================================================================================================
Atmp=A

CALL DGETRF(dimA, dimA, Atmp, dimA, ipiv, info)

IF(info.NE.0)THEN
   CALL abort(__STAMP__,&
              'Matrix is numerically singular!')
END IF
!upper diagonal part
U=0.
DO i=1,dimA
  DO j=i,dimA
    U(i,j)=Atmp(i,j)
  END DO
END DO
!lower diagonal part, known: 1 on the diagonal
L=0.
DO i=1,dimA
  L(i,i)=1.
  DO j=1,i-1
    L(i,j)=Atmp(i,j)
  END DO
END DO

IF(PRESENT(P))THEN
  P=0.
  DO i=1,dimA
    P(i,i)=1.
  END DO

  !pivoting of rows in L is pivoting of columns in P
  DO i=1,dimA
    DO j=1,dimA
      tmp=P(j,i)
      P(j,i)=P(j,ipiv(i))
      P(j,ipiv(i))=tmp
    END DO
  END DO

  !CHECK
  IF(SUM(ABS(MATMUL(P,MATMUL(L,U))-A)).GT.1e-12*dimA*dimA) THEN
    CALL abort(__STAMP__,&
               'A=P*L*U decomposition not correct')
  END IF
ELSE
  !pivoting of rows in L, backwards
  DO i=dimA,1,-1
    DO j=1,dimA
      tmp=L(i,j)
      L(i,j)=L(ipiv(i),j)
      L(ipiv(i),j)=tmp
    END DO
  END DO

  !CHECK
  IF(SUM(ABS(MATMUL(L,U)-A)).GT.1e-12*dimA*dimA) THEN
    CALL abort(__STAMP__,&
               'A=L*U decomposition not correct')
  END IF
END IF

END SUBROUTINE getLU