lapackers

[julie <julie@icl.utk.edu>]

Hi,

As you may have seen, I committed the small fixes suggested by Ed to our LAPACK trunk.

For the other changes, I would prefer to have a discussion about them.

Thank you

Julie

Evaluation of LAPACK 3.5.0

Edward Anderson

January 26, 2015

I downloaded LAPACK 3.5.0 from netlib and installed it on a Red Hat Linux 6 system using the Intel 15.0, PGI 14.7, and Gnu 4.4.7 compilers.  Here are my comments on this release.

General comments

Since LAPACK development was started in the 1980s, the first releases had the look and feel of software developed at that time.   More recent developments in the Fortran standard and the OpenMP standard for expressing shared memory parallelism have made it possible to develop Fortran code that is more portable and easier to maintain.  Unfortunately, LAPACK 3.5.0 uses none of these features, even though the way forward has already been described and implemented in related packages such as LAPACK95 and LAPACK3E.  My previous suggestions for simplifying the algorithms in LAPACK and modernizing the software design are documented in “Performance Improvements to LAPACK for the Cray Scientific Library” (LAWN 126, with Mark Fahey, 1997) and “LAPACK3E:  A Fortran 90-enhanced version of LAPACK” (LAWN 158, 2002).

KIND and other constants

LAPACK 3.5.0 provides most subroutines in four different precisions, 32-bit and 64-bit real and 32-bit and 64-bit complex.  Declarations of REAL variables are changed to DOUBLE PRECISION and from COMPLEX to COMPLEX*16 for 64-bit, constant expressions are modified from 32-bit expressions to 64-bit expressions (for example, 1.0E+0 to 1.0D+0), intrinsic function calls are changed from 32-bit names to 64-bit names (for example, SQRT to DSQRT), and calls to auxiliary subroutines and functions are changed from their 32-bit equivalents to 64-bit equivalents (for example, SGEMM to DGEMM).  In the 1980s, these code transformations were automated by use of NAG’s Toolpack package, but for maintainers without access to Toolpack, keeping the 32-bit and 64-bit versions in sync became a challenge.

Innovations in the Fortran 90 standard (and subsequent updates) made most of these differences obsolete.  The KIND parameter allows declarations and expressions to be standardized between 32-bit and 64-bit, so it is only necessary to change the value of a single parameter.  All standard intrinsic functions were made generic, so that it is not necessary to use type-specific names such as DSQRT.  And modules were added to provide a standard way of including common definitions in multiple program units and for specifying generic interfaces for user-created subroutines.

LAPACK95, which should be familiar because it was developed with the participation of several LAPACK authors, used most of these features by parameterizing the KIND and providing interface modules as wrappers to the standard LAPACK routines.  However, it still provided new routines in four different precisions.  LAPACK3E went a step further, putting all the commonly used constants in a module, using generic interfaces for all subroutine calls, and providing include files with compile-time directives to allow the compiler to create 32-bit and 64-bit versions from a common source file.  (John Lewis, then at Boeing, urged me to do common source for real and complex too, but this proved to be unworkable.) The interface modules are a debugging aide because the compiler must match the number and type of arguments at compile time.

I still think that the LAPACK3E approach is the way to go.  If you give up the interface specifications, you give up compile time argument checking of your subroutine calls, and you lose the ability to fix a poor choice of arguments (for example, CPTSVX should have had an UPLO argument, and we can add one and still support LAPACK-style calls through the generic interface).  If you give up on the common source for 32-bit and 64-bit, you have to maintain a tool to convert from 32-bit to 64-bit and remember to regenerate the 64-bit version every time you make a change to the 32-bit version.  If you don’t put all your constants in a module, you have to pay the overhead of calling SLAMCH every time you want to use one, which leads to bad programming practices like storing the constants in a data area for use on subsequent calls.

Of course, there are tradeoffs to the LAPACK3E design.  The interface modules limit the amount of array reshaping you can do at call time, and the common source for 32-bit and 64-bit requires large include files of remapping specifications, which make it hard to extract a single module from LAPACK.

At the very least, you should parameterize the KIND and/or include the constants module.  Then 32-bit routines would have as their first lines after the SUBROUTINE or FUNCTION specification:

                  USE LA_CONSTANTS32

and 64-bit routines would begin with

                  USE LA_CONSTANTS

using the constants modules from LAPACK3E.  If you want to be really explicit, you could list the constants you intend to use, for example:

                  USE LA_CONSTANTS, ONLY: WP

