Investigate matrix size question in ropp_1dvar

[Ian Culverwell]

Joe Nielsen (DMI) reports the following issue.

Ian, I do not expect you to know the details of this, but maybe you 
can give me a hint or suggestion? The problem is the use of 
State1DFM varaibles, especially the  state % cov property/variable;

In ropp_fm_types it is defined as a compressed matrix as I 
understand it:

     TYPE(matrix_pp)                 :: cov                     ! Error covariance

 In ropp_1dvar/common/ropp_1dvar_diagnostics.f90 the state variables 
is declared like this:

 TYPE(State1DFM),   INTENT(inout)      :: state

 ....

And the further down it is used like this:

! 2.2 Calculate error covariance of the solution

! ----------------------------------------------

     Bm1 = matrix_invert(state % cov)
     Rm1 = matrix_invert(obs % cov)
     Kt  = TRANSPOSE(K)
     CALL matrix_toast(Rm1, Kt, KtRm1K)

     S = Bm1 + KtRm1K
     state % cov = matrix_invert(S)

 (all the matrices are allocated as full matrices)

How can that work? matrix_invert outputs a full matrix, and 
matrix_toast understands only full matrices. Looks like an error, 
but I hope there is an explanation....

 I am suddenly getting a bit worried here because I copy/pasted this 
part into our RE1 code at some point (only for postprocessing 
though).

 thanks

 joe

[Ian Culverwell]

I replied

Anyway, it certainly does look a bit odd.  matrix_invert 
should mean matrx_invert_gen in ropp_1dvar/math/matrix/matrix_invert.f90 
if its input is a full matrix like S.  What are SHAPE(S) and SHAPE(state%cov) 
after you've set the latter to matrix_invert(S)?

[Ian Culverwell]

Joe replied

Hi Ian

It seems like ropp_1dvar_covar  returns a full covariance 
matrix - which does not make sense since the supposedly 
overlayed ropp_1dvar_covar_bg code explicitly writes a 
packed cov matrix. I have dumped B like this

   IF ( bg%state_ok ) CALL ropp_1dvar_covar (bg, config)

     Bmat=bg%cov
     write(*,*) 'shape Bmat', shape(Bmat)
     write(*,*) 'Bmat start'
     write(*,'(E15.7)') Bmat
     write(*,*) 'Bmat end'
     write(*,*) 'shape Rm', shape(Rm)

- and plottet Bmat;

Bmat is full (and correct) no doubt.

Look, I really need to move on with the GPAC 2.4 which is 
due Friday (last week).  So what I do for now is that I 
code under the assumption that bg%cov is full, which 
seems to be style other places in  ROPP. I will let you 
decide whether this needs a ROPP-ticket.

I noted this: In my first attempt  I failed to allocate 
Bmat which resulted in an runtime error saying that the 
matrix needs to be at least 275x275 (ie. full) hinting 
that the compiler thinks that the output of 
ropp_1dvar_covar contains a full cov matrix.

Obvious next move would be to check whether 
ropp_1dvar_covar_bg is ever reached. It looks to me as if 
it is somehow bypassed...

joe

[Ian Culverwell]

While we're here, ropp_1dvar/math/matrix/matrix_invert.f90 has lots of variable declarations like

  real(wp), dimension((int(sqrt(8.*size(A%d))+1)-1)/2, &
                         (int(sqrt(8.*size(A%d))+1)-1)/2) :: x

where size(A%d) = the number of elements in the packed arrays = n(n+1)/2 for an n x n unpacked array. The bit of algebra (int(sqrt(8.*size(A%d))+1)-1)/2 is presumably designed to return n. And so it does - but rather indirectly:

(int(sqrt(8.*size(A%d))+1)-1)/2 = (int(sqrt(4n(n+1))+1)-1)/2

But

4n² < 4n(n+1) < 4(n+1)²

=> 2n < sqrt(4n(n+1)) < 2(n+1)

=> 2n + 1 < sqrt(4n(n+1)) + 1 < 2(n+1) + 1

=> int(sqrt(4n(n+1))+1) = 2n + 1

=> (int(sqrt(4n(n+1))+1)-1)/2 = n

OK. But the alternative expression (nint(sqrt(8.*size(A%d)+1))-1)/2 does the job much more directly:

(nint(sqrt(8.*size(A%d)+1))-1)/2

(nint(sqrt(4n(n+1)+1))-1)/2

(nint(sqrt((2n+1)²)-1)/2

(nint(2n+1)-1)/2

n

Perhaps, since we still need a nint to handle rounding errors, there's not much in it. But at least this note might help future developers who worry that the ROPP expression is incorrect.

[Ian Culverwell]

(Better still: 4n² < 4n(n+1) < 4(n+1/2)²

=> 2n + 1 < sqrt(4n(n+1)) + 1 < 2(n+1/2) + 1 = 2n + 2

=> int(sqrt(4n(n+1))+1) = 2n + 1 as before.)

[Ian Culverwell]

It's all a result of operator overloading.

When the code says

diag % bg_refrac % cov (of type matrix_pp)  =  R (of type real(:, :))

it is ‘overloading’ the operator ‘=’.

The entry

interface assignment(=)

in ropp_1dvar/math/matrix/matrix.f90 shows that in this case we are calling

subroutine matrix_full2ppt(pp, data)

in ropp_1dvar/math/matrix/matrix_assign.f90, which in turn calls

subroutine matrix_full2pp_alloc_double(array, packed, uplo)

in ropp_1dvar/math/matrix/matrix_full2pp_alloc.f90, which finally does the packing of the elements of R into the %d component of diag % bg_refrac % cov.

(matrix_full2ppt actually calls matrix_full2pp_alloc, but this calls matrix_full2pp_alloc_double for these variables, as a result of the interface block in matrix.f90.)

It’s involved, but logical and very powerful, because another chain of calls would allow us to reverse the assignment and say

R (of type real(:, :)) = diag % bg_refrac % cov (of type matrix_pp)

Basically, the meaning of '=' in 'LHS = RHS' depends on the Fortran types of LHS and RHS.

[Ian Culverwell]

Overloading operators and generic interfaces are powerful features of fortran90, and actually we use them all over the place in ROPP, e.g. you can call ropp_io_free on an ROprof structure, a covariance matrix, a vlist of ‘extra data’ etc, all as a result of the ropp_io_free interface block in ropp_io/ropp/rop_io.f90. And we often write roprof1 = roprof2, which only makes sense because ropp_io_assign_rotype copies each element of the roprof2 structure to the corresponding element of roprof1. You can write anything1 = anything2 as long as ‘=’ has been overwritten for the fortran types of anything1 and anything2. So you can probably do it for any two ROPP matrices in ROPP – I doubt that Chris would have left holes. The compiler will let you know pretty quickly if the underlying infrastructure isn’t in place.

[Ian Culverwell]

I think this answers the question. Closing the ticket.