Bug was sent by nmozarto on Jan 27th (see forum topic 2156)
Problem in new function ?SYTRI2 was found: the part of A below the diagonal is changed in the case UPLO='U' .
But in the description of arguments
If UPLO = 'U', the upper triangular part of the
inverse is formed and the part of A below the diagonal is not
referenced; if UPLO = 'L' the lower triangular part of the
inverse is formed and the part of A above the diagonal is
not referenced.
These elements zeroized after calling ?GEMM function in ?SYTRI2X.
CALL SGEMM('T','N',NNB,NNB,CUT,ONE,A(1,CUT+1),LDA,
$ WORK,N+NB+1, ZERO, A(CUT+1,CUT+1), LDA)
traite les matrices avec des lda.
Exemple:
> gfortran timer_dsytri2.c /Users/langou/Documents/lapack-dev-trunk.svn/lapack_LINUX.a -framework vecLib
> ./a.out -U -lda 1079 -n 100
DSYTRI2 info=0, check=1.891597e-14
Dispo dans ldlt.git. Marche bien. Upper ou Lower. LDA different de N ou pas. Le
check est
|| A . inv(A) - I || / || A ||
(Je pense que il faut diviser par || A || * || inv(A) ||, pas tres important.)
2) Du coup, dans ce commit, je rectifie tous les codes. Cela compile. Je n'ai
pas verifie les testings. Je n'ai pas ecrit de driver pour verifier. (Le driver
n'est que pour la version DOUBLE.)
JL
Je n'ai meme pas essaye de compiler en fait. En gros mettre LDA en parametre,
definir A comme etant A( LDA, N ) et pour l'increment du swap mettre LDA a la
place de N quand tu travailles sur les lignes. Donc c'est tout con.
Donc il faut ecrire un driver qui teste ca avec une LDA differente de N ... Et
appliquer les modifications aux autres routines CHE, CSY, ZHE, ZSY et SSY.
Sinon pourquoi dans DSYSWAPR, est ce que tu utilises DSWAP de temps en temps,
et de temps en temps tu le fais a la main? Je n'ai pas regarde les details.
JL.
Message from Zlatko:
The following changes are made in the
current version of the code:
1. in dgejsv, sgejsv:
- a typo in checking the parameters LSVEC.OR.LSVEC
has been changed to LSVEC.OR.RSVEC
- the length of WORK, LWORK, its description
and the minimal length for different JOBs
have been revised and corrected
- in a call to xGESVJ with WORK(N+1), the length
of the workspace is set to the correct value of
LWORK-N, instead of the incorrect LWORK
- a missing RETURN after a call to XERBLA
has been inserted
- In the case of zero matrix on input, IWORK(3)
is set to ZERO, to correspond to the description
of IWORK in other nontrivial cases
2. in all routines
- simple editing so that the single and
the corresponding double routines have matching
lines in the source codes.
2. Added CheckFortranIntSize.cmake to auto-determine the default integer size
3. Added more auto-detection features using integer size and compiler checks in FindBLAS.cmake
- Chuck Atkins (Kitware)
2. Added libblas to ESSL libs to cover missing functionality from ESSL
3. Renamed ESSL_6464, ESSL_SMP, and ESSL_SMP_6464 to ESSL6464, ESSLSMP, and ESSLSMP6464 respectively to mimic the actual library names.
4. Added copyright notice to CheckLAPACKCompilerFlags
- Chuck Atkins (Kitware)
2. Changed LOOK_FOR_OPTIMIZED_BLAS to USE_OPTIMIZED_BLAS per CMake convention
3. Changed BLAS_VENDORS to use library names instead of manufacturer names:
From: BLAS_VENDORS=AMD;IBM;INTEL;HP;SUN etc
To: BLAS_VENDORS=ACML;ESSL;MKL;VECLIB;PERFLIB
This seems to be a more intuitive interface
4. Added support for Intel MKL single dynamic library
- Chuck Atkins (Kitware)
Email from rencang:
"I have to clarify that the code does not do strict bisection. It uses bisection as the fall back when the computed approximation somehow gets out of the interval that the approximation should be in. Then one step of bisection is taken and then it goes back to the usual Newton-like method again (because that's faster). The interval is updated every step and thus shrinking. Xlasd4 was handled similarly.
Setting MAXIT=20 was too aggressive by my side as we know now. MAXIT=64 may be too conservative. Perhaps we should set MAXIT=40. If that breaks in the future, I'd really like to revisit the code. But if we'd like to play it safe, your suggested 64 would be the best bet."
julie