This is what I am doing in the routines I have modified from LAPACK 3.5.0.

Reference BLAS

The Level 1 BLAS do not handle negative increments consistently.  Subroutines with two vectors interpret a negative increment as starting at the end of the vector and stepping back towards x(1), while subroutines with one vector are basically broken.  Rather than fix them, in the early 90s someone made things worse by having the reference BLAS return without an error condition if incx < 0.  In the Cray Scientific Library, and some other vendor libraries, negative increments are handled correctly, and some LAPACK routines require them to work that way in order to work as documented and to pass the unit tests I wrote.  It is necessary to replace xNRM2, xSCAL, xASUM, and IxAMAX if using the reference BLAS or any vendor library that follows the standard implied by the reference BLAS.

Software design

Besides converting from single precision to double precision automatically, Toolpack also formatted the LAPACK source code so that capitalization, spacing, indentation and variable declarations were consistent.  This is obviously not being done any more because the newest software is not consistent.  If Toolpack is no longer available, then an equivalent script should be written to perform these functions.

Structured programming was all the rage when I took my first programming class in 1980, and that style is reflected in the Level 2 and 3 BLAS, which contain no GOTO statements at all.  There were some GOTOs in LAPACK in subroutines that were largely unchanged from EISPACK.  However, LAPACK has regressed on this front as the latest release has 291 subroutines in 32-bit alone with GOTO statements in them.  Many of these could easily be eliminated with a better design or by using constructs such as EXIT.

Some of the newest routines have variable names that begin with LA_ or SLA_ or BLAS_.  The LA_ prefix is used to indicate a generic interface in LAPACK95 and should not be used for regular variables.  In fact, any short prefix following by an underscore should be avoided for regular variables to avoid conflicts with other libraries.  As a rule, I don’t like straying from the Fortran conventions of a leading I-N for integers and A-H, O-Z for floating point variables.  The BLAS and LAPACK have a style and it should be followed for consistency.

The new routines xGESVJ contain in the comments:

Unfortunately, some BLAS implementations compute SNRM2(M,A(1,p),1) as SQRT(SDOT(M,A(1,p),1,A(1,p),1)), which may cause the result to overflow for ||A(:,p)||_2 > SQRT(overflow_threshold), and to underflow for ||A(:,p)||_2 < SQRT(underflow_threshold).

Of course, we knew this in the 1980s, which is why we wrote SLASSQ, and why I have been insisting since the first LAPACK release that SLASSQ needs to be fixed.  Until we get our own routine right, we can’t whine about the BLAS.  Observations like this are inappropriate in the comments.

Iterative scaling

Several LAPACK routines do iterative scaling to try to get an input value into a desired range.  This will result in an infinite loop if the input value is NaN or infinity.  Iterative scaling with no stopping point is a very bad idea.  I have replaced these instances as I find them with either a single scaling or iterative scaling with a stopping point.

SLABAD

This function was introduced in LAPACK during the development phase to avoid testing problems near overflow and underflow on systems with an unscaled complex divide.  The “BAD” in the name was a political statement intended to disparage Cray, although in my view it is a warning of incomplete testing by the developers.  I fixed some minor scaling bugs and removed all calls to SLABAD before putting LAPACK in the Cray Scientific Library in 1992, and removed them again in LAPACK3E in 2002, but I no longer have access to a system with Cray arithmetic for testing.  This subroutine doesn’t do anything and could be removed any time.

LAPACK3E improvements

I have unit tests for subroutines that were changed significantly in LAPACK3E, including xLASSQ, xLARFG, xLARTG, and xGEBAL.  Recently I have revised xLASSQ to use multiple accumulators and a single pass through the vector as was done in James Blue’s original 1976 TOMS paper, acting on a suggestion from Damian McGuckin of Pacific Engineering Systems International.  We suspect that the Level 1 BLAS routine xNRM2 was not written by Blue because the algorithm in Blue’s original paper is cleaner.

The LAPACK version of SLASSQ fails 613 out of 3072 tests because it does not scale the vector values if they are smaller than the input value of SCL, even if the input value of SUMSQ is zero, and it does not handle negative increments as I expect.  A very reasonable way to initialize a sum of squares is with SCL = 1.0 and SUMSQ = 0.0, but LAPACK routines that call SLASSQ are forced to use SCL = 0.0 and SUMSQ = 1.0 or else they won’t be scaled properly.  In addition, my test program for SLASSQ assumes a negative increment means you want to traverse the vector in reverse order as is done in the Level 2 BLAS.  The LAPACK version starts with element 1 and goes into negative territory, which is unlikely to ever be what you want.

