Skip to content

Bug in estimation of test error in [c,d,s,z]gebal with usage SNRM2/SCNRM2 #42

Open
@elivanova

Description

@elivanova

There is the following difference in /LAPACK/[s,d]gebal.f between 3.6.0 and 3.5.0:
In LAPACK-3.5.0, C and R are calculated as follows:

       DO 200 I = K, L
         C = ZERO
         R = ZERO

          DO 150 J = K, L
             IF( J.EQ.I )
     $         GO TO 150
             C = C + ABS( A( J, I ) )
             R = R + ABS( A( I, J ) )
  150  CONTINUE
. . . . . .
  200  CONTINUE

I.e. | Aii | is not added to C and R.

In LAPACK-3.6.0, C and R are calculated as follows:

        DO 200 I = K, L

           C = SNRM2( L-K+1, A( K, I ), 1 )
           R = SNRM2( L-K+1, A( I, K ), LDA )
. . . .
 200 CONTINUE

See /BLAS1/snrm2.f

_    SNRM2 := sqrt( x'_x )
Naturally, there will be another result, because SNRM2 is the euclidean norm and Aii elements are included in calculations.

The same situation is in the case of LAPACK/cgebal.f
In LAPACK-3.5.0, to calculate C and R (estimation of norm) they use a special function
*     .. Statement Function definitions ..
      CABS1( CDUM ) = ABS( REAL( CDUM ) ) + ABS( AIMAG( CDUM ) )
In LAPACK-3.5.0, C and R are calculated as follows:

      DO 200 I = K, L
         C = ZERO
         R = ZERO

         DO 150 J = K, L
            IF( J.EQ.I )
     $         GO TO 150
            C = C + CABS1( A( J, I ) )
            R = R + CABS1( A( I, J ) )
  150    CONTINUE
. . . . .
  200 CONTINUE
Also   CABS1( Aii ) are not added to C and R.
In LAPACK-3.6.0, C and R are calculated as follows:
      DO 200 I = K, L

         C = SCNRM2( L-K+1, A( K, I ), 1 )
         R = SCNRM2( L-K+1, A( I , K ), LDA )
. . . . .
200 CONTINUE

See /BLAS1/scnrm2.f.
*  SCNRM2 returns the euclidean norm of a vector via the function
*  name, so that

_     SCNRM2 := sqrt( conjg( x' )_x )

SCNRM2 should work correctly, because it is successfully used by a lot ot other routines.
Any case, Aii are used in calculations of C and R.

Why have you changed the working algorithm?
I think, balancing works correctly, because those changes were connected only with estimations of
test error.

The difference in results, for example, in cbal.out
In LAPACK-3.5.0:
.. test output of CGEBAL ..
 value of largest test error            =    0.000E+00
 example number where info is not zero  =    0
 example number where ILO or IHI wrong  =    0
 example number having largest error    =    0
 number of examples where info is not 0 =    0
 total number of examples tested        =   13

In LAPACK-3.6.0:
 .. test output of CGEBAL ..
 value of largest test error            =    0.875E+00
 example number where info is not zero  =    0
 example number where ILO or IHI wrong  =    0
 example number having largest error    =    6
 number of examples where info is not 0 =    0
 total number of examples tested        =   13

This is definitely because Aii are not calculated in one case (in LAPACK-3.5.0) and are calculated in another (in LAPACK-3.6.0).
I did the following experiment:
I used cgebal.f from LAPACK-3.5.0, but commented the IF condition:
      DO 200 I = K, L
         C = ZERO
         R = ZERO
*
         DO 150 J = K, L
*            IF( J.EQ.I )
*     $         GO TO 150
            C = C + CABS1( A( J, I ) )
            R = R + CABS1( A( I, J ) )
  150    CONTINUE
. . . . .
  200 CONTINUE

And i've received the same result, as in the case of LAPACK-3.5.0.

.......
.. test output of CGEBAL ..
 value of largest test error            =    0.875E+00
 example number where info is not zero  =    0
 example number where ILO or IHI wrong  =    0
 example number having largest error    =    6
 number of examples where info is not 0 =    0
 total number of examples tested        =   13

If that is a bug and get fixed in 3.6, should the test also be changed? Or is the "value of largest test error  =  0.875E+00" acceptable as passed?

Metadata

Metadata

Assignees

Type

No type

Projects

No projects

Milestone

No milestone

Relationships

None yet

Development

No branches or pull requests

Issue actions