forums.silverfrost.com

[qt]

Hi Paul,

here is another error Prof. Lassmann encounters:

The demonstration program Test_where leads to Error 112: Reference to undefined variable, array element or function result (/UNDEF).
It seems that this error is caused by the intrinsic function matmul, which might be the cause that some programs from Numerical Recipes do not run correctly.


Winapp

Program Test_where

Implicit None

Real , Dimension (1:2,1:2) :: u
Real , Dimension (1:2 ) :: b, w, tmp

u = +0.
u (1,1) = -1.
u (2,2) = -1.

b (1 ) = 1.
b (2 ) = 0.

w (1 ) = 1.
w (2 ) = 2.

Write (*,*) 'u = ', u
Write (*,*) 'b = ', b
Write (*,*) 'w = ', w

w = w - 1.

Write (*,*) 'w = ', w

! The following statement gives
! Error 112, Reference to undefined variable, array element
! or function result (/UNDEF)

where (w /= 0.0)
tmp=matmul(b,u)/w
elsewhere
tmp=0.0
end where

Write (*,*) 'tmp = ', tmp (1:2)

End Program Test_where

[JohnCampbell]

I'd check Prof. Lassmann 's version of ftn95, as my Ver 4.9.1 does not produce the errors you are finding.

However, I do get " tmp = 0.00000 6.611829E-37" , rather than tmp = 0.0 0.0, so I'm not sure what is happening there.

Why not calculate MATMUL before the WHERE ?

[mecej4]

[small]John Campbell wrote:[/small]
I'd check Prof. Lassmann 's version of ftn95, as my Ver 4.9.1 does not produce the errors you are finding.

However, I do get " tmp = 0.00000 6.611829E-37" , rather than tmp = 0.0 0.0, so I'm not sure what is happening there.


To reproduce the "undefined" error, compile with the /undef option.

That this is a compiler bug is reinforced by running the program (with the /undef option used) under SDBG. The line with the WHERE statement is where an array element with an undefined value is reported to be accessed. However, inspecting the variables shows that none of the elements of w are undefined, nor are any elements of the variables on the right side of the assignment on the next line (b, u,w).


[small]John Campbell wrote:[/small]
Why not calculate MATMUL before the WHERE ?

That makes the bug go away! Working around compiler bugs that way is feasible only for toy programs.


MECEJ4

[DrTip]

MECEJ4

Although I do agree this *might* be a compiler bug and if it is should be fixed.

However in my humble opinion the previous suggestion of a work around is the non toy solution. In my company we have some programs with 3 hour+ runtimes, which rapidly goes up with when they are chained together in batches. They are mostly legacy FTN77 code compiled to FTN95. In this commercial environment which I guess is a non toy environment you go with what works, and you test every single component for numerical validity otherwise it bites you later.

I guess I am just saying not many commercial users have the time to wait for a new Compiler to be released!

Personally I also thought that the work around was actually more readable with no additional variables required.

I came up with exactly the same work around before John posted when I played with the code.

the reason I say it *might* be a compiler bug is that, does the code actually make sense?

I don't have my FTN95 standard handy but I rarely use the where construct but when I do I see it as a shorthand for a do/if combination. Expanding the where out you must do the matmul() first, ie outside the loop.


tmp = matmul(a,b)

do i=1,2
if (abs(w(i) > 0.0)) then
tmp(i) = tmp(i)/w(i)
else
tmp(i) = 0.0
end if
end do

I might have also gone for somethign like this


where (abs(w) > 0.0)
w = 1/w
end where

tmp = w*(matmul(a,b))


Another possible source of problems is the use of of (w /= 0.0) rather (abs(w) > 0.0) my compiler will give a warning in an if statement about unpredictable results with equalities and reals. (this examples seems to fit the bill to me)











Carl

[JohnCampbell]

The probable cause of the error is that the compiler may calculate all values of MATMUL(a,b)/w, rather than only where W > 0.
A point I would also like to make is that in the 10 years that WHERE and FORALL have been available, I am yet to find a practical use for these constructs.
There was a spate of FORALL bugs a few years ago, which was the only time I tried to use FORALL. I gave up and used a DO loop and a few if statements. If the arrays become large, stack size limitations can make the temporary array values a problem for memory management, producing very inefficient solutions. There appears to be negligable performance improvement available, just the reduction in a few lines of DO constructs.
All WHERE appears to provide in this case is an inefficient way of programming that should be avoided.

