Correct bug 0076 reported by Intel Team on lapack forum

Problem in ?(sy/he)tri2 was found.

For matrices with small N subroutine pass parameter NBMAX to ?(sy/he)tri2x which could be large than N.
In test for this functionality this is hide. Local test subroutine ilaenv.f pass proper values of NB. But if use ordinary ilaenv we will get mistake.

If we implemented

IF ( NBMAX .GE. N ) THEN
MINSIZE = N
ELSE
MINSIZE = (N+NBMAX+1)*(NBMAX+3)
END IF

and will call 

IF( NBMAX .GE. N ) THEN
CALL SSYTRI( UPLO, N, A, LDA, IPIV, WORK, INFO )
ELSE
CALL SSYTRI2X( UPLO, N, A, LDA, IPIV, WORK, NBMAX, INFO )
END IF

the problem will be solved.


Some other minor changes on workspace query.

10 files changed, 93 insertions, 37 deletions

diff --git a/SRC/chetri2.f b/SRC/chetri2.f
index f8c55ee9..6a496b5e 100644
--- a/SRC/chetri2.f
+++ b/SRC/chetri2.f
@@ -7,6 +7,8 @@
 *
 *  -- Written by Julie Langou of the Univ. of TN    --
 *
+* @generated c
+*
 *     .. Scalar Arguments ..
       CHARACTER          UPLO
       INTEGER            INFO, LDA, LWORK, N
@@ -19,10 +21,10 @@
 *  Purpose
 *  =======
 *
-*  CHETRI2 computes the inverse of a complex hermitian indefinite matrix
+*  CHETRI2 computes the inverse of a COMPLEX hermitian indefinite matrix
 *  A using the factorization A = U*D*U**T or A = L*D*L**T computed by
 *  CHETRF. CHETRI2 set the LEADING DIMENSION of the workspace
-*  before calling CHETRI2X that actually compute the inverse.
+*  before calling CHETRI2X that actually computes the inverse.
 *
 *  Arguments
 *  =========
@@ -94,7 +96,11 @@
       LQUERY = ( LWORK.EQ.-1 )
 *     Get blocksize
       NBMAX = ILAENV( 1, 'CHETRF', UPLO, N, -1, -1, -1 )
-      MINSIZE = (N+NBMAX+1)*(NBMAX+3)
+      IF ( NBMAX .GE. N ) THEN
+         MINSIZE = N
+      ELSE
+         MINSIZE = (N+NBMAX+1)*(NBMAX+3)
+      END IF
 *
       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
          INFO = -1
@@ -113,13 +119,17 @@
          CALL XERBLA( 'CHETRI2', -INFO )
          RETURN
       ELSE IF( LQUERY ) THEN
-         WORK(1)=(N+NBMAX+1)*(NBMAX+3)
+         WORK(1)=MINSIZE
          RETURN
       END IF
       IF( N.EQ.0 )
      $   RETURN
       
-      CALL CHETRI2X( UPLO, N, A, LDA, IPIV, WORK, NBMAX, INFO )
+      IF( NBMAX .GE. N ) THEN
+         CALL CHETRI( UPLO, N, A, LDA, IPIV, WORK, INFO )
+      ELSE
+         CALL CHETRI2X( UPLO, N, A, LDA, IPIV, WORK, NBMAX, INFO )
+      END IF
       RETURN
 *
 *     End of CHETRI2
diff --git a/SRC/csysv.f b/SRC/csysv.f
index 26c761b5..faa99004 100644
--- a/SRC/csysv.f
+++ b/SRC/csysv.f
@@ -145,8 +145,8 @@
          IF( N.EQ.0 ) THEN
             LWKOPT = 1
          ELSE
-            NB = ILAENV( 1, 'CSYTRF', UPLO, N, -1, -1, -1 )
-            LWKOPT = N*NB
+            CALL CSYTRF( UPLO, N, A, LDA, IPIV, WORK, -1, INFO )            
+            LWKOPT = WORK(1)
          END IF
          WORK( 1 ) = LWKOPT
       END IF
diff --git a/SRC/csytri2.f b/SRC/csytri2.f
index 9d74fe10..0bb7b50a 100644
--- a/SRC/csytri2.f
+++ b/SRC/csytri2.f
@@ -7,6 +7,7 @@
 *
 *  -- Written by Julie Langou of the Univ. of TN    --
 *
+* @generated c
 *     .. Scalar Arguments ..
       CHARACTER          UPLO
       INTEGER            INFO, LDA, LWORK, N
@@ -19,7 +20,7 @@
 *  Purpose
 *  =======
 *
-*  CSYTRI2 computes the inverse of a complex symmetric indefinite matrix
+*  CSYTRI2 computes the inverse of a COMPLEX hermitian indefinite matrix
 *  A using the factorization A = U*D*U**T or A = L*D*L**T computed by
 *  CSYTRF. CSYTRI2 sets the LEADING DIMENSION of the workspace
 *  before calling CSYTRI2X that actually computes the inverse.
