Age | Commit message (Collapse) | Author | Files | Lines |
|
SRC/sgeev.f
SRC/dgeevx.f
SRC/dgeev.f
SRC/sgeevx.f
Reported by Konstantinos Kafoysas (Beta CAE Systems S.A., Greece) on Thu Aug
9th, 2012 through LAPACK mailing list.
> In the comments of dgeev function
>
> * The left eigenvector u(j) of A satisfies
> * u(j)**T * A = lambda(j) * u(j)**T
> * where u(j)**T denotes the transpose of u(j).
>
> u is supposed to satisfy u(j)**H * A = lambda(j) * u(j)**H
|
|
Routine affected: [CZ]HEGST et [CZ]HEGS2,
B is IN in comment, but in the code, B is actually changed, but after restored.
Here is the code portion that is to blame. [CZ]HEGS2
00203 IF( K.LT.N ) THEN
00204 CALL ZDSCAL( N-K, ONE / BKK, A( K, K+1 ), LDA )
00205 CT = -HALF*AKK
00206 CALL ZLACGV( N-K, A( K, K+1 ), LDA )
00207 CALL ZLACGV( N-K, B( K, K+1 ), LDB )
00208 CALL ZAXPY( N-K, CT, B( K, K+1 ), LDB, A( K, K+1 ),
00209 $ LDA )
00210 CALL ZHER2( UPLO, N-K, -CONE, A( K, K+1 ), LDA,
00211 $ B( K, K+1 ), LDB, A( K+1, K+1 ), LDA )
00212 CALL ZAXPY( N-K, CT, B( K, K+1 ), LDB, A( K, K+1 ),
00213 $ LDA )
00214 CALL ZLACGV( N-K, B( K, K+1 ), LDB )
00215 CALL ZTRSV( UPLO, 'Conjugate transpose', 'Non-unit',
00216 $ N-K, B( K+1, K+1 ), LDB, A( K, K+1 ),
00217 $ LDA )
00218 CALL ZLACGV( N-K, A( K, K+1 ), LDA )
|
|
All the routines from the SRC folder have been updated to integrate the current Doxygen layout.
Everything seems to be fine, all tests passed without problem.
|
|
|
|
previously, this was (usually) only the case for the 'F' norm
|
|
Followed recommendation, use the existing sorting code.
Report sent to LAPACK mailing list on June 26th 2012
From Joseph:
> There appears to be an inconsistency and possible bug in the dstemr
> implementation. When calculating the eigenvalues of a matrix, the
> returned eigenvalues are supposed to be returned in ascending order.
> Although this appears to be the case for N >= 3, it does not appear to
> be the case for N=2. I believe this happens because the dstemr routine
> has special cases for N=0,1, and 2, which immediately return after their
> computation. Because these cases return immediately, they do not call
> the sorting routines around line 723 (in LAPACK version 3.4.1). As
> such, a simple fix would be to have the N=2 case call this sorting code
> rather than returning.
|
|
|
|
Intel.
|
|
Confirmed and Corrected by Julie on May 21st 2012
====================================
EMAIL:
We would like to report a bug in {S,D}SYEVR and {C,Z}HEEVR routines. This bug causes writing two different data in the same address in IWORK, potentially producing wrong answers.
The bug description is provided in the bottom of the email. Please let us know, if you have any questions.
Thank you,
Line 543 in {S,D}SYEVR:
Line 717 on {C,Z}HEEVR:
=== Code Starts here==
INDISP = INDIBL + N
* IWORK(INDIFL:INDIFL+N-1) stores the indices of eigenvectors
* that corresponding to eigenvectors that fail to converge in
* DSTEIN. This information is discarded; if any fail, the driver
* returns INFO > 0.
INDIFL = INDISP + N
* INDIWO is the offset of the remaining integer workspace.
INDIWO = INDISP + N <- -- It is suspicious.
=====End of the Code=========
I think this should be INDIWO = INDIFL+N. Otherwise, subsequent {S,D}STEIN call takes the same address for IWORK(INDIWO) and INWROK(INDIFL).
CALL DSTEIN( N, WORK( INDD ), WORK( INDE ), M, W,
$ IWORK( INDIBL ), IWORK( INDISP ), Z, LDZ,
$ WORK( INDWK ), IWORK( INDIWO ), IWORK( INDIFL ),
$ INFO )
====================================
|
|
One more minor notice:
- In sgsvj0.f and sgsvk1.f: EPS and SFMIN become INTEGER in description whereas these are actually REAL.
|
|
on April 23th 2012
It seems that this one has been around forever
============
hi.
i found a subtle difference between dsyevd.f and ssyevd.f when using lapack 3.4.1.
dsyevd updates LOPT after calling dsytrd and dlacpy. but those codes are not visible in ssyevd.f, cheevd.f and zheevd.f.
pls refer to line 329, 344 in dsyevd.f.
i wonder whether those codes are necessary, because dsyevd has precalculated LOPT with at least 1+6*N+2*N**2.
if those codes must be there, why not for ssyevd?
thanks in advance.
Yao Tao
|
|
appearing on some platforms, most notably IBM/XLF, causing an infinite loop in xLARTG (which also should be fixed by enforcing a max iteration count when computing SCALE)
|
|
Just taking a quick look on the new release I noticed some minor issues there, which probably caused by a generation script:
- If you explore diff for files (c/z)gbrfsx, cgbsvx, cgbsvxx, cheequb, clanhf, zheequb you could find that number of datatype descriptions in documentation were changed to incorrect one (DOUBLE PRECISION instead of COMPLEX, or COMPLEX*16 instead of just COMPLEX).
+ cunbdb has CMPLX instead of COMPLEX for X12 parameter.
+ zggevx RWORK become REAL instead of correct DOUBLE PRECISION
- Number of files got following string "/ output)", which seems meaningless. Just grep for it.
|
|
|
|
|
|
|
|
Got a new bug for you! Mick Pont found this problem in DLAED6. The code in question is the 40 loop - when DSCALE(I)=TAU you get a divide by zero (rare in practice). This can cause some compilers to immediately stop, e.g. the Sun compiler.
Mick proposed solution is below:
DO 40 I = 1, 3
IF (DSCALE( I ).NE.TAU) THEN
TEMP = ONE / ( DSCALE( I )-TAU )
TEMP1 = ZSCALE( I )*TEMP
TEMP2 = TEMP1*TEMP
TEMP3 = TEMP2*TEMP
TEMP4 = TEMP1 / DSCALE( I )
FC = FC + TEMP4
ERRETM = ERRETM + ABS( TEMP4 )
DF = DF + TEMP2
DDF = DDF + TEMP3
ELSE
* On rare occasions dscale(i) can be exactly equal to
* tau, leading to division by zero. If no trap occurs,
* there is no problem; the quantities above all overflow
* and the test on abs(f) below sends you to the end
* with good results. If a trap occurs, though, the
* user program will stop. Avoid that happening by
* jumping directly out.
GO TO 60
END IF
40 CONTINUE
This seems to work OK in our testing.
|
|
(see http://netlib.org/lapack/Errata/vrac/lapack_known_issues.html)
|
|
|
|
Fix bug bug0088 reported by Mike Pont from NAG on the forum
(see http://icl.cs.utk.edu/lapack-forum/viewtopic.php?f=13&t=2893)
Actually there were a lot of problems regarding arguments checking.
I tried to correct most of them.
Apply the fix propose to all x[he/sy]rfsx.f routines
- Use IGNORE_CWISE as suggested to prevent use of unitialize variable PARAMS
I also appied the previous TESTING fix to all routines
INFO has new value in ZHERFSX (see description of INFO between ZHESVX and ZHESVXX)
This is set on line 634 (IF ( INFO .LE. N ) INFO = N + J) of zherfsx.f
And this is not handled by the testing LIN/zdrvhex.f
I just add .AND. INFO.LE.N at line 638 to avoid raising an error when INFO = N + J
At the moment, I would recommand a further look at those routines.
ZHE, ZSY led to 182 Tests failing to pass the threshold
and the same for complex
better than before, but still....
|
|
the case with the changes made to xLARFT
|
|
|
|
|
|
when the single (diagonal) element is negative
see: http://icl.cs.utk.edu/lapack-forum/archives/lapack/msg01195.html
|
|
|
|
|
|
DHGEQZ (QZ iteration failed)
* bug report by Hong Bo Peng Sandgren, on 03-19-2012.
* See link:http://icl.cs.utk.edu/lapack-forum/archives/lapack/msg01257.html[LAPACK Mailing list msg 01257]
I am doing some work with DGGEV. When I check the return msg and the actual code, I found something may be wrong. Here is part of comments in the header of DGGEV.F.
* INFO (output) INTEGER
* = 0: successful exit
* < 0: if INFO = -i, the i-th argument had an illegal value.
* = 1,...,N:
* The QZ iteration failed. No eigenvectors have been
* calculated, but ALPHAR(j), ALPHAI(j), and BETA(j)
* should be correct for j=INFO+1,...,N.
* > N: =N+1: other than QZ iteration failed in DHGEQZ.
* =N+2: error return from DTGEVC.
When INFO = 1...N, there is an error in DHGEQZ (QZ iteration failed). From the code, we can see it jumps to label 110 then set WORK(1) and return.
But in case of we scaled the matrix, we still need to undo scale for the output array ALPHAR, ALPHAI and BETA for those values j=INFO+1,...,N.
In DGEEVX, we can see that it jumps to label 50 in case of DHSEQR failure and then undo scale before return.
|
|
|
|
|
|
|
|
Now we can generate dll for LAPACK and LAPACKE directly for Mingw so that FORTRAN compiler is longer needed.
Because LAPACKE contains some routines from MATGEN (for PLASMA), LAPACKE will requires the tmglib library.
Add some LAPACK 3.4.0 routines were missing in the CMAKE LAPACK build.
|
|
|
|
|
|
|
|
shift, and set a d to zero if it is insignificant compared to the cumulative shift and the current shift is zero. This should guarantee that it always finishes in a linear number of iterations, without hurting accuracy at all.
|
|
|
|
documentation of cheequb.f (resp. zheequb.f)
|
|
(see http://icl.cs.utk.edu/lapack-forum/viewtopic.php?f=13&t=2893)
This is related to the LAPACK-XBLAS routine: zherfsx.f
Here is what I did:
- Introduce IINFO so that INFO is not overwritten
- Use IGNORE_CWISE as suggested to prevent use of unitialize variable PARAMS
But this did not fix the problem reported.
INFO has new value in ZHERFSX (see description of INFO between ZHESVX and ZHESVXX)
This is set on line 634 (IF ( INFO .LE. N ) INFO = N + J) of zherfsx.f
And this is not handled by the testing LIN/zdrvhex.f
I just add .AND. INFO.LE.N at line 638 to avoid raising an error when INFO = N + J
Please send feedback as I am not sure this is the best way to fix the issue.
I will commit other precision once fix approved.
Thanks
Julie
|
|
Routien Name was modified before.
|
|
|
|
|
|
ssytf2_rook.o, slasyf_rook.o, ssytrf_rook.o, ssytrs_rook.o, ssytri_rook.o, ssycon_rook.o, ssysv_rook.o to SLASRC entry;
dsytf2_rook.o, dlasyf_rook.o, dsytrf_rook.o, dsytrs_rook.o, dsytri_rook.o, dsycon_rook.o, dsysv_rook.o to DLASRC entry;
csytf2_rook.o, clasyf_rook.o, csytrf_rook.o, csytrs_rook.o, csytri_rook.o, csycon_rook.o, csysv_rook.o to CLASRC entry;
zsytf2_rook.o, zlasyf_rook.o, zsytrf_rook.o, zsytrs_rook.o, zsytri_rook.o, zsycon_rook.o, zsysv_rook.o to ZLASRC entry;
|
|
algorithm for for symmetric indefinite matrices: zlasyf_rook.f zsytri_rook.f zsytrs_rook.f zsytrf_rook.f zsysv_rook.f zsycon_rook.f zsytf2_rook.f
|
|
*SYTRF by *SYTRF_ROOK
|
|
N+1 ), WORK, IWORK, AINVNM, KASE, ISAVE ) \n to CLACN2( N, WORK( N+1 ), WORK, AINVNM, KASE, ISAVE )
|
|
clasyf_rook.f and csytf2_rook.f
|
|
for symmetric indefinite matrices: slasyf_rook.f ssytri_rook.f ssytrs_rook.f ssytrf_rook.f ssysv_rook.f ssycon_rook.f ssytf2_rook.f
|
|
|
|
|
|
|