From b047e96899773e807f290667a9aea387d8a2f842 Mon Sep 17 00:00:00 2001 From: Brian Eaton Date: Thu, 13 Apr 2023 13:50:28 -0400 Subject: [PATCH 1/2] reimplement Invert_Matrix using LAPACK DGESV --- src/utils/zonal_mean_mod.F90 | 137 +++++------------------------------ 1 file changed, 19 insertions(+), 118 deletions(-) diff --git a/src/utils/zonal_mean_mod.F90 b/src/utils/zonal_mean_mod.F90 index 3c41a4b4ec..25e3f8564a 100644 --- a/src/utils/zonal_mean_mod.F90 +++ b/src/utils/zonal_mean_mod.F90 @@ -1318,105 +1318,32 @@ 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 @@ -1424,41 +1351,15 @@ subroutine Invert_Matrix(I_Mat,Nbas,O_InvMat) 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 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 From d0f1303ff0c3cbbcda36870870eb486d98aa4da6 Mon Sep 17 00:00:00 2001 From: Brian Eaton Date: Tue, 18 Apr 2023 10:55:54 -0400 Subject: [PATCH 2/2] update ChangeLog --- doc/ChangeLog | 79 +++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 79 insertions(+) diff --git a/doc/ChangeLog b/doc/ChangeLog index 025b32dff6..36ce62b27e 100644 --- a/doc/ChangeLog +++ b/doc/ChangeLog @@ -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