In this program I assume that n3 can be 1 or equal to n2 so, everything should work. This subroutine works fine if n3>1, but if n3=1 I have problems and I cannot figure out why. In particular, the values of x that goes from (ista+1:iend) were set to zero and I don’t know why.
I expect that if n3=1, then tmp3(1:n3) is converted to a scalar and ALL the values of x(ista:iend) will be set. Instead, it seems that the compiler performs a strange conversion between ista:iend and 1:n3.
This will never happen. If n3=1, then tmp3(1:n3) will be an array of length 1. The lhs of the statement x(ista:iend) = c2*tmp3(1:n3)**2 must be conformant (the same size) as the rhs (length 1 in this case). If not, it is a programmer error. Your compiler probably has options to print a run time error message if that is not satisfied.
As I think you know, a scalar expression would also be conformant, and if the rhs of the assignment were a scalar then it would be broadcast to all the array elements of the lhs. But an array of length 1 is still an array, not a scalar, so that does not occur. The language does not define what occurs if ista:iend on the lhs is not conformant with 1:n3 on the rhs, it is just a programmer error.
Related to this, do you know if there is a way to automatize the assignment? I don’t want each time to check if n3=1 or n3>1 and then create 2 different assignments. I was wondering if there was a fast way.
Another option might be to use allocation/reallocation on assignment for x BUT be aware this has the potential for triggering a performance degradation if the final size of x is large (as in the thousands or millions of elements). I’ve found that this techniques works well for arrays sizes of around 100 or 200 but would look at other options if x is larger. Basically you do the following
Make x allocatable
Use allocation on assignment to extend x by whatever size n3 is
The code fragment for your sub1 subroutine would look something like
subroutine sub1(n1,n2,x,y)
real(kind=8), allocatable, intent(inout) :: x(:)
real(kind=8) :: y(n1*n2)
real(kind=8) :: tmp1(n2),tmp2(n2)
real(kind=8) :: tmp3(n3)
integer :: i,ista,iend
allocate(x(0)) ! not really necessary but I had to do this to avoid segfaults once
do i=1,n1
ista = 1+(i-1)*n2
iend = i*n2
CALL CALC(tmp3(1:n3),y(ista:iend)) !calculate tmp3
x = [x, tmp3(1:n3)]
enddo
end subroutine
This technique is basically a poor mans C++ vector type.
Just make sure x is also allocatable in the routine that calls sub1. Try this but again with the caveat that you will probably take a performance hit if the size of x becomes large. If you see a large slowdown try profiling sub1 to see if its really the allocation on assignment thats causing the problem.
I’m not sure what is the question. Why would you need different assignments?
You know upon entry to sub1() that n1*n3 elements of x(:) will be set. In your original code there are instead n1*n2 elements that are referenced. Maybe the problem is that there is an inconsistency between n2 and n3? Should the output x(:) be referenced contiguously with n1*n3 elements, or should it be referenced sparsely with n1*n2 elements? In the sparse case, do the unreferenced elements need to be initialized or defined?
Here are some suggestions to consider. Don’t hardwire the real kind values in the declarations. It just makes the code difficult to understand and, if necessary, to modify. Instead, define a parameter constant. It looks like the arrays x and y are really being referenced as 2D arrays. If so, then why not just declare them that way? You will end up with something like this:
I’ve assumed a few things here, like the argument intents. Also, I assumed that your undeclared c2 was really supposed to be the declared but unreferenced c1 (the implicit none statement would have warned you about that). Your tmp1(:) and tmp2(:) were unreferenced, so I just removed them.
Hopefully your actual problem in your actual code can be determined from these comments and modifications.
Hi, yes, x,y are two 2D arrays of that shape. Just a point: if in input of the subroutine I give a flattened array, then, can I reshape it as you have done?
Yes, with explicit shape dummy argument arrays you can go back and forth. This uses storage sequence association rules, so the behavior is well defined within the language. One warning about this is that if the actual argument array is not contiguous, then the compiler must use copy-in/copy-out association for the dummy array argument, so this can result in extra overhead effort during the call.
I wanted to add one other comment. Notice in the statement x(:,i) = c1 * tmp3**2 I removed the array stride. It was originally written as tmp3(1:n3). One reason for doing that is just in case the variable n3 is changed within the CALC() call. If that were to happen, then the lhs and rhs array exprssions would no longer be the same size. That would be a programmer error, and one that might be difficult to track down. By removing n3 from that expression, it is now guaranteed from the declarations x(n3,n1) and tmp3(n3) that the lhs and rhs are always the same size. That statement could also be written as x(:,i) = c1 * tmp3(:)**2 and the same goal could be accomplished. It is a matter of taste which to use. This is because n3 is a module variable that can be legally changed elsewhere, e.g. somewhere within CALC().