The LAPACK version of xLARFG does iterative scaling and will go into an infinite loop if BETA is negative infinity (although my tests do not use any exceptional values).  Both the real and complex xLARFG require negative increments to be fixed in xNRM2 and xSCAL in order to pass my tests.

The LAPACK version of SLARTG fails 40 out of 676 tests because the returned value of C is not always nonzero.  The complex case in LAPACK passes the test but does iterative scaling which makes it susceptible to going into an infinite loop.  This subroutine was supposed to be replaced with a new version described in a paper by Bindel and Demmel, but apparently the software never made it to release.

The LAPACK version of SGEBAL does about 3 times as much work as it has to and is susceptible to going into an infinite loop in the presence of NaNs.  For example, consider the following scaling loop for the array column sum C and row sum R, where CA is the maximum absolute value in the column and RA the maximum absolute value in the row.  The iterative scaling loop in LAPACK is

         G = R / SCLFAC

         F = ONE

160  CONTINUE

         IF( C.GE.G .OR. MAX( F, C, CA ).GE.SFMAX2 .OR.

     $       MIN( R, G, RA ).LE.SFMIN2 )GO TO 170

         F = F*SCLFAC

         C = C*SCLFAC

         CA = CA*SCLFAC

         R = R / SCLFAC

         G = G / SCLFAC

         RA = RA / SCLFAC

         GO TO 160

In LAPACK, SCLFAC = 2.  Now suppose some element in the column is NaN, so that both C and CA are NaN.  The first _expression_ in the IF statement will be false because C is NaN.  The second _expression_ is indeterminate in IEEE arithmetic, but if IEEE arithmetic were sane, it would also be NaN and the _expression_ would be false.  The third _expression_ need not involve any NaNs, and could be constructed to be false.  So we will do the scalings and return to the 160 statement with C and CA still NaN, and loop until R is very small, or forever if R is also NaN.

 

The same loop in LAPACK3E is implemented as

               F = R*SCL

               G = ONE

               DO WHILE( C < F .AND. CA < G*SFMAX .AND. RA*G >  SFMIN )

                  F = F*SCL2

                  G = G*SCL

               END DO

               F = ONE / G

In this context, SCL = 1/8 and SCL2 = SCL*SCL.  The SCL2 factor takes the place of scaling the row by SCL and the column by 1/SCL and follows the algorithm in EISPACK.  Now if C and CA were NaN, at least two of the three DO loop expressions would be false, and no scaling would be done.  This is the correct behavior, there is no point in scaling a matrix with NaNs in it.

My improvements to SGEBAL are described in LAWNs 126 and 158 and I attempted to contribute my version and new test code as early as 1995.  The toy problems in the old LAPACK test code are inadequate to exercise the code.

I replaced xLASSQ with the latest revision and xLARFG, xLARTG, xGEBAL, xRSCL, and xLAPY2 with better versions from LAPACK3E.  I also wrote a replacement for xNRM2 calling xLASSQ and added replacements for the Level 1 BLAS routines xSCAL, xASUM, and IxAMAX to LAPACK/SRC to fix the negative increments problem.  My new version of xLASSQ uses constants not found anywhere else; currently they are hard coded as parameters, if generally useful they could be added to the LA_CONSTANTS modules.

Other improvements

[CORRECTED 1523]
STGEX2 can call Level 3 BLAS and I did this in LAPACK3E.  I replaced STGEX2 with the LAPACK3E version (when I say STGEX2, I mean DTGEX2 too) and also fixed an incorrect leading dimension in DTGEX2.

[OK?]
When testing with the Gnu compiler, I had one test failure in the complex Hermitian generalized eigenvalue problem in which CHEGVX returned INFO = 1 because CSTEIN failed to converge.  The xSTEIN routines were extensively modified in LAWNs 126 and 158, although LAWN 126 describes a change in the orthogonalization procedure that Nick Higham pointed out amounted to replacing MGS with regular Gram-Schmidt.  I went back to regular Gram-Schmidt and lower level BLAS in LAPACK3E.  Substituting the LAPACK3E version, including auxiliary routines xLAEIF, xLAEIL, xLAEIU, and xCLCPY, solved this convergence failure (although I also changed the initialization of the random number seed from (1,1,1,1) to (0,0,0,1) which could have influenced the results. In case of failure, it might be worth trying a new value here before giving up.)

[OK?]

