Skip to content

Commit

Permalink
Merge pull request #788 from brian-eaton/lapack_invert
Browse files Browse the repository at this point in the history 
cam6_3_107: Reimplement zonal_mean_mod::Invert_Matrix using LAPACK DGESV
  • Loading branch information
@brian-eaton
brian-eaton authored Apr 18, 2023
2 parents 0470324 + d0f1303 commit fc3cc80
Show file tree
Hide file tree
Showing 2 changed files with 98 additions and 118 deletions.
79 changes: 79 additions & 0 deletions doc/ChangeLog
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,84 @@
===============================================================

Tag name: cam6_3_107
Originator(s): eaton
Date: Tue Apr 18 10:27:45 AM EDT 2023
One-line Summary: Reimplement zonal_mean_mod::Invert_Matrix using LAPACK DGESV
Github PR URL: https://github.com/ESCOMP/CAM/pull/788

Purpose of changes (include the issue number and title text for each relevant GitHub issue):

The Invert_Matrix subroutine in module zonal_mean_mod has been
reimplemented using the LAPACK subroutine DGESV.

Resolves:
. Replace "Invert_Matrix" subroutine in "zonal_mean_mod.F90" with LAPACK version #736
(https://github.com/ESCOMP/CAM/issues/736)
. Bug in zonal mean "Invert_Matrix" subroutine #745
(https://github.com/ESCOMP/CAM/issues/745)

Describe any changes made to build system: none

Describe any changes made to the namelist: none

List any changes to the defaults for the boundary datasets: none

Describe any substantial timing or memory changes: not tested

Code reviewed by: fvitt, cacraigucar, peverwhee, patcal

List all files eliminated: none

List all files added and what they do: none

List all existing files that have been modified, and describe the changes:

src/utils/zonal_mean_mod.F90
. modify subroutine Invert_Matrix to use LAPACK DGESV routine.

If there were any failures reported from running test_driver.sh on any test
platform, and checkin with these failures has been OK'd by the gatekeeper,
then copy the lines from the td.*.status files for the failed tests to the
appropriate machine below. All failed tests must be justified.

cheyenne/intel/aux_cam:
ERP_Ln9_Vnuopc.ne30pg3_ne30pg3_mg17.FW2000climo.cheyenne_intel.cam-outfrq9s_wcm_ne30 (Overall: FAIL) details:
- pre-existing failure

ERS_Ld3_Vnuopc.f10_f10_mg37.F1850.cheyenne_intel.cam-outfrq1d_14dec_ghg_cam_dev (Overall: DIFF) details:
- the zonal mean (*zm) fields in the h1 file have 2-3 significant digits of agreement

izumi/nag/aux_cam:
DAE_Vnuopc.f45_f45_mg37.FHS94.izumi_nag.cam-dae (Overall: FAIL) details:
- pre-existing failure

SMS_D_Ln6_Vnuopc.ne5_ne5_mg37.QPWmaC4.izumi_nag.cam-outfrq3s_physgrid_tem (Overall: DIFF) details:
- the zonal mean (*zm) fields in the h2 file have 2-3 significant digits of agreement

izumi/gnu/aux_cam:
ERC_D_Ln9_Vnuopc.ne5_ne5_mg37.QPC4.izumi_gnu.cam-outfrq3s_nudging_ne5_L26 (Overall: DIFF) details:
- diffs in all fields in h0 and h1 files. This is expected because the
the nudging uses the zonal_mean_mod code if Nudge_ZonalFilter is true,
which it is in this test.

CAM tag used for the baseline comparison tests if different than previous
tag:

Summarize any changes to answers:

. The only changes to answers are for runs using nudging with
Nudge_ZonalFilter set to true.

. The zonal mean diagnostic output only agrees with previous output to a
couple of significant figures. That's because previous output was
produced with a bug in the Invert_Matrix subroutine (issue #745).
Independent testing showed that fixing that bug and comparing with the
new version of Invert_Matrix yields roundoff level differences (15-16
significant digits of agreement).

===============================================================
===============================================================

Tag name: cam6_3_106
Originator(s): cacraig, fvitt, eaton
Date: Thu Apr 6 07:10:24 PM EDT 2023
Expand Down
137 changes: 19 additions & 118 deletions src/utils/zonal_mean_mod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1318,147 +1318,48 @@ subroutine Invert_Matrix(I_Mat,Nbas,O_InvMat)
!
! Invert_Matrix: Given the NbasxNbas matrix, calculate and return
! the inverse of the matrix.
!
! Implemented with the LAPACK DGESV routine.
!
!====================================================================
real(r8),parameter:: TINY = 1.e-20_r8
!
! Passed Variables
!------------------
real(r8),intent(in ):: I_Mat (:,:)
integer ,intent(in ):: Nbas
real(r8),intent(out):: O_InvMat(:,:)
real(r8), intent(inout) :: I_Mat(:,:) ! input matrix contains P*L*U
! decomposition on output
integer, intent(in) :: Nbas
real(r8), intent(out) :: O_InvMat(:,:)
!
! Local Values
!-------------
real(r8),allocatable:: Mwrk(:,:),Rscl(:)
integer ,allocatable:: Indx(:)
real(r8):: Psgn,Mmax,Mval,Sval
integer :: ii,jj,kk,ndx,i2,ii_max, astat
integer, allocatable :: Indx(:) ! pivot indices
integer :: astat, ii
character(len=*), parameter :: subname = 'Invert_Matrix'
character(len=80) :: msg

external DGESV

! Allocate work space
!---------------------
allocate(Mwrk(Nbas,Nbas), stat=astat)
call handle_allocate_error(astat, subname, 'Mwrk')
allocate(Rscl(Nbas), stat=astat)
call handle_allocate_error(astat, subname, 'Rscl')
allocate(Indx(Nbas), stat=astat)
call handle_allocate_error(astat, subname, 'Indx')

! Copy the Input matrix so it can be decomposed
!-------------------------------------------------
Mwrk(1:Nbas,1:Nbas) = I_Mat(1:Nbas,1:Nbas)

! Initailize Row scales
!----------------------
Psgn = 1._r8
do ii=1,Nbas
Mmax = 0._r8
do jj=1,Nbas
if(abs(Mwrk(ii,jj))>Mmax) Mmax = abs(Mwrk(ii,jj))
end do
if(Mmax==0._r8) then
call endrun('Invert_Matrix: Singular Matrix')
endif
Rscl(ii) = 1._r8/Mmax
end do

! Decompose the matrix
!-----------------------
do jj=1,Nbas

if(jj>1) then
do ii=1,(jj-1)
Sval = Mwrk(ii,jj)
if(ii>1) then
do kk=1,(ii-1)
Sval = Sval - Mwrk(ii,kk)*Mwrk(kk,jj)
end do
Mwrk(ii,jj) = Sval
endif
end do
endif

Mmax = 0._r8
do ii=jj,Nbas
Sval = Mwrk(ii,jj)
if(jj>1) then
do kk=1,(jj-1)
Sval = Sval - Mwrk(ii,kk)*Mwrk(kk,jj)
end do
Mwrk(ii,jj) = Sval
endif
Mval = Rscl(ii)*abs(Sval)
if(Mval>=Mmax) then
ii_max = ii
Mmax = Mval
endif
end do

if(jj/=ii_max) then
do kk=1,Nbas
Mval = Mwrk(ii_max,kk)
Mwrk(ii_max,kk) = Mwrk(jj,kk)
Mwrk(jj,kk) = Mval
end do
Psgn = -Psgn
Rscl(ii_max) = Rscl(jj)
endif

Indx(jj) = ii_max
if(jj/=Nbas) then
if(Mwrk(jj,jj)==0._r8) Mwrk(jj,jj) = TINY
Mval = 1._r8/Mwrk(jj,jj)
do ii=(jj+1),Nbas
Mwrk(ii,jj) = Mwrk(ii,jj)*Mval
end do
endif

end do ! jj=1,Nbas

if(Mwrk(Nbas,Nbas)==0._r8) Mwrk(Nbas,Nbas) = TINY

! Initialize the inverse array with the identity matrix
!-------------------------------------------------------
O_InvMat(:,:) = 0._r8
do ii=1,Nbas
O_InvMat(ii,ii) = 1._r8
end do

! Back substitution to construct the inverse
!---------------------------------------------
do kk=1,Nbas

i2 = 0
do ii=11,Nbas
ndx = Indx(ii)
Sval = O_InvMat(ndx,kk)
O_InvMat(ndx,kk) = O_InvMat(ii,kk)
if(i2/=0) then
do jj=i2,(ii-1)
Sval = Sval - Mwrk(ii,jj)*O_InvMat(jj,kk)
end do
elseif(Sval/=0._r8) then
i2 = ii
endif
O_InvMat(ii,kk) = Sval
end do
call DGESV(Nbas, Nbas, I_Mat, Nbas, Indx, O_InvMat, Nbas, astat)

do ii=Nbas,1,-1
Sval = O_InvMat(ii,kk)
if(ii<Nbas) then
do jj=(ii+1),Nbas
Sval = Sval - Mwrk(ii,jj)*O_InvMat(jj,kk)
end do
endif
O_InvMat(ii,kk) = Sval/Mwrk(ii,ii)
end do

end do ! kk=1,Nbas
if (astat < 0) then
write(msg, '(a, i1, a)') 'argument # ', abs(astat), ' has an illegal value'
call endrun(subname//': DGESV error return: '//msg)
else if (astat > 0) then
call endrun(subname//': DGESV error return: matrix is singular')
end if

! Clean up this mess
!---------------------
deallocate(Mwrk)
deallocate(Rscl)
deallocate(Indx)

end subroutine Invert_Matrix
Expand Down

0 comments on commit fc3cc80

Please sign in to comment.