bugs in ?GGHD3 : B,LDA-->B,LDB and zeroing of too large a submatrix

Quick and dirty bug fix.

Bug report and bug fix from Dan Kortschak.

See:
http://icl.cs.utk.edu/lapack-forum/viewtopic.php?t=4810
https://github.com/xianyi/OpenBLAS/issues/642

This is a quick bug fix, this only checks the triangular part of the trapezoid,
the rectangular part is not checked. Another commit should follow up soon.

Adding [un/or]csd2by1 routines to LAPACKE

Adding xlascl to LAPACKE

 Updating Copyright with current year

Adding xsytrf_rook, xsytrs_rook, xhetrf_rook and xhetrs_rook to LAPACKE

Adding xlahrd to DEPRECATED (problem reported by zerothi)

Fixing comments in [sd]lasd6

Issue : Parameter IDXQ of function SLASD6 is described as [out], should be described as [in,out].
Reported by dbakshee from INTEL on the forum on Sep 18th 2015
Link: http://icl.cs.utk.edu/lapack-forum/viewtopic.php?f=13&t=4799

Handle Deprecated routines in LAPACKE Makefile

Issues and fixes submitted by Christoph Conrads on Sept 14th by email
(1st commit)

-> File modified: DOCS/Doxyfile --
- Doxygen documentation
The Doxygen "project version" is still 3.4.2. The attached patch bumps
this to 3.6.0.
-> let's be ready!
- Doxygen documentation
The Doxygen configuration file uses multiple obsolete tags and warnings
are shown right after starting Doxygen. The warnings are easy to see if
you set "QUIET = YES" in DOCS/Doxyfile.
-> set "QUIET = YES"