In SLALN2, SLASCL, and STRSYL, I modified the definition of SMLNUM to scale more conservatively, consistent with other uses of SMLNUM.

[OK?]

In SLASY2, I copied the fix for an illegal call to SCOPY (with INCX=0) from LAPACK3E.

[OK? How shall we signal failure of if not converging]

In SLASD8, I removed a call to XERBLA because it should only be used to signal failures in the error exits.

[OK?]

Using the Gnu gfortran compiler, version 4.4.7, I found that the BLAS routine crotg would not pass my unit tests of clartg because the cabs (or abs) intrinsic function did not scale enough.  I replaced cabs with my own subroutine scabs2 in crotg.  A separate unit test of the abs intrinsic shows that the bug only occurs for a few values near overflow.  LAPACK doesn’t call crotg, but this is also an issue for LAPACK because there are many places where the abs function is called with a complex argument.

[OK?]

When testing with the PGI compiler, I replaced some instances of ABS with CABS1 CABS2 and the others with SCABS2 in CLAHQR (or DZABS2 in ZLAHQR) to pass tests scaled near underflow.

[OK? – line commented; BIG=ONE/SFMIN was added in rev 891 by Zlatko – there is also a note in the code saying that you need to make sure SLAMCH is working properly on your machine ]    

I fixed a bad bug in SGESVJ caused by setting BIG to SLAMCH( ‘Overflow’ ), which is then used to construct a scaling constant that is later inverted.  There is no guarantee that SLAMCH( ‘Overflow’ ) is invertible.  The smallest number that you can safely invert is SLAMCH( ‘Safe minimum’ ), that is why it is called the “safe minimum”.  There is no inquiry function for ‘Safe maximum’, it is computed as the reciprocal of the safe minimum, so BIG should be ONE / SFMIN, which was in the code commented out.  It was a bad bug because the tests passed with the Intel compiler and failed with the PGI compiler (all the tests numbered 15 and 16 failed).  Even with this change, SGESVJ is still operating at the ragged edge of the floating point range and should probably be reworked to a safer range.  SGESVJ in general looks very rough, like a code under development.

LAPACK testing bug fixes

[CORRECTED 1519]
In TESTING/LIN, I fixed an old bug in the count of the number of tests in xdrvge and xdrvgb.

[CORRECTED 1520]
In TESTING/LIN, the test routines xchkqr kept track of the number of tests performed in the variable NT, which could be 2, 6, or 7.  Recently another test was added at the end and the loop that checks for failing tests was changed to run from 1 to NT to 1 to NTESTS, including tests that were never performed.  As a result, with THRESH = 0., the number of tests failing exceeds the number of tests performed (for example, with n = 1 as the only N value and THRESH = 0.0 you get “604 out of 514 tests failed to pass the threshold”.)  The right way to fix this would be to renumber the tests so that the new test is test #3 and the loop can run from 1 to NT again, but I just changed the count of tests performed from NT to NTESTS.  For consistency between schkqr and schkrq and other test routines in the QR family, it would have been better for the new tests to go in a separate test module.

[CORRECTED 1521]
In TESTING/LIN/schkrq.f and schklq.f (only those, not the other precisions), the loop to look for test failures runs from 1 to NTESTS, it should be 1 to NT.

[OK?]
In TESTING/EIG, I added my BAL/BAK test code to s/d/c/z chkee yet again and restored s/d/c/z ckbal, s/d/c/z blt01, and s/d/c/z blt02. 

[FILE From Ed come from 3.0]

Also in TESTING/EIG, I replaced illegal BLAS calls (using INCX=0) with inline code in s/d get34 and s/d/c/z get37.

[CORRECTED 1522 and 1524]

Also in TESTING/EIG, I fixed two bugs in s/d chkhs; the diff of schkhs is:

790c790

<             CALL SLACPY( ' ', N, N, U, LDU, UZ, LDU )

---

>             CALL SLACPY( ' ', N, N, U, LDU, UZ, LDA )

831c831

<      $                 ABS( WI1( J )-WI3( J ) ) )

---

>      $                 ABS( WR1( J )-WR3( J ) ) )

In further testing, I was able to produce a test case where Test 8 fails comparing the eigenvalues from xHSEQR(‘E’,’N’) with the eigenvalues from xHSEQR(‘S’,’V’).  The eigenvalues were the same, but in a different order.  It appears that xHSEQR does some rearrangement of the Schur matrix and so we should not expect the eigenvalues to be in the same order, although I am trying to confirm this with the authors.  In the meantime, I have modified xCHKEE and xCHKHS to compute the eigenvalues from xHSEQR(‘S’,’N’) in a separate array W2 (or WR2 and WI2) so that test 8 can compare the eigenvalues from xHSEQR(‘S’,’N’) with the eigenvalues from xHSEQR(‘S’,’V’).