@@ -94,7 +95,11 @@
       LQUERY = ( LWORK.EQ.-1 )
 *     Get blocksize
       NBMAX = ILAENV( 1, 'CSYTRF', UPLO, N, -1, -1, -1 )
-      MINSIZE = (N+NBMAX+1)*(NBMAX+3)
+      IF ( NBMAX .GE. N ) THEN
+         MINSIZE = N
+      ELSE
+         MINSIZE = (N+NBMAX+1)*(NBMAX+3)
+      END IF
 *
       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
          INFO = -1
@@ -113,13 +118,17 @@
          CALL XERBLA( 'CSYTRI2', -INFO )
          RETURN
       ELSE IF( LQUERY ) THEN
-         WORK(1)=(N+NBMAX+1)*(NBMAX+3)
+         WORK(1)=MINSIZE
          RETURN
       END IF
       IF( N.EQ.0 )
      $   RETURN
       
-      CALL CSYTRI2X( UPLO, N, A, LDA, IPIV, WORK, NBMAX, INFO )
+      IF( NBMAX .GE. N ) THEN
+         CALL CSYTRI( UPLO, N, A, LDA, IPIV, WORK, INFO )
+      ELSE
+         CALL CSYTRI2X( UPLO, N, A, LDA, IPIV, WORK, NBMAX, INFO )
+      END IF
       RETURN
 *
 *     End of CSYTRI2
diff --git a/SRC/dsysv.f b/SRC/dsysv.f
index 2c8a43e5..65ba671c 100644
--- a/SRC/dsysv.f
+++ b/SRC/dsysv.f
@@ -145,8 +145,8 @@
          IF( N.EQ.0 ) THEN
             LWKOPT = 1
          ELSE
-            NB = ILAENV( 1, 'DSYTRF', UPLO, N, -1, -1, -1 )
-            LWKOPT = N*NB
+            CALL DSYTRF( UPLO, N, A, LDA, IPIV, WORK, -1, INFO )            
+            LWKOPT = WORK(1)
          END IF
          WORK( 1 ) = LWKOPT
       END IF
diff --git a/SRC/dsytri2.f b/SRC/dsytri2.f
index ed32d543..6e1fa5bd 100644
--- a/SRC/dsytri2.f
+++ b/SRC/dsytri2.f
@@ -7,6 +7,7 @@
 *
 *  -- Written by Julie Langou of the Univ. of TN    --
 *
+* @generated d
 *     .. Scalar Arguments ..
       CHARACTER          UPLO
       INTEGER            INFO, LDA, LWORK, N
@@ -19,7 +20,7 @@
 *  Purpose
 *  =======
 *
-*  DSYTRI2 computes the inverse of a real symmetric indefinite matrix
+*  DSYTRI2 computes the inverse of a DOUBLE PRECISION hermitian indefinite matrix
 *  A using the factorization A = U*D*U**T or A = L*D*L**T computed by
 *  DSYTRF. DSYTRI2 sets the LEADING DIMENSION of the workspace
 *  before calling DSYTRI2X that actually computes the inverse.
@@ -94,7 +95,11 @@
       LQUERY = ( LWORK.EQ.-1 )
 *     Get blocksize
       NBMAX = ILAENV( 1, 'DSYTRF', UPLO, N, -1, -1, -1 )
-      MINSIZE = (N+NBMAX+1)*(NBMAX+3)
+      IF ( NBMAX .GE. N ) THEN
+         MINSIZE = N
+      ELSE
+         MINSIZE = (N+NBMAX+1)*(NBMAX+3)
+      END IF
 *
       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
          INFO = -1
@@ -113,13 +118,17 @@
          CALL XERBLA( 'DSYTRI2', -INFO )
          RETURN
       ELSE IF( LQUERY ) THEN
-         WORK(1)=(N+NBMAX+1)*(NBMAX+3)
+         WORK(1)=MINSIZE
          RETURN
       END IF
       IF( N.EQ.0 )
      $   RETURN
       
-      CALL DSYTRI2X( UPLO, N, A, LDA, IPIV, WORK, NBMAX, INFO )
+      IF( NBMAX .GE. N ) THEN
+         CALL DSYTRI( UPLO, N, A, LDA, IPIV, WORK, INFO )
+      ELSE
+         CALL DSYTRI2X( UPLO, N, A, LDA, IPIV, WORK, NBMAX, INFO )
+      END IF
       RETURN
 *
 *     End of DSYTRI2
diff --git a/SRC/ssysv.f b/SRC/ssysv.f
index 915fcae7..9743dd8a 100644
--- a/SRC/ssysv.f
+++ b/SRC/ssysv.f
@@ -145,8 +145,8 @@
          IF( N.EQ.0 ) THEN
             LWKOPT = 1
          ELSE