[PaulLaidler]

Joerg

I have just got back after a brief holiday and I will check out these items as soon as I can.
There are probably bugs to fix here and I can confirm this after one or two tests.

On the general point, WHERE and FORALL are designed for high performance Fortran (parallel processing) and will almost certainly be slower than alternative constructions under FTN95 which is designed for a single processor PC. There are known bugs in this area and although they need to be fixed (so that teachers can demonstrate their use), they have a relatively low priority. For most practical purposes WHERE and FORALL should be avoided when using FTN95.

[DrTip]

How true.

However getting rid of a few do loops can be a good thing the matrix = scalar construct is a great invention for initialisation. The older code I maintain has massive 3 page subroutines given over to just to setting some multi rank arrays to 0. written in ftn95 these would be a lot more readable as A = 0.0 etc.

I am a big fan of readablity as I am picking up the pieces of people who weren't! With the newer constructs I decide whether the code is more or less readble with them and use them accordingly. Often they are less readable I would suggest the example here is an example of less readable.

Right I am goign to do some work



Carl

[PaulLaidler]

Joerg

I have looked at this and although technically the code reveals a bug, I do not think that it is a serious one.

From the point of view of a) efficiency b) clarity and c) writing compiler friendly code, it would be better if matmul was applied before the WHERE statement.

t1 = matmul(b,u)
where (w /= 0.0)
tmp=t1/w
elsewhere
tmp=0.0
end where

When viewed as parallel processing statements this is essentially a two stage process and should be programmed as such.

[KL]

Dear Paul,

I am the author of the "WHERE" problem. However, the program lines stem from the SVD algorithm of Numerical Recipes. When I found that this algorithm did not run, I extracted the error and came to the "WHERE" problem.

In my opinion it is not acceptable to recommend not to use an important Fortran 95 concept such as the WHERE construct. Other programs that use this construct should of course work with FTN95.

This is not the only problem that I found when I ran test problems from Numerical Recipes. I extracted several problems and Mr. Kuthe send these problems to the forum. Including the last 3 problems which I have sent to Mr. Kuthe I found 5 or 6 severe cases. In the case of the WHERE construct I would like to ask you to fix this problem.

Best regards

Klaus Lassmann



K. Lassmann

[PaulLaidler]

Klaus

We can certainly have another look at this problem however, the reasons I suggested that the issue was not a serious one are
a) the problem only occurs when using /UNDEF and
b) the matrix multiplication should be placed before the WHERE construct otherwise it is likely that the compiler will calculate the whole vector result for each scan of the construct even though it will only use one element from the vector each time. In other words the calculation will be highly repetitive unless the compiler is clever enough to take the matrix product out of the construct for you. What you need inside the construct is essentially a dot product to return a scalar result rather than a full matrix multiplication that returns a vector but I could not think a way to program this.

[PaulLaidler]

Klaus

We have done some more work on this and this is how things appear to stand at the moment.

If you make a simple call to MATMUL within a WHERE construct then you will get the right answer (with FTN95) but the resulting assembly code is inefficient.

If you include MATMUL within an assignment expression (within a WHERE construct), as in your code, then FTN95 will probably give the wrong answer (as well as producing it in an inefficient way!).

It may take a while to fix this bug because the construction is highly complex. When it is fixed, the problem of inefficiency will not necessarily go away.

It seems natural to me to say that MATMUL returns an array whatever the context and, if the programmer only wants a particular element from the result, then they should provide the corresponding code. Having said this, it may still be possible to make FTN95 work in the way that you expect.

[KL]

Dear Paul,

I am not sure whether you have received my more detailed message
on the WHERE problem. Here it is again, forget it if you had received it
already.

Dear Paul,

thank you very much for your answer. I have placed the matrix multiplication before
the where construct and it works. With the same modification, the SVD algorithm
from Numerical Recipes also runs. So we could consider the problem as being solved.

Nevertheless, I would like to ask you or your coworker to have another look to the
problem:

According to the Fortran 95 standard, the intrinsic function MATMUL is elemental.
(Of course all following arguments do not hold if this condition is not fulfilled).
If in the expression matmul (b, u) b has the shape m and u has shape m,n,
then the result has shape n. Element j has the value sum (b*u(:,j)) and it should
be possible
to devide this element j by w(j). So the algorithm of Numerical Recipes seems to
be correct. It is also efficient since it avoids the calculation
of element j when w(j) is zero. I guess that this is the reason why this algorithm
has been chosen in Numerical Recipes.

I have tested this problem with the little program "WHERE" below. Part 1 is the
FORTRAN 77 algorithm, that acts as a reference. Part 2 is
the modified and unmodified algorithm of Numerical
Recipes and part 3 is a simulation of the elemental treatment which should be
automatically included in Fortran 95.

It is interesting to note that the original algorithm from Numerical
Recipes seems to work correctly for n>m. So I guess that there is only a
minor error for n<=m.

It is also interesting to note that the modified algorithm leads to the
error "Access Violation" if m>n (This error has been reported already to the
forum and your answer to J. Kuthe was that this bug has been fixed).

To summarise: I am quite sure that the algorithm of
Numerical Recipes is correct and efficient (they normally think in
large system of equations) and that there might be a minor error in FTN95.

I have worked with FTN77 /95 since many years and I do hope that such a
discussion helps to increase the quality of FTN95.

Below
is the program. I am using FTN95 Version 4.7 and I have investigated
the problem with

Del *.lis
Del *.exe
Ftn95 Where.f95 /Checkmate /Link >> comp.lis
sdbg where.exe

!----------------------------------------------------------------------------------------
Winapp
Program Where

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

Implicit none

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

Integer , Parameter :: m = 5
Integer , Parameter :: n = 3

Real , Dimension (m,n) :: u
Real , Dimension (m ) :: b
Real , Dimension (n ) :: w
Real , Dimension (n ) :: tmp1, tmp2, tmp3
Real :: s
Integer :: i, j

u = 1.
b = 1.
w = 0.

w (1) = 1.

! ------------------------------------------------------------------
! FORTRAN 77 style

Do j = 1,n
s = 0.
if ( w (j) .ne. 0.) Then
do i = 1,m
s = s + u (i,j) * b (i)
end do
s = s / w (j)
end if
tmp1 (j) = s
End Do

Write (*,*) 'tmp1 = ', tmp1

! ------------------------------------------------------------------
! Modified and unmodified Fortran 95 style (see Numerical Recipes svbksb)

tmp2 = matmul (b,u)

where (w /= 0.0)
! Original statement from Numerical Recipes
! tmp2 = matmul (b,u)/w
tmp2 = tmp2/w
else where
tmp2=0.0
end where

! Write (*,*) 'tmp2 = ', tmp2 gives "Access Violation"
Write (*,*) 'tmp2 = ', tmp2 (1:n)

! ------------------------------------------------------------------
! "

[KL]

Dear Paul,

as I can see from the forum, my message was incomplete. I try it again:

thank you very much for your answer. I have placed the matrix multiplication before
the where construct and it works. With the same modification, the SVD algorithm
from Numerical Recipes also runs. So we could consider the problem as being solved.

Nevertheless, I would like to ask you or your coworker to have another look to the
problem:

According to the Fortran 95 standard, the intrinsic function MATMUL is elemental.
(Of course all following arguments do not hold if this condition is not fulfilled).
If in the expression matmul (b, u) b has the shape m and u has shape m,n,
then the result has shape n. Element j has the value sum (b*u(:,j)) and it should
be possible
to devide this element j by w(j). So the algorithm of Numerical Recipes seems to
be correct. It is also efficient since it avoids the calculation
of element j when w(j) is zero. I guess that this is the reason why this algorithm
has been chosen in Numerical Recipes.

I have tested this problem with the little program "WHERE" below. Part 1 is the
FORTRAN 77 algorithm, that acts as a reference. Part 2 is
the modified and unmodified algorithm of Numerical
Recipes and part 3 is a simulation of the elemental treatment which should be
automatically included in Fortran 95.