-> File modified: SRC/*larfgp.f
It says "non-negatibe" in the documentation of xLARFGP.
-> Fixed typo

-> File modified: SRC/*orcsd2by1.f SRC/*uncsd2by1.f
The parameter M in the xORCSD2BY1 is not the number of columns of X and
there is no X12 (parameter P).
-> Fixed

-> File modified: LAPACKE/src/CMakeLists.txt
The C interfaces for the deprecated GSVD routines are still built with
cmake -DBUILD_DEPRECATED=OFF. This causes a linker error with cmake
-DBUILD_TESTING=ON as soon as the Fortran functions are not found in
liblapack.so.
-> Add a if(BUILD_DEPRECATED)

-> File modified: trunk/LAPACKE/include/lapacke.h
The C declaration (lapacke.h:11840-11843) for xGGSVD3 is in a block
of #defines for LAPACK 3.2.X (see lapacke.h:12089,12122).
-> moved them down

Bug Fix for bug0135
Bug reported by kortschak on Sep 10th
Link: http://icl.cs.utk.edu/lapack-forum/viewtopic.php?f=13&t=4793
Fix typo substituting O with 0

Fixing bug 133

Bug reported by the GONUM team member: btracey on LAPACK forum
Link: http://icl.cs.utk.edu/lapack-forum/viewtopic.php?f=13&t=4771
Description:  the number of columns of A depends on the Side, so the test should compare against m or n.

Fix by Julie.

Minor edits in the BLAS header to better describe the functionality.
See: http://icl.cs.utk.edu/lapack-forum/viewtopic.php?t=4786

Committing Nico's change about build system improvements for LAPACK

"One of the improvements with the patches is that SOVERSION and VERSION are now properly set and the proper symlinks are created:
```
liblapack.so -> liblapack.so.3
liblapack.so.3 -> liblapack.so.3.5.0
liblapack.so.3.5.0
```
Since BLAS is shipped with LAPACK and no separate version number is given for BLAS, I applied the same there, too."

Tested on Julie's Mac

TO DO: To test under Windows
====

Updated OSX RPATH settings
In response to CMake 3.0 generating warnings regarding policy CMP0042,
the OSX RPATH settings have been updated per recommendations found
in the CMake Wiki:
  http://www.cmake.org/Wiki/CMake_RPATH_handling#Mac_OS_X_and_the_RPATH

A is input/output for the subroutines ?sytrs2.f

Correcting Bug 129
See http://icl.cs.utk.edu/lapack-forum/viewtopic.php?f=13&t=4747

 Modify LAPACK_?syconv interface  - Work is used to return values

 (1) Remove the local work variable and its allocation and deallocation
 (2) Add work in LAPACK_?syconv interface
 (3) Update lapacke.h

correcting bug 128
See: http://icl.cs.utk.edu/lapack-forum/viewtopic.php?t=4746

When LAPACK_NAN_CHECK is enabled, in input of xLASET, we now allow if the matrix A contains NaN.

Since the goal of xLASET is to set values in A, we consider it OK if A contains NaN in input.

In other words, the rationale is that we do not consider A has an input/output. We consider A as an output only.

Note: We allow neither alpha nor beta to be a NaN.

(So we do check for NaN in alpha and in beta and return error if NaN.)

So this commit
(1) removes the NaN check on matrix A in LAPACKE wrappers: lapacke_?laset.c
(2) declares that A is OUTPUT in LAPACK subroutines: ?laset.f
(3) adds a thanks to Victor Kostin from Intel

Propagate comment added in rev 1568 to make.inc with gfortran compiler

changed variable names from LAM_MAX to SING_MAX and from LAM_MIN to SING_MIN to reflect the chnage in the algorithm for TEST4 for rook pivotig routines, since we now use singular values instead of eigenvalues in TEST4 to compute the norm for 2-by-2 marices in TEST 4. Also corrected the bug in for Hermitian case in TEST4, where during the constraction of Hermitian 2-by-2 matrices the operation of complex conjugation was missing. Changes made to: cchkhe_rook.f cchksy_rook.f dchksy_rook.f schksy_rook.f zchkhe_rook.f zchksy_rook.f

Moved deprecated SGVD routines into DEPRECATED directory.
Added a BUILD_DEPRECATED option for make and cmake, which build
deprecated routines into the LAPACK library (OFF by default).
This completes r1563.

TODO: document the BUILD_DEPRECATED option (i.e. update documentation).

Removed unused and obsolete statement functions.

Added LAPACKE wrappers for GSVD (complements r1572).

Added BLAS3 routines for generalised SVD.

TODO: LAPACKE wrappers.

Fixed BUG number 123, that is testing of 'rook' routines, specifically TEST 4 for UPLO='U' case that test bounds of block diagonal matrices. TEST 4 wasfailing for UPLO='U', since 2-by-2 was passed incorrectly to the -EVX routine.  Also, during the fix the  TEST 4 for both UPLO='U' and 'L' was modified. Instead of computing the ratio of eigenvalues for 2-by-2 matrices, using the -EVX routines, now we compute the ratio of singular values using *GESVD routines. Affected test routines: TESTING/LIN/schksy_rook.f TESTING/LIN/dchksy_rook.f TESTING/LIN/cchksy_rook.f TESTING/LIN/zchksy_rook.f TESTING/LIN/cchkhe_rook.f TESTING/LIN/zchkhe_rook.f

Removed declaration of unused variables.

Replaced archaic intrinsics AMIN1, DMIN1, MIN0, AMAX1, DMAX1, MAX0 by MIN/MAX.

adding a comment in make.inc.example saying the following:

#
#  Note: During a regular execution, LAPACK might create NaN and Inf
#  and handle these quantities appropriately. As a consequence, one
#  should not compile LAPACK with flags such as -ffpe-trap=overflow.
#

Julien.

Fix bugs reported by Vitorrio Zecca on August 10, 2015.

Fix bug reported by Vittorio Zecca.

Bux fix from Rodney James about implicit typing.

Some variables were undeclared, some were mistyped, and others were unused.

Move all but one deprecated routines to a new separate directory.
They are compiled and included in the liblapack.a archive when
the MAKEDEPRECATED option in the make.inc file is set.

TODO: Add the corresponding behavior to the cmake files.

NOTE: The xGEQPF routines are deprecated but still used by some
non-deprecated routines. They have not been moved to the DEPRECATED
directory.

Remove deprecated xGEQPF from the test suite.

WARNING: xGEQPF is still used in xGGSVP routines.

Remove all but one deprecated routines from the test suite.

NOTE: The tests are renumbered in xDRVLS and xCHKTZ.
TODO: remove deprecated xGEQPF, when it is replaced by xGEQP3 in xGGSVP.

missing EXTERNAL

Commit bug reported by Alexander Kobotov (INTEL) on LAPACK Forum
Link: http://icl.cs.utk.edu/lapack-forum/viewtopic.php?f=13&t=4762&p=11439#p11439
Recent commit r1533 introduced out of bound access which could lead to a crash.
It results scaling of a vector instead of a single number.
Fix is simple:

zhgeqz.f:669
- CALL ZSCAL( ILAST, SIGNBC, H( ILAST, ILAST ), 1 )
+ CALL ZSCAL( 1, SIGNBC, H( ILAST, ILAST ), 1 )

Similar for chgeqz on line 668

Revert commit 1474: files ?lansy.f are back as they were in revision 1473.

(minor) Consistency of comments between the four precisions.

Adding error test to each element of S must be positive. - Thank you Elena (Oracle)

Finish removing unnecessary code - TMP = DBLE(M) - see 1547 - Thank you Elena

Fix [CDSZ]QRT17.

This is a routine to test DGELS which solves solves full rank overdetermined ( linear least squares ) or underdetermined (minimum norm) linear systems.

The check is, as follows,

(1) If the system is consistent, (smallresidual,) we check that
|| A^T ( b – A * x ) || / (||A||*||b||*max(M,N,NRHS)*eps)
is small (say less than 30).

(2) If the system is not consistent, (large residual,) we check that
|| A^T ( b – A * x ) || / (||A||*|| b - A*x ||*max(M,N,NRHS)*eps)
is small (say less than 30).

DQRT17 was checking in the not consistent case (large residual) that
|| A^T ( b – A * x ) || / (||A||*||x||*max(M,N,NRHS)*eps)
is small. It is a mistake and we need to change the ||x|| by a ||r||.

This fix is consistent with what the documentation says and makes more sense.

Typo in the comments of [S,C,D,Z]LATM1.
The matrix D needs to be of size N.

These are two fixes in [S,C,D,Z]-LAGGE.

(1) The routine was writing ZEROs out of bound in some specific cases. This is
corrected. Thanks to Yi Meng (CU Denver) for finding the bug.

(2) Quick exit when KL=KU=0 in this case, the user wants a diagonal matrix,
there is probably no reason to call this routine if one simply wants a diagonal
matrix, but that said we still ought to return a correct answer. A quick exit
makes the trick. The routine was not working otherwise.

bug in [c,d,s,z]get02.f when calling the norm computation for A

typos

typos in these testing functions

IDIST is an INTEGER, change the comment

remove the COMPLEX*16 TMP from DLAHILB
this has nothing to do here, probably a bad copy-paste from Z to D

typos in the comment

This partially fixes bug 061 reported by Victor Liu.

Some compilers allocate local arrays on the heap when their size is
above a particular threshold.
This leads to wrong results when multiple threads call the same routine.

The bug fix consists in using a larger workspace, as proposed by Victor Liu
(Tue Nov 13, 2012).
Some routines still have large local arrays and cannot be fixed that way
because they have no workspace parameter: xGBTRF, xBPTRF, and xHSEQR.

Applying patch provided by Julien Schueller (Phimeca) on March 20th 2015 via Lapack mailing list to fix some warnings

Update version number and release date on rook routines - those routines were part of the 3.5.0 release - reported by Elena (Oracle)

Correct typo in variable name - reported by Elena (Oracle)

Propagate fix to bug111 to lapacke - minimum work space - Issue reported by Andrew James on LAPACK mailing list on April 30th 2015

Bug fix from Osni Marques, Beresford Parlett and Jim Demmel.
This fixes bugs 032 and 056.
Edits in cstein.f, dstein.f, sstein.f and zstein.f

From Osni on Monday, May 4:

As we discussed in our last conference call, I am attaching a new version of
_STEIN, which fixes bugs 032 and 056.  (The bugs were reported for DSTEIN, but
I have propagated the fix to the other versions.) Also, in the process of
testing the fix for those two bugs with more difficult cases, we stumbled upon
a matrix that led SSTEIN to return NaNs in the eigenvectors ... We fixed that
too and, in summary, these are the changes in _STEIN:

1) The assignment after the GO TO 60 needs to be GPIND = J1 (instead of GPIND =
B1).

2) In 'Normalize and scale the righthand side vector Pb' (after GO TO 100 or GO
TO 120) we have replaced _ASUM with I_AMAX.

Correct info testing for SYSVXX routines

in ?SYSVXX
 INFO = -10 needs to be changed to INFO = -11
The code in this section of the files is about to check and make sure all S
elements are positive. And S is argument 11 not 10.

in ?ERRVXX
 Add info testing for input parameter 10 and 11

Problem reported by Elena (Oracle)

Some arguments return the wrong Info number
For example LDX is not the 14th argument, but 15th.

Correct corresponding testing routines

**** Correct bug 125 ****

Change the maximum number of QR iterations (ITMAX) for xLAHQR from ITMAX=30 to
ITMAX = 30 * MAX( 10, NH ).

Bug reported by Yoshihide Okimi on July 16, 2014
https://icl.cs.utk.edu/lapack-forum/viewtopic.php?f=13&t=4566
and, independently, by Marco Caliari (Università di Verona) on February 11,
2015.

Involved in bug fix: Yoshihide Okimi, Marco Caliari (Università di Verona),
Meiyue Shao (LBL), Julien Langou, and Daniel Kressner (EPFL).

** Bug **

xLAHQR v3.5 does not converge on some matrices while xLAHQR v2 converges

** Note **

xLAHQR v2 is available through ARPACK so xLAHQR v2 is still in "wide" use

** Description **

xLAQR0 v3.5 and xLAQR4 v3.5 have ITMAX set to:
    ITMAX = MAX( 30, 2*KEXSH )*MAX( 10, ( IHI-ILO+1 ) )

xLAHQR v2 has ITMAX set to:
    ITMAX = 30 * NH

xLAHQR v3.5 has ITMAX set to:
    ITMAX = 30

Yoshihide Okimi reports that, in his application, (1) xLAHQR v3.5 does not
converge, (2) xLAHQR v2 converges, (3) changing ITMAX to 30*NH in xLAHQR v3.5
enables convergence.

Marco Caliari reports that, in his application, using xLAHQR v2.0 converges,
while xLAHQR v3.5 does not converge.

** Bug fix **

Set ITMAX in xLAHQR to
    ITMAX = 30 * MAX( 10, NH )

Email from Meiyue: I feel that setting ITMAX proportional to NH makes sense.
Perhaps we can change ITMAX to ITMAX*MAX( 10, NH ) so that it is also safe for
ill-conditioned tiny matrices.

**** Fix bug 124 ****
**** (I also fixed the documentation of the six subroutines ****

Bug reported by Victor_K on the forum.

See forum http://icl.cs.utk.edu/lapack-forum/viewtopic.php?f=13&t=4643

Bug fix by Julien Langou. Involved as well Nick Higham, Sven Hammarling and Craig Lucas.

** Bug **

xPSTF2 and xPSTRF returns NaNs when the matrix in input has N negative diagonal entries.

** Note **

A matrix with N negative diagonal entries is a highly non valid input, but, yet no reason to create some NaNs, better to have a cleaner solution.

** Description **

xPSTF2 and xPSTRF first scan the diagonal of the matrix to find the greatest element on the diagonal.

During this pass on the diagonal, xPSTF2/xPSTRF checks as well for ZERO or NAN on the diagonal and aborts if it finds a ZERO or a NAN on the diagonal.

*
*     Compute stopping value
*
      PVT = 1
      AJJ = A( PVT, PVT )
      DO I = 2, N
         IF( A( I, I ).GT.AJJ ) THEN
            PVT = I
            AJJ = A( PVT, PVT )
         END IF
      END DO
      IF( AJJ.EQ.ZERO.OR.DISNAN( AJJ ) ) THEN
         RANK = 0
         INFO = 1
         GO TO 170
      END IF

Since the input matrix is filled with a negative diagonal (and no ZERO and no NAN), xPSTF2/xPSTRF finds a maximum element (which is the largest negative elements on the diagonal), and then this (negative) element is used to compute DSTOP, the stopping value. But then we get a negative DSTOP,

      IF( TOL.LT.ZERO ) THEN
         DSTOP = N * DLAMCH( 'Epsilon' ) * AJJ
      ELSE
         DSTOP = TOL
      END IF

and then I stop even thinking on what would happen next. It seems that the code continues and does produce NAN in the first column. Why not. But I think we should catch this case and try to handle it in a different way. Maybe.

** Bug fix **

Replace
      IF( AJJ.EQ.ZERO.OR.DISNAN( AJJ ) ) THEN
With
      IF( AJJ.LE.ZERO.OR.DISNAN( AJJ ) ) THEN
So that if there is any nonpositive element or any NAN on the diagonal of A, we return right away with RANK=0 and INFO=1.

** Email from Nick Higham **

[ Julien wrote ] So that if there is any nonpositive element or any NAN on the diagonal of A, we return right away with RANK=0 and INFO=1?

[ Nick wrote ]

I think so.
This is justified by the documentation, which says

   INFO is INTEGER
   If INFO = -K, the K-th argument had an illegal value,
   = 0: algorithm completed successfully, and
   > 0: the matrix A is either rank deficient with computed rank
   as returned in RANK, or is indefinite.

We are in the indefinite case if any diagonal element is negative.
Actually, that is not quite true: we are in the "not positive definite"
case.  We can only say "indefinite" if there are both positive *and* negative
diagonal entries!  However, if that is going to be said then it would be
necessary to check that all other LAPACK routines are similarly precise in
their wording, so it's simplest, and not too misleading,
to stay with "indefinite".

Correcting Bug111 reported by Kevin Squire on LAPACK Mailing - August 15th 2013

From Kevin: We were following the recommendation for RWORK size calculation from
the ZGESDD header here ( http://netlib.org/lapack/complex16/zgesdd.f), and were
finding that when INFO=N, the size was insufficient, and generally led to a segmentation fault.

The current definition says that "If JOBZ = N, LRWORK >= 5*min(M,N)", but we found
that the minimum size needed to be at least 7*min(M,N). This was through experimentation
 and comparison with code from Numpy and Octave, so you may wish to verify this through the code.

JL: bug confirmed and reproduced - correction in comments
If JOBZ = N, LRWORK >= 7*min(M,N)

Correct issue reported by Elena (Oracle)
Add SIGMA and EPS in a list of arguments and describe them in the "Arguments" section

Correction reported by Elena (Oracle)

LAPACK 3.5.0 [cz]hgeqz.f compute eigenvalues differently depending on whether
eigenvectors are requested or not. Since the computation is different, eigenvalues
are not bit-by-bit identical whether eigenvectors are requested or not.

This fix enables to have [cz]hgeqz.f behaving more of the same whether eigenvectors
are requested or not.

Note: The fact that eigenvalues are bit-by-bit identical whether eigenvectors are
requested or not is checked in the LAPACK test suite.

Correct two ". GT." to " .GT.". Bug reported by email ny Willem Vermin "Minor
issue in dorbdb.f and sorbdb.f" on Friday, February 27, 2015 at 4:03 AM.

Fix wrong update when the right Schur vector is not computed.

Add back cggglm to Makefile...build goes fine

Add xGGHD3: blocked Hessenberg reduction, code from Daniel Kressner.
Add xGGES3 and xGGEV3: computation of the Schur form, the Schur vectors, and
   the generalized eigenvalues using the blocked Hessenberg reduction.

Fix typo in comments.

Fix typo in comments.

Fix typo.

Systematically return with INFO=1 when the root finder failed.
This corresponds to the documentation of [d/s]lasdx routines and fixes bug121
reported Justin Si.

Ed Anderson Bug Fix - Jan 26 2015

related to 1522
Correct call to [SD]chkhs to add new arguments.

Ed Anderson Improvement - Jan 26 2015

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.

Ed Anderson Bug Fix - Jan 26 2015

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’).

Ed Anderson Bug Fix - Jan 26 2015

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.

Ed Anderson Bug Fix - Jan 26 2015

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.

Ed Anderson Bug Fix - Jan 26 2015 - I fixed an old bug in the count of the number of tests in xdrvge and xdrvgb

Applied Jakub's CBLAS patch  sent on Jan 9th

Netlib CBLAS can be put on top of BLAS with 64-bit integers.
It can be accomplished using the flag "WeirdNEC" :)
It basically redefines int to long, but at the same time, it had multiple bugs.

Apply patch from Pavel Holoborodko.
See http://icl.cs.utk.edu/lapack-forum/viewtopic.php?f=13&t=4624

Please thanks to
  Pavel Holoborodko. Multiprecision Computing Toolbox for MATLAB. Advanpix LLC.
Below is a copy-paste from Pavel's post on the forum.

ZLARFT:
ZGEMM is not listed in EXTERNAL section (but extensively used in the code).

DLADIV:
This file has three different routines - DLADIV, DLADIV1, DLADIV2.

DLADIV2 has "if" in lower case:
         if( BR.NE.ZERO ) THEN

DGGBAL/BAK, ZGGBAL/BAK:
There are few places in these routines where DOUBLE PRECISION variables are implicitly converted to INTEGER.

For example:

         IR = LSCALE( I ) + SIGN( HALF, LSCALE( I ) )

         JC=RSCALE(I)+SIGN(HALF,RSCALE(I))

         K=RSCALE(I)

         K=LSCALE(I)

R/LSCALE(I) is double precision, IR, JC and K are integers.  For the sake of
correct semantic and F77 rules we have to use INT to show conversion
explicitly.  Also this would make conversion to other languages much easier.
Interestingly, only those routines are missing the INT.

 Fix typo reported by Elena (Oracle)

Bug Fix form Elena (Oracle) - Correct test of error exits

Bux fix from elena Ivanova (Oracle)  Rearrange If statements for not overwriting info value

Change INFO value when testing N - N is first parameter - INFO should be -1 - bug reported by Elena Ivanova (oracle)

Set default value for INFO - bug reported by Elena Ivanova (oracle) and David Hough

Comment Fix from Jim: Update of Purpose section of Jacobi SVD about accuracy

Comment fix from Jim: argument A in xsyconv.f should be in,out

From Jim: Routine description improvement - Adding guarantees that each R(i,i) >= 0 in comment and reference to lawn 203 in Furhter details section

correct typo while editing Makefile with VI

Fixing folder uppercase / lower case issue - Thank you Don

Commit fix from Osni Marquez for Bug 115 reported from Duncan Po (Mathworks)

Description of bug http://icl.utk.edu/lapack-forum/viewtopic.php?f=13&t=4391,
It turns out to be related to bug http://icl.cs.utk.edu/lapack-forum/viewtopic.php?f=5&t=4498 as well. Further, I suspect that the problem reported in an e-mail sent by Justin Weiguang Si to lapack@cs.utk.edu on Sep 22 is also related to those, although Justin has not yet replied to my request for the offending matrix...

The bug has been traced to [d/s]laed6, which computes the root closest to the origin of a secular equation and is used in the D&C tridiagonal eigensolver (DSTEDC) and D&C least squares solver (DGELSD). We have interacted with Ren-Cang Li (the original developer of [d/s]laed6) about possible fixes, and I am attaching a new version of [d/s]laed6. This version has been tested with our 'torture cases' (more on this below) in:

- single and double precision
- Intel Xeon Westmere with gnu (gfortran) and intel (ifort) compilers
- AMD MagnyCours with pgi, cray, intel and gnu compilers

The bug was related to a too stringent tolerance convergence criterion, line 390 in laed6. The fix is

390          IF( ( ABS( F ).LE.FOUR*EPS*ERRETM ) .OR.
391      $      ( (UBD-LBD).LE.FOUR*EPS*ABS(TAU) )  )

to replace

390          IF( ABS( F ).LE.EPS*ERRETM )

corrected bugs: 115 and 121

Fix problem in Info return value when checking input parameter

Definition of xLAED7( ICOMPQ, N, QSIZ, TLVLS, CURLVL, CURPBM, D, Q, [ etc..]

QSIZ is the third argument, not the fourth

       ELSE IF( ICOMPQ.EQ.1 .AND. QSIZ.LT.N ) THEN
-         INFO = -4
+         INFO = -3

Bug reported by Elena Ivanova (Oracle) directly to Julie on Oct 2nd 2014
Bug labelled as Bug122

Correction to bug120

Fix workspace definition in comment following bug report from Rasmus Munk Larsen, Google Knowledge & Research
See forum https://icl.cs.utk.edu/lapack-forum/viewtopic.php?f=13&t=4587

See revision 1349, and 649 for history

Typos in the comments of [s/d]tpmqrt.f reported by Andreas Noack Jensen on Sept. 15 2014.
See LAPACK forum: http://icl.cs.utk.edu/lapack-forum/viewtopic.php?f=13&t=4583

Clean CMAKE build for includes

Update make.inc's following the integration of CBLAS

Commit to integrate CMAKE and MAKEFILE build system for CBLAS
Also modify the way the CBLAS mangling was generated...now CBLAS follows LAPACKE model

Details on MAKEFILE build
  To compile CBLAS library from the LAPACK folder:
    - make cblaslib --> create the cblas lib
    - make cblas_testing --> compile BLAS Ref if necessary, build exe, then run exe

Details on CMAKE build
  Now you have an option to select CBLAS
  If CBLAS is check the Library is built and if BUILD_TESTING is enable, test and example are built and run

following
-->  modification of all make.inc
-->  more clean up
-->  debugging on other platforms

** PLEASE TRY and SEND FEEDBACK - Thanks - Julie **

See post on the forum from Pauli Virtanen.
http://icl.cs.utk.edu/lapack-forum/viewtopic.php?t=4560

Documentation of xSTEVR says that VL and VU are not referenced if RANGE = 'A'
or 'I', however the code does reference VL and VU for any value of RANGE:

      VLL = VL
      VUU = VU

So simply protect these two lines of code:

      IF( VALEIG ) THEN
         VLL = VL
         VUU = VU
      END IF

Julien.

Edited the code of xGBMV, xGEMV and xGEMM by removing the test for an entry of
B (or X) being zero to skip a loop so as to enable better NaN (or Inf)
propagation.

Taking DGEMM as an example, if you use notation:
C(I,J) = C(I,J) + alpha * A(I,L) * B(L,J),
the reference BLAS DGEMM is using the JLI version of matrix matrix multiply and
is checking, for each J (from 1 to N) and for each L (from 1 to K), whether
B(L,J) is zero (or not) to save (or not) the 2M following operations.

(See the "IF (B(L,J).NE.ZERO) THEN" in the code below.)

The snippets of code is as follows

              DO 90 J = 1,N
                  DO 80 L = 1,K
                      IF (B(L,J).NE.ZERO) THEN
                          TEMP = ALPHA*B(L,J)
                          DO 70 I = 1,M
                              C(I,J) = C(I,J) + TEMP*A(I,L)
   70                     CONTINUE
                      END IF
   80             CONTINUE
   90         CONTINUE

This induces some non NaN-propagation in a pretty ad-hoc way. For better NaN
propagation, this patch removes the above IF statement.

The snippet of code now becomes

              DO 90 J = 1,N
                  DO 80 L = 1,K
                      TEMP = ALPHA*B(L,J)
                      DO 70 I = 1,M
                          C(I,J) = C(I,J) + TEMP*A(I,L)
   70                 CONTINUE
   80             CONTINUE
   90         CONTINUE

This enables correct NaN propagation for this piece of code.

Rationale: BLAS does not correctly propagate all NaNs (and Infs). We still have
no NaN propagation where for example ALPHA=0, etc. The goal of this commit is
to have correct NaN propagation no matter what the entries of the input
matrices/vectors (A, B, C, X, etc.) are.  BLAS do not correctly propagate NaNs
and Infs based on some values of the scalars (ALPHA, BETA, etc.).

See below the email from Tom Callaway from RedHat, sent on July 9th to
lapack@cs.utk.edu.

Hello LAPACK people,

Martyn & Lejeczek (on CC) reported an issue to Fedora relating to R using our
system copy of BLAS (from LAPACK).

As noted in the R administration and Installation Manual, "R relies on ISO/IEC
60559 compliance of an external BLAS. This can be broken if for example the
code assumes that terms with a zero factor are always zero and do not need to
be computed - whereas x*0 can be NaN. This is checked in the test suite."

In the stock BLAS, DGBMV, DGEMM, and DGEMV fail this. R has been patching their
bundled BLAS to resolve this issue since 2010, but Fedora now uses the system
BLAS.

Attached is a patch (from upstream R) to fix this issue in the LAPACK BLAS.
Please consider applying it.

Thanks,

~tom

fix lapack_testing.py  for windows platforms - Fix provided by amran.illoul@ensam.eu on July 9th

Bug in comments on array size reported by a Sergey Kuznetsov from Intel.

The size of array RWORK in {c,z}dbsqr was given incorrectly.

Previous version for cbdsqr.f read
          RWORK is REAL array, dimension (2*N)
          if NCVT = NRU = NCC = 0, (max(1, 4*N-4)) otherwise

Now reads
          RWORK is REAL array, dimension (4*N)

bug found by Themos Tsikas fron NAG Ltd

LAPACK FORUM topic #4558

Hello

In zlanhf.f, first executable statements
Code: Select all
      IF( N.EQ.0 ) THEN
         ZLANHF = ZERO
         RETURN
      ELSE IF( N.EQ.1 ) THEN
         ZLANHF = ABS(A(0))
         RETURN
      END IF

When N=1, we should be ignoring any imaginary part (as we do for general diagonal elements later on) and so it should read
Code: Select all
      IF( N.EQ.0 ) THEN
         ZLANHF = ZERO
         RETURN
      ELSE IF( N.EQ.1 ) THEN
         ZLANHF = ABS(DBLE(A(0)))
         RETURN
      END IF

Themos Tsikas
NAG Ltd