xLANGE underflows for finite 1x2 matrix

[christoph-conrads]

Description

SLANGE underflows for 1x2 matrices of the form [μ, 0], where μ is the smallest normal floating-point value. The overflow occurs when SCOMBSSQ sums the Euclidean norms of the individual columns. I expect SLANGE not to underflow for the given input. (Curiously enough, overflows in SLANGE were recently fixed.)

// Copyright 2021 Christoph Conrads

#include <cmath>
#include <cstddef>
#include <cstdint>
#include <cstdio>

#include <limits>

using Integer = std::int32_t;

extern "C" {
    float slange_(
        char* norm, Integer* m, Integer* n, const float* A, Integer* lda,
        float* work, std::size_t norm_len);
}



int main() {
    using Real = float;

    auto normid = 'F';
    auto m = Integer{1};
    auto n = Integer{2};
    auto min = std::numeric_limits<Real>::min();
    Real a[] = {min, 0};
    auto lda = m;
    auto work = std::numeric_limits<Real>::quiet_NaN();
    auto norm = slange_(&normid, &m, &n, a, &lda, &work, 1);

    if(!std::isfinite(norm) || norm != min) {
        std::fprintf(
            stderr, "expected norm %8.2e, got %e\n", min, norm);
    }
}

Output on Ubuntu 20.04:

/tmp/build.lapack $ g++ -Wextra -Wall -pedantic -std=c++11 ~/slange-underflow.cpp -Llib -llapack -g
/tmp/build.lapack $ env LD_LIBRARY_PATH="$PWD/lib" ./a.out 
expected norm 1.18e-38, got 0.000000e+00
/tmp/build.lapack $ git -C ~/lapack/ rev-parse HEAD
1c8d81aae95c6bb5d563d7e585f15af569a67ed0

Checklist

• I've included a minimal example to reproduce the issue
• I'd be willing to make a PR to solve this issue

[langou]

Thanks @christoph-conrads

Same issue with numpy.linalg.norm

import numpy

print(numpy.version.version)

# Frobenius-norm of the 1x2 matrix [[1e-170, 0]]
print(numpy.linalg.norm(numpy.array([[1e-170, 0]]),'fro'))

# 2-norm of the row vector [1e-170, 0]
print(numpy.linalg.norm(numpy.array([1e-170, 0]),2))

Output for me is

1.19.5
0.0
0.0

Matlab works just fine though

f = realmin('double'); norm([f,0],'fro')

output:

2.2251e-308

[weslleyspereira]

Hi. I did some tests on my laptop.

My system: Ubuntu 18.04.5 LTS (5.4.0-70-generic). GNU compilers, GCC version 7.5.0

LAPACK was compiled with

CC = gcc
CFLAGS = -O0
FC = gfortran
FFLAGS = -O0 -ggdb3 -frecursive -fimplicit-none -fcheck=all

1. Output using slange-underflow.cpp and current master (6ace55e) has the same behavior reported by @christoph-conrads:

weslley@lmpnote03:~/repositories/ballistic/565 - issue$ g++ -Wextra -Wall -pedantic -std=c++11 slange-underflow.cpp ~/repositories/lapack/liblapack.a -g -lgfortran
weslley@lmpnote03:~/repositories/ballistic/565 - issue$ ./a.out 
expected norm 1.18e-38, got 0.000000e+00
weslley@lmpnote03:~/repositories/ballistic/565 - issue$ git -C ~/repositories/lapack/ rev-parse HEAD
6ace55e9626bbf98a54e0c1199a171b4b645b572

2. Output using slange-underflow.cpp and master branch 77a0ceb (one commit before 9352192) works as expected:

weslley@lmpnote03:~/repositories/ballistic/565 - issue$ g++ -Wextra -Wall -pedantic -std=c++11 slange-underflow.cpp ~/repositories/lapack/liblapack.a -g -lgfortran -o slange-underflow-77a0ceb6c.out
weslley@lmpnote03:~/repositories/ballistic/565 - issue$ ./slange-underflow-77a0ceb6c.out 
weslley@lmpnote03:~/repositories/ballistic/565 - issue$ git -C ~/repositories/lapack/ rev-parse HEAD77a0ceb6c9a3c757e7039027f1f324e2811e6ee0

3. Output using snrm2-underflow.cpp (same code, but changes changes slange by snrm2) works as expected in the current master branch (6ace55e):

weslley@lmpnote03:~/repositories/ballistic/565 - issue$ g++ -Wextra -Wall -pedantic -std=c++11 snrm2-underflow.cpp ~/repositories/lapack/libblas.a -g -lgfortran -o snrm2-underflow-6ace55e.out
weslley@lmpnote03:~/repositories/ballistic/565 - issue$ ./snrm2-underflow-6ace55e.out
weslley@lmpnote03:~/repositories/ballistic/565 - issue$ git -C ~/repositories/lapack/ rev-parse HEAD
6ace55e9626bbf98a54e0c1199a171b4b645b572

I will prepare more tests.

[weslleyspereira]

Hi. I found the problem and proposed a solution (#569). The problem is that sCOMBSSQ updates the scaled sum V1 even if V2(2) is zero. Since a null sum can have whatever scaling value to be well represented in the machine, I thought the better place to apply the change would be in xCOMBSSQ.

[weslleyspereira]

Note:

In the master (6ace55e), the norm of a matrix like

[ min  0 ]
[ 0    0 ]

would also be wrongly computed. This is because:

1. SLANGE calls SLASSQ for each column.
2. The first column generates a non-null scaled sum with a low scaling factor
3. The second column generates a null scaled sum with scaling factor 1.
4. Then, SCOMBSSQ messes the sum up.

When we call xLANGE with a single column vector, SLASSQ is called once, so no problems with SCOMBSSQ in the old version.