It is interesting to note that the original algorithm from Numerical
Recipes seems to work correctly for n>m. So I guess that there is only a
minor error for n<=m.

It is also interesting to note that the modified algorithm leads to the
error "Access Violation" if m>n (This error has been reported already to the
forum and your answer to J. Kuthe was that this bug has been fixed).

To summarise: I am quite sure that the algorithm of
Numerical Recipes is correct and efficient (they normally think in
large system of equations) and that there might be a minor error in FTN95.

I have worked with FTN77 /95 since many years and I do hope that such a
discussion helps to increase the quality of FTN95.

Below
is the program. I am using FTN95 Version 4.7 and I have investigated
the problem with

Del *.lis
Del *.exe
Ftn95 Where.f95 /Checkmate /Link >> comp.lis
sdbg where.exe

!----------------------------------------------------------------------------------------
Winapp
Program Where

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

Implicit none

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

Integer , Parameter :: m = 5
Integer , Parameter :: n = 3

Real , Dimension (m,n) :: u
Real , Dimension (m ) :: b
Real , Dimension (n ) :: w
Real , Dimension (n ) :: tmp1, tmp2, tmp3
Real :: s
Integer :: i, j

u = 1.
b = 1.
w = 0.

w (1) = 1.

! ------------------------------------------------------------------
! FORTRAN 77 style

Do j = 1,n
s = 0.
if ( w (j) .ne. 0.) Then
do i = 1,m
s = s + u (i,j) * b (i)
end do
s = s / w (j)
end if
tmp1 (j) = s
End Do

Write (*,*) 'tmp1 = ', tmp1

! ------------------------------------------------------------------
! Modified and unmodified Fortran 95 style (see Numerical Recipes svbksb)

tmp2 = matmul (b,u)

where (w /= 0.0)
! Original statement from Numerical Recipes
! tmp2 = matmul (b,u)/w
tmp2 = tmp2/w
else where
tmp2=0.0
end where

! Write (*,*) 'tmp2 = ', tmp2 gives "Access Violation"
Write (*,*) 'tmp2 = ', tmp2 (1:n)

! ------------------------------------------------------------------
! "Simulation of Elemental Style"


Do j = 1,n

If ( w (j) .ne. 0. ) then

[KL]

Del *.lis
Del *.exe
Ftn95 Where.f95 /Checkmate /Link >> comp.lis
sdbg where.exe

!----------------------------------------------------------------------------------------
Winapp
Program Where

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

Implicit none

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

Integer , Parameter :: m = 5
Integer , Parameter :: n = 3

Real , Dimension (m,n) :: u
Real , Dimension (m ) :: b
Real , Dimension (n ) :: w
Real , Dimension (n ) :: tmp1, tmp2, tmp3
Real :: s
Integer :: i, j

u = 1.
b = 1.
w = 0.

w (1) = 1.

! ------------------------------------------------------------------
! FORTRAN 77 style

Do j = 1,n
s = 0.
if ( w (j) .ne. 0.) Then
do i = 1,m
s = s + u (i,j) * b (i)
end do
s = s / w (j)
end if
tmp1 (j) = s
End Do

Write (*,*) 'tmp1 = ', tmp1

! ------------------------------------------------------------------
! Modified and unmodified Fortran 95 style (see Numerical Recipes svbksb)

tmp2 = matmul (b,u)

where (w /= 0.0)
! Original statement from Numerical Recipes
! tmp2 = matmul (b,u)/w
tmp2 = tmp2/w
else where
tmp2=0.0
end where

! Write (*,*) 'tmp2 = ', tmp2 gives "Access Violation"
Write (*,*) 'tmp2 = ', tmp2 (1:n)

! ------------------------------------------------------------------
! "Simulation of Elemental Style"


Do j = 1,n

If ( w (j) .ne. 0. ) then
tmp3 (j) = sum ( b * u (:,j) ) / w (j)
Else
tmp3 (j) = 0.
End If

End Do

Write (*,*) 'tmp3 = ', tmp3 ! gives "Access Violation"
Write (*,*) 'tmp3 = ', tmp3 (1:n)

End Program Where

[PaulLaidler]

Thank you. This problem is already in the pipeline.