cam6_3_107: Reimplement zonal_mean_mod::Invert_Matrix using LAPACK DGESV

[brian-eaton]

The Invert_Matrix subroutine in module zonal_mean_mod has been reimplemented using the LAPACK subroutine DGESV. This change causes roundoff level differences in the diagnostic zonal mean fields compared to the original version of Invert_Matrix with the fix from issue #745 applied.

closes #736
fixes #745

[peverwhee]

Looks good to me!

[fvitt]

Pleasantly surprised that the differences in the zonal means are only roundoff level differences

[patcal]

Adding the dimensions to the passed arrays provides clarity and is a good safeguard against a compiler that may not like the (:,:) usage.

[brian-eaton]

@patcal, I've been thinking about your comment on using explicit shape dummy args. After a bit of searching on pros/cons of the assumed shape vs explicit shape dummy args I think I'm going to stick with assumed shape. It seems to be a more flexible option, and we could potentially move the Invert_Matrix routine into another module to be used elsewhere in the code.

[gold2718]

Adding the dimensions to the passed arrays provides clarity and is a good safeguard against a compiler that may not like the (:,:) usage.

For the record, adding the dimensions in the declaration statement disables shape checks so if the array passed into that routine does not match the declared size*, the array references will be different than in the calling routine which usually leads to unexpected results.
In my view, this makes Brian's approach safer.

*in at least one dimension that is not the last one.

[patcal]

Great, I am just happy that everyone decided that I had it coded correctly in the first place.

…

On Tue, Apr 18, 2023 at 1:53 PM goldy ***@***.***> wrote: Adding the dimensions to the passed arrays provides clarity and is a good safeguard against a compiler that may not like the (:,:) usage. For the record, adding the dimensions in the declaration statement disables shape checks so if the array passed into that routine does not match the declared size*, the array references will be different than in the calling routine which usually leads to unexpected results. In my view, this makes Brian's approach safer. *in at least one dimension that is not the last one. — Reply to this email directly, view it on GitHub <#788 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AJMQHAONQSXXAJ3OZZHESB3XB3WJ3ANCNFSM6AAAAAAW5N5D4Q> . You are receiving this because you were mentioned.Message ID: ***@***.***>

[sjsprecious]

Adding the dimensions to the passed arrays provides clarity and is a good safeguard against a compiler that may not like the (:,:) usage.

For the record, adding the dimensions in the declaration statement disables shape checks so if the array passed into that routine does not match the declared size*, the array references will be different than in the calling routine which usually leads to unexpected results. In my view, this makes Brian's approach safer.

*in at least one dimension that is not the last one.

Thanks @gold2718 for this additional information, which I did not know before. I just wanted to add another point here: based on David Appelhans's talk at GTC23, using assumed-shape variables could easily hurt the GPU performance, compared to using the variables with explicit size.

[gold2718]

I just wanted to add another point here: based on David Appelhans's talk at GTC23, using assumed-shape variables could easily hurt the GPU performance, compared to using the variables with explicit size.

This seems like a classic tension (engineering for flexibility and testability vs. engineering for performance).
@sjsprecious, do you have any data on this? It would be great to know what the penalty is on real current GPU architectures. How big is the penalty?

[sjsprecious]

According to the examples in David's talk, the performance penalty could be between 2x and 40x on NVIDIA's A100 GPU when using the assumed shape variables. You could find more details by searching his GTC23 talk titled "Best Practices for Programming GPUs using Fortran, OpenACC, and CUDA". I could send you a copy of the slides too in case it is no longer accessible for the public.

[brian-eaton]

I think in the current context using the assumed shape arrays is fine because Invert_Matrix is just a wrapper for calling DGESV. The dummy arrays in the DGESV interface are assumed size arrays which means the compiler is going to create a copy to provide contiguous storage for the array references that are passed down.

[sjsprecious]

Thanks @brian-eaton. I agreed with you that the assumed shape arrays might not be an issue here since you were calling an LAPACK interface and not using any GPU. Just want to mention this potential performance issue somewhere for record but probably I leave it at the wrong place (apologize if that is the case).

In addition, you raised another good point and I want to make it clear: according to David's talk, using assumed size arrays won't hurt the GPU performance but using assumed shape arrays will.