?lan(gb/ge/gt/hs/sb/sf/sp/st/sy/tb/tp/tr) behavior with NaN

[zotkevich]

Dear lapackers,

I'd like to discuss behavior of the functions (see Subject) when input matrix A contains NaN at some point.

It seems natural that the functions should return NaN for such input, but it isn’t.

For instance ?lange could return real value even through input matrix contains NaN or NaNs when input parameter NORM is equal to 'M' or 'm', '1', 'O' or 'o', 'I' or 'i'.

It'd be nice if the functions from the list defined above behave more predictable for faulty input data. Functions could return NaN when input matrix contains NaNs.

[Julien Langou]

Hi,

This relates to the unpredictable behavior of the intrinsic MAX function in FORTRAN when one has a NaN as argument.
We can speak about this. Of course, we can preprocess the data to check for NaNs. (We have a DISNAN function so this should not be too hard.)
I think we could also use an INTRINSIC for this.

FYI, I found a nice conversation about this issue at:
http://objectmix.com/fortran/375103-wha ... t-nan.html
It's a mess indeed.

I am not sure what the MAX( NaN, 1) should be in the general case but
I think in our case, you are right, it makes more sense to have the
norm of a matrix with NaNs to be NaN.

Cheers,
Julien.

[zotkevich]

Do you have plans to add the preprocessing?

What do you think about 2 variants for instance for max(abs(A(i,j))):

1. Variant returns NaN if matrix contains NaNs and max(abs(A(i,j))) otherwise (modified variant ‘M’).
2. Variant returns max(abs(A(i,j))) passing NaNs.

[Julien Langou]

Hi Alex,

Let me try to involve more people in this thread.

( I just edited my post. ) I am not too much in favor of having the output of the matrix norm depends on a compilation flag.
I.e. having the behavior of the library depending on a compilation flag. I am not sure this is what you were proposing though.

However if it depends on an input argument, why not. We can offer other behaviors through other flags though. This can be
done at the Fortran level actually. There are plenty of characters different from M, F, I, O, 1, m, f, i, o. ( ?? )

If we need to decide for one behavior, I am in favor of having NaN returned for the norm (any norm) when the matrix has
some Nans.

Cheers,
Julien.

[Julien Langou]

And just as a follow up, I am posting Rodney's experience with the C++ std library. (Thanks Rodney.)

Code: Select all

   std::max(NaN,1.0) = NaN
   std::max(1.0,NaN) = 1.0


So actually, max( 1, NaN) can even return both: NaN and 1.
I did not think it about that one. (max being non-commutative.)

[zotkevich]

Hi Alex,

...If we need to decide for one behavior, I am in favor of having NaN returned for the norm (any norm) when the matrix has
some Nans.

Cheers,
Julien.

Hi, Julien,

At the moment it looks fine to have one behavior and I strongly support your suggestion.

Anyway a variant returning max(abs(A(i,j))) omitting NaNs could be added in the future if needed.

With reference to
std::max(NaN,1.0) = NaN
std::max(1.0,NaN) = 1.0
: yes I met similar issue. Result in this example depends on an implementation of MAX routine; logical operations involving NaN always return false.

[rodney]

I have modfied all of the 54 norm functions in LAPACK to consistently return a NaN when the input matrix contains a NaN. The new versions of the norm functions can be obtained from the trunk of the LAPACK svn repository:

svn co https://icl.cs.utk.edu/svn/lapack-dev/lapack/trunk

--Rodney

[zotkevich]

Hi Rodney,

Your solution is fine.

I have only one little remark about ?lassq: it looks like NaNs are being omitted in the routine.

WBR,
Alex.

[rodney]

Alex,

Thanks for pointing out the lassq routines. I have added a NaN check to those as well, and the new versions are available in the LAPACK SVN repository.

--Rodney

[Mark Hoemmen]

On a related note: Jack Poulson found that ZLARTG has an infinite loop related to NaN .EQ. 0. It's the line below "ELSE IF( SCALE.LE.SAFMN2 )" between 10 and 20. When G is NaN, then the loop goes forever.

[rodney]

Mark,

Thanks for pointing this out -- this has now been fixed for the next release, v3.5.

--Rodney

[asiakon]

Thanks for pointing this out -- this has now been fixed for the next release, v3.5.

[asiakon]

You cannot make another post so soon after your last.