Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

cam6_3_107: Reimplement zonal_mean_mod::Invert_Matrix using LAPACK DGESV #788

Merged
merged 2 commits into from
Apr 18, 2023
Merged

cam6_3_107: Reimplement zonal_mean_mod::Invert_Matrix using LAPACK DGESV #788

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
b047e96
reimplement Invert_Matrix using LAPACK DGESV
brian-eaton Apr 13, 2023
d0f1303
update ChangeLog
brian-eaton Apr 18, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
79 changes: 79 additions & 0 deletions doc/ChangeLog
View file
Open in desktop
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
Open in desktop
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(:,:)
cacraigucar marked this conversation as resolved.
Show resolved Hide resolved
!
! 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