This routine is a modification of the LINPACK routine DGEFA and performs an LU decomposition of an upper Hessenberg matrix A. There are two options available:
(1) performing a fresh factorization
(2) updating the LU factors by adding a row and a
column to the matrix A.
DHEFA factors an upper Hessenberg matrix by elimination.
the matrix to be factored.
the leading dimension of the array A .
the order of the matrix A .
JOB
JOB = 1 means that a fresh factorization of the
matrix A is desired.
JOB .ge. 2 means that the current factorization of A
will be updated by the addition of a row
and a column.
an upper triangular matrix and the multipliers which were used to obtain it.
The factorization can be written A = L*U where L is a product of permutation and unit lower triangular matrices and U is upper triangular.
an integer vector of pivot indices.
INFO
= 0 normal value.
= k if U(k,k) .eq. 0.0 . This is not an error
condition for this subroutine, but it does
indicate that DHESL will divide by zero if called.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| real(kind=dp), | intent(inout) | :: | A(Lda,*) | |||
| integer, | intent(in) | :: | Lda | |||
| integer, | intent(in) | :: | N | |||
| integer, | intent(inout) | :: | Ipvt(*) | |||
| integer, | intent(out) | :: | Info | |||
| integer, | intent(in) | :: | Job |
| Type | Visibility | Attributes | Name | Initial | |||
|---|---|---|---|---|---|---|---|
| integer, | public | :: | j | ||||
| integer, | public | :: | k | ||||
| integer, | public | :: | km1 | ||||
| integer, | public | :: | kp1 | ||||
| integer, | public | :: | l | ||||
| integer, | public | :: | nm1 | ||||
| real(kind=dp), | public | :: | t |
subroutine dhefa(A,Lda,N,Ipvt,Info,Job) ! integer , intent(in) :: Lda real(kind=dp) , intent(inout) :: A(Lda,*) integer , intent(in) :: N integer , intent(inout) :: Ipvt(*) integer , intent(out) :: Info integer , intent(in) :: Job ! integer :: j , k , km1 , kp1 , l , nm1 real(kind=dp) :: t ! if ( Job>1 ) then ! ! The old factorization of A will be updated. A row and a column ! has been added to the matrix A. ! N-1 is now the old order of the matrix. ! nm1 = N - 1 ! ! Perform row interchanges on the elements of the new column, and ! perform elimination operations on the elements using the multipliers. ! if ( nm1>1 ) then do k = 2 , nm1 km1 = k - 1 l = Ipvt(km1) t = A(l,N) if ( l/=km1 ) then A(l,N) = A(km1,N) A(km1,N) = t endif A(k,N) = A(k,N) + A(k,km1)*t enddo endif ! ! Complete update of factorization by decomposing last 2x2 block. ! Info = 0 ! ! Find L = pivot index ! l = idamax(2,A(nm1,nm1),1) + nm1 - 1 Ipvt(nm1) = l ! ! Zero pivot implies this column already triangularized ! if ( A(l,nm1)==0.0D0 ) then Info = nm1 else ! ! Interchange if necessary ! if ( l/=nm1 ) then t = A(l,nm1) A(l,nm1) = A(nm1,nm1) A(nm1,nm1) = t endif ! ! Compute multipliers ! t = -1.0D0/A(nm1,nm1) A(N,nm1) = A(N,nm1)*t ! ! Row elimination with column indexing ! t = A(l,N) if ( l/=nm1 ) then A(l,N) = A(nm1,N) A(nm1,N) = t endif A(N,N) = A(N,N) + t*A(N,nm1) endif Ipvt(N) = N if ( A(N,N)==0.0D0 ) Info = N else ! ! A new facorization is desired. This is essentially the LINPACK ! code with the exception that we know there is only one nonzero ! element below the main diagonal. ! ! Gaussian elimination with partial pivoting ! Info = 0 nm1 = N - 1 if ( nm1>=1 ) then do k = 1 , nm1 kp1 = k + 1 ! ! Find L = pivot index ! l = idamax(2,A(k,k),1) + k - 1 Ipvt(k) = l ! ! Zero pivot implies this column already triangularized ! if ( A(l,k)==0.0D0 ) then Info = k else ! ! Interchange if necessary ! if ( l/=k ) then t = A(l,k) A(l,k) = A(k,k) A(k,k) = t endif ! ! Compute multipliers ! t = -1.0D0/A(k,k) A(k+1,k) = A(k+1,k)*t ! ! Row elimination with column indexing ! do j = kp1 , N t = A(l,j) if ( l/=k ) then A(l,j) = A(k,j) A(k,j) = t endif call daxpy(N-k,t,A(k+1,k),1,A(k+1,j),1) enddo endif enddo endif Ipvt(N) = N if ( A(N,N)==0.0D0 ) Info = N return endif end subroutine dhefa