Using Intel 15.0, I had to compile the following routines without optimization:

SRC:  slasd4, dlasd4

TESTING/EIG:  sget31, dget31, sget53, dget53

Comments

A new routine SLARFGP has appeared with an identical calling sequence and functionality to SLARFG, but with scaling to make beta positive.  It fails my unit tests of SLARFG because forcing beta to be positive means Re(tau) is no longer between 1 and 2.  I don’t understand why this routine is needed.

SORBDB6 has parameters REALONE and REALZERO which are the same as ONE and ZERO.  This sort of inconsistency would be fixed by using a module to define the constants as in LAPACK3E.

As in past LAPACK releases, xGGESX is not tested for any value of N greater than 2 and xGGEVX is not tested for any value of N greater than 5.  This makes me think these routines don’t really work.

There are many calls to xLARFG that are technically illegal because the same address is passed to two different arguments that would have to be specified as INOUT if you were to construct an interface.  For example, in SGELQ2 we have

      K = MIN( M, N )

      DO 10 I = 1, K

CALL SLARFG( N-I+1, A( I, I ), A( I, MIN( I+1, N ) ), LDA, TAU( I ) )

…

If N = 1 then K = 1, so on the one and only time through this loop the call becomes

CALL SLARFG( 1, A(1,1), A(1,1), LDA, TAU(1) )

The first A(1,1) goes to a scalar ALPHA while the second one goes to a vector X.  Both are labeled as input/output in the comments, however, if N = 1 as in this case, all SLARFG does is set TAU to zero and return.  So I think it should work, but if you had constructed an interface describing both ALPHA and X as INTENT(INOUT), the compiler would have a right to complain about it.  In other similar calls, this problem is avoided by doing something like

                  ALPHA = A(1,1)

                  CALL SLARFG( N-I+1, ALPHA, A( I, MIN( I+1, N ) ), LDA, TAU( I ) )

Alternatively, you could just skip the call and set TAU(I) to 0 if I = N.  It would mean an extra IF test in the inner loop, but you could make up for it by eliminating the MIN.

On Jan 29, 2015, at 10:03 AM, ICL <julie@icl.utk.edu> wrote:

Ed,

I am one of the keeper of the LAPACK source tree indeed.

I talked briefly with UC Denver about the changes you sent. They will have to be reviewed by utk, Denver and Berkeley - I am only committing bug fixes - any algo change needs to be discussed with the entire team
From your document, it seems that you have a patch for the routines you fixed.

If you sent it over to me, I would be happy to include your fixes in our source tree.

Regards,

Julie

On Jan 29, 2015, at 8:20 AM, Anderson, Edward <Anderson.Edward@epa.gov> wrote:

Julie,

 

So, are you the keeper of the LAPACK source tree?  I have some changes that should be submitted to someone.

 

   --Ed

 

Edward Anderson

Lockheed Martin, supporting the U.S. EPA

919-541-0299, anderson.edward@epa.gov

 

From: Jack Dongarra [mailto:dongarra@eecs.utk.edu]
Sent: Tuesday, January 27, 2015 4:32 PM
To: Julie Langou
Cc: Anderson, Edward
Subject: Fwd: Updated Evaluation of LAPACK

 

Julie,

Here are some comments from Ed on LAPACK.

Jack

**********************************************************
Prof. Jack Dongarra; Innovative Computing Laboratory; EECS Department;
1122 Volunteer Blvd; University of Tennessee; Knoxville TN 37996-3450;
+1-865-974-8295; dongarra@eecs.utk.edu; http://www.cs.utk.edu/~dongarra/

Begin forwarded message:

 

From: "Anderson, Edward" <Anderson.Edward@epa.gov>

To: Jack Dongarra <dongarra@cs.utk.edu>, "demmel@cs.berkeley.edu" <demmel@cs.berkeley.edu>

Subject: Updated Evaluation of LAPACK

Date: January 26, 2015 at 5:31:43 PM EST

 

Jack, Jim,

Here is an updated version of my "Evaluation of LAPACK" document.  Sorry, no change bars, for anyone who read the previous version.  If you'd like my software changes, please let me know who is managing the package these days.

  --Ed

Edward Anderson
Anderson.Edward@epa.gov
919-541-0299