-            NB = ILAENV( 1, 'SSYTRF', UPLO, N, -1, -1, -1 )
-            LWKOPT = N*NB
+            CALL SSYTRF( UPLO, N, A, LDA, IPIV, WORK, -1, INFO )            
+            LWKOPT = WORK(1)
          END IF
          WORK( 1 ) = LWKOPT
       END IF
diff --git a/SRC/ssytri2.f b/SRC/ssytri2.f
index 21867f47..0d488bcc 100644
--- a/SRC/ssytri2.f
+++ b/SRC/ssytri2.f
@@ -7,19 +7,20 @@
 *
 *  -- Written by Julie Langou of the Univ. of TN    --
 *
+* @generated s
 *     .. Scalar Arguments ..
       CHARACTER          UPLO
       INTEGER            INFO, LDA, LWORK, N
 *     ..
 *     .. Array Arguments ..
       INTEGER            IPIV( * )
-      REAL              A( LDA, * ), WORK( * )
+      REAL               A( LDA, * ), WORK( * )
 *     ..
 *
 *  Purpose
 *  =======
 *
-*  SSYTRI2 computes the inverse of a real symmetric indefinite matrix
+*  SSYTRI2 computes the inverse of a REAL hermitian indefinite matrix
 *  A using the factorization A = U*D*U**T or A = L*D*L**T computed by
 *  SSYTRF. SSYTRI2 sets the LEADING DIMENSION of the workspace
 *  before calling SSYTRI2X that actually computes the inverse.
@@ -94,7 +95,11 @@
       LQUERY = ( LWORK.EQ.-1 )
 *     Get blocksize
       NBMAX = ILAENV( 1, 'SSYTRF', UPLO, N, -1, -1, -1 )
-      MINSIZE = (N+NBMAX+1)*(NBMAX+3)
+      IF ( NBMAX .GE. N ) THEN
+         MINSIZE = N
+      ELSE
+         MINSIZE = (N+NBMAX+1)*(NBMAX+3)
+      END IF
 *
       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
          INFO = -1
@@ -113,13 +118,17 @@
          CALL XERBLA( 'SSYTRI2', -INFO )
          RETURN
       ELSE IF( LQUERY ) THEN
-         WORK(1)=(N+NBMAX+1)*(NBMAX+3)
+         WORK(1)=MINSIZE
          RETURN
       END IF
       IF( N.EQ.0 )
      $   RETURN
       
-      CALL SSYTRI2X( UPLO, N, A, LDA, IPIV, WORK, NBMAX, INFO )
+      IF( NBMAX .GE. N ) THEN
+         CALL SSYTRI( UPLO, N, A, LDA, IPIV, WORK, INFO )
+      ELSE
+         CALL SSYTRI2X( UPLO, N, A, LDA, IPIV, WORK, NBMAX, INFO )
+      END IF
       RETURN
 *
 *     End of SSYTRI2
diff --git a/SRC/zhetri2.f b/SRC/zhetri2.f
index 7938b898..750dfaeb 100644
--- a/SRC/zhetri2.f
+++ b/SRC/zhetri2.f
@@ -7,6 +7,8 @@
 *
 *  -- Written by Julie Langou of the Univ. of TN    --
 *
+* @precisions normal z -> c
+*
 *     .. Scalar Arguments ..
       CHARACTER          UPLO
       INTEGER            INFO, LDA, LWORK, N
@@ -22,7 +24,7 @@
 *  ZHETRI2 computes the inverse of a COMPLEX*16 hermitian indefinite matrix
 *  A using the factorization A = U*D*U**T or A = L*D*L**T computed by
 *  ZHETRF. ZHETRI2 set the LEADING DIMENSION of the workspace
-*  before calling ZHETRI2X that actually compute the inverse.
+*  before calling ZHETRI2X that actually computes the inverse.
 *
 *  Arguments
 *  =========
@@ -94,7 +96,11 @@
       LQUERY = ( LWORK.EQ.-1 )
 *     Get blocksize
       NBMAX = ILAENV( 1, 'ZHETRF', UPLO, N, -1, -1, -1 )
-      MINSIZE = (N+NBMAX+1)*(NBMAX+3)
+      IF ( NBMAX .GE. N ) THEN
+         MINSIZE = N
+      ELSE
+         MINSIZE = (N+NBMAX+1)*(NBMAX+3)
+      END IF
 *
       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
          INFO = -1
@@ -113,13 +119,17 @@
          CALL XERBLA( 'ZHETRI2', -INFO )
          RETURN
       ELSE IF( LQUERY ) THEN
-         WORK(1)=(N+NBMAX+1)*(NBMAX+3)
+         WORK(1)=MINSIZE
          RETURN
       END IF
       IF( N.EQ.0 )
      $   RETURN
       
-      CALL ZHETRI2X( UPLO, N, A, LDA, IPIV, WORK, NBMAX, INFO )
+      IF( NBMAX .GE. N ) THEN
+         CALL ZHETRI( UPLO, N, A, LDA, IPIV, WORK, INFO )
+      ELSE
+         CALL ZHETRI2X( UPLO, N, A, LDA, IPIV, WORK, NBMAX, INFO )
+      END IF
       RETURN
 *
 *     End of ZHETRI2
diff --git a/SRC/zsysv.f b/SRC/zsysv.f
index f1f04a78..d5e81c8d 100644
--- a/SRC/zsysv.f
+++ b/SRC/zsysv.f
@@ -145,8 +145,8 @@
          IF( N.EQ.0 ) THEN
             LWKOPT = 1
          ELSE
-            NB = ILAENV( 1, 'ZSYTRF', UPLO, N, -1, -1, -1 )
-            LWKOPT = N*NB
+            CALL ZSYTRF( UPLO, N, A, LDA, IPIV, WORK, -1, INFO )            
+            LWKOPT = WORK(1)
          END IF
          WORK( 1 ) = LWKOPT
       END IF
diff --git a/SRC/zsytri2.f b/SRC/zsytri2.f
index 63fad25d..f4c06bae 100644
--- a/SRC/zsytri2.f
+++ b/SRC/zsytri2.f
@@ -7,19 +7,20 @@
 *
 *  -- Written by Julie Langou of the Univ. of TN    --
 *
+* @precisions normal z -> s d c
 *     .. Scalar Arguments ..
       CHARACTER          UPLO
       INTEGER            INFO, LDA, LWORK, N
 *     ..
 *     .. Array Arguments ..
       INTEGER            IPIV( * )
-      DOUBLE COMPLEX     A( LDA, * ), WORK( * )
+      COMPLEX*16         A( LDA, * ), WORK( * )
 *     ..
 *
 *  Purpose
 *  =======
 *
-*  ZSYTRI2 computes the inverse of a complex symmetric indefinite matrix
+*  ZSYTRI2 computes the inverse of a COMPLEX*16 hermitian indefinite matrix
 *  A using the factorization A = U*D*U**T or A = L*D*L**T computed by
 *  ZSYTRF. ZSYTRI2 sets the LEADING DIMENSION of the workspace
 *  before calling ZSYTRI2X that actually computes the inverse.
@@ -36,7 +37,7 @@
 *  N       (input) INTEGER
 *          The order of the matrix A.  N >= 0.
 *
-*  A       (input/output) DOUBLE COMPLEX array, dimension (LDA,N)
+*  A       (input/output) COMPLEX*16 array, dimension (LDA,N)
 *          On entry, the NB diagonal matrix D and the multipliers
 *          used to obtain the factor U or L as computed by ZSYTRF.
 *
@@ -54,7 +55,7 @@
 *          Details of the interchanges and the NB structure of D
 *          as determined by ZSYTRF.
 *
-*  WORK    (workspace) DOUBLE COMPLEX array, dimension (N+NB+1)*(NB+3)
+*  WORK    (workspace) COMPLEX*16 array, dimension (N+NB+1)*(NB+3)
 *
 *  LWORK   (input) INTEGER
 *          The dimension of the array WORK.
@@ -94,7 +95,11 @@
       LQUERY = ( LWORK.EQ.-1 )
 *     Get blocksize
       NBMAX = ILAENV( 1, 'ZSYTRF', UPLO, N, -1, -1, -1 )
-      MINSIZE = (N+NBMAX+1)*(NBMAX+3)
+      IF ( NBMAX .GE. N ) THEN
+         MINSIZE = N
+      ELSE
+         MINSIZE = (N+NBMAX+1)*(NBMAX+3)
+      END IF
 *
       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
          INFO = -1
@@ -113,13 +118,17 @@
          CALL XERBLA( 'ZSYTRI2', -INFO )
          RETURN
       ELSE IF( LQUERY ) THEN
-         WORK(1)=(N+NBMAX+1)*(NBMAX+3)
+         WORK(1)=MINSIZE
          RETURN
       END IF
       IF( N.EQ.0 )
      $   RETURN
       
-      CALL ZSYTRI2X( UPLO, N, A, LDA, IPIV, WORK, NBMAX, INFO )
+      IF( NBMAX .GE. N ) THEN
+         CALL ZSYTRI( UPLO, N, A, LDA, IPIV, WORK, INFO )
+      ELSE
+         CALL ZSYTRI2X( UPLO, N, A, LDA, IPIV, WORK, NBMAX, INFO )
+      END IF
       RETURN
 *
 *     End of ZSYTRI2