From 318f35554ef359cb5caafd5fcfcec4d30ce030d8 Mon Sep 17 00:00:00 2001
From: julie <julielangou@users.noreply.github.com>
Date: Fri, 28 Jan 2011 23:04:40 +0000
Subject: [PATCH] Correct bug0069 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)
---
 SRC/chetri2x.f | 18 +++++++++++++++---
 SRC/csytri2x.f | 18 +++++++++++++++---
 SRC/dsytri2x.f | 18 ++++++++++++++++--
 SRC/ssytri2x.f | 17 +++++++++++++++--
 SRC/zhetri2x.f | 18 +++++++++++++++---
 SRC/zsytri2x.f | 17 +++++++++++++++--
 6 files changed, 91 insertions(+), 15 deletions(-)

diff --git a/SRC/chetri2x.f b/SRC/chetri2x.f
index 566fede49..6a01d5d28 100644
--- a/SRC/chetri2x.f
+++ b/SRC/chetri2x.f
@@ -273,11 +273,17 @@
 *
         CALL CTRMM('L','U','C','U',NNB, NNB,
      $             CONE,A(CUT+1,CUT+1),LDA,WORK(U11+1,1),N+NB+1)
+*
+         DO I=1,NNB
+            DO J=I,NNB
+              A(CUT+I,CUT+J)=WORK(U11+I,J)
+            END DO
+         END DO
 *
 *          U01'invD*U01->A(CUT+I,CUT+J)
 *
          CALL CGEMM('C','N',NNB,NNB,CUT,CONE,A(1,CUT+1),LDA,
-     $              WORK,N+NB+1, ZERO, A(CUT+1,CUT+1), LDA)
+     $              WORK,N+NB+1, ZERO, WORK(U11+1,1), N+NB+1)
 *
 *        U11 =  U11T*invD1*U11 + U01'invD*U01
 *
@@ -438,13 +444,19 @@
 *
         CALL CTRMM('L',UPLO,'C','U',NNB, NNB,
      $             CONE,A(CUT+1,CUT+1),LDA,WORK(U11+1,1),N+NB+1)
-
+*
+         DO I=1,NNB
+            DO J=1,I
+              A(CUT+I,CUT+J)=WORK(U11+I,J)
+            END DO
+         END DO
+*
         IF ( (CUT+NNB) .LT. N ) THEN
 *
 *          L21T*invD2*L21->A(CUT+I,CUT+J)
 *
          CALL CGEMM('C','N',NNB,NNB,N-NNB-CUT,CONE,A(CUT+NNB+1,CUT+1)
-     $             ,LDA,WORK,N+NB+1, ZERO, A(CUT+1,CUT+1), LDA)
+     $             ,LDA,WORK,N+NB+1, ZERO, WORK(U11+1,1), N+NB+1)
        
 *
 *        L11 =  L11T*invD1*L11 + U01'invD*U01
diff --git a/SRC/csytri2x.f b/SRC/csytri2x.f
index 5812e6d14..a10661baa 100644
--- a/SRC/csytri2x.f
+++ b/SRC/csytri2x.f
@@ -271,11 +271,17 @@
 *
         CALL CTRMM('L','U','T','U',NNB, NNB,
      $             ONE,A(CUT+1,CUT+1),LDA,WORK(U11+1,1),N+NB+1)
+*
+         DO I=1,NNB
+            DO J=I,NNB
+              A(CUT+I,CUT+J)=WORK(U11+I,J)
+            END DO
+         END DO
 *
 *          U01'invD*U01->A(CUT+I,CUT+J)
 *
          CALL CGEMM('T','N',NNB,NNB,CUT,ONE,A(1,CUT+1),LDA,
-     $              WORK,N+NB+1, ZERO, A(CUT+1,CUT+1), LDA)
+     $              WORK,N+NB+1, ZERO, WORK(U11+1,1), N+NB+1)
 *
 *        U11 =  U11T*invD1*U11 + U01'invD*U01
 *
@@ -436,13 +442,19 @@
 *
         CALL CTRMM('L',UPLO,'T','U',NNB, NNB,
      $             ONE,A(CUT+1,CUT+1),LDA,WORK(U11+1,1),N+NB+1)
-
+*
+         DO I=1,NNB
+            DO J=1,I
+              A(CUT+I,CUT+J)=WORK(U11+I,J)
+            END DO
+         END DO
+*
         IF ( (CUT+NNB) .LT. N ) THEN
 *
 *          L21T*invD2*L21->A(CUT+I,CUT+J)
 *
          CALL CGEMM('T','N',NNB,NNB,N-NNB-CUT,ONE,A(CUT+NNB+1,CUT+1)
-     $             ,LDA,WORK,N+NB+1, ZERO, A(CUT+1,CUT+1), LDA)
+     $             ,LDA,WORK,N+NB+1, ZERO, WORK(U11+1,1), N+NB+1)
        
 *
 *        L11 =  L11T*invD1*L11 + U01'invD*U01
diff --git a/SRC/dsytri2x.f b/SRC/dsytri2x.f
index 742a81740..d0481f56f 100644
--- a/SRC/dsytri2x.f
+++ b/SRC/dsytri2x.f
@@ -270,11 +270,18 @@
 *
         CALL DTRMM('L','U','T','U',NNB, NNB,
      $             ONE,A(CUT+1,CUT+1),LDA,WORK(U11+1,1),N+NB+1)
+*
+         DO I=1,NNB
+            DO J=I,NNB
+              A(CUT+I,CUT+J)=WORK(U11+I,J)
+            END DO
+         END DO         
 *
 *          U01'invD*U01->A(CUT+I,CUT+J)
 *
          CALL DGEMM('T','N',NNB,NNB,CUT,ONE,A(1,CUT+1),LDA,
-     $              WORK,N+NB+1, ZERO, A(CUT+1,CUT+1), LDA)
+     $              WORK,N+NB+1, ZERO, WORK(U11+1,1), N+NB+1)
+        
 *
 *        U11 =  U11T*invD1*U11 + U01'invD*U01
 *
@@ -436,12 +443,19 @@
         CALL DTRMM('L',UPLO,'T','U',NNB, NNB,
      $             ONE,A(CUT+1,CUT+1),LDA,WORK(U11+1,1),N+NB+1)
 
+*
+         DO I=1,NNB
+            DO J=1,I
+              A(CUT+I,CUT+J)=WORK(U11+I,J)
+            END DO
+         END DO
+*
         IF ( (CUT+NNB) .LT. N ) THEN
 *
 *          L21T*invD2*L21->A(CUT+I,CUT+J)
 *
          CALL DGEMM('T','N',NNB,NNB,N-NNB-CUT,ONE,A(CUT+NNB+1,CUT+1)
-     $             ,LDA,WORK,N+NB+1, ZERO, A(CUT+1,CUT+1), LDA)
+     $             ,LDA,WORK,N+NB+1, ZERO, WORK(U11+1,1), N+NB+1)
        
 *
 *        L11 =  L11T*invD1*L11 + U01'invD*U01
diff --git a/SRC/ssytri2x.f b/SRC/ssytri2x.f
index 29168dbad..9c1630531 100644
--- a/SRC/ssytri2x.f
+++ b/SRC/ssytri2x.f
@@ -270,11 +270,17 @@
 *
         CALL STRMM('L','U','T','U',NNB, NNB,
      $             ONE,A(CUT+1,CUT+1),LDA,WORK(U11+1,1),N+NB+1)
+*
+         DO I=1,NNB
+            DO J=I,NNB
+              A(CUT+I,CUT+J)=WORK(U11+I,J)
+            END DO
+         END DO         
 *
 *          U01'invD*U01->A(CUT+I,CUT+J)
 *
          CALL SGEMM('T','N',NNB,NNB,CUT,ONE,A(1,CUT+1),LDA,
-     $              WORK,N+NB+1, ZERO, A(CUT+1,CUT+1), LDA)
+     $              WORK,N+NB+1, ZERO, WORK(U11+1,1), N+NB+1)
 *
 *        U11 =  U11T*invD1*U11 + U01'invD*U01
 *
@@ -436,12 +442,19 @@
         CALL STRMM('L',UPLO,'T','U',NNB, NNB,
      $             ONE,A(CUT+1,CUT+1),LDA,WORK(U11+1,1),N+NB+1)
 
+*
+         DO I=1,NNB
+            DO J=1,I
+              A(CUT+I,CUT+J)=WORK(U11+I,J)
+            END DO
+         END DO
+*
         IF ( (CUT+NNB) .LT. N ) THEN
 *
 *          L21T*invD2*L21->A(CUT+I,CUT+J)
 *
          CALL SGEMM('T','N',NNB,NNB,N-NNB-CUT,ONE,A(CUT+NNB+1,CUT+1)
-     $             ,LDA,WORK,N+NB+1, ZERO, A(CUT+1,CUT+1), LDA)
+     $             ,LDA,WORK,N+NB+1, ZERO, WORK(U11+1,1), N+NB+1)
        
 *
 *        L11 =  L11T*invD1*L11 + U01'invD*U01
diff --git a/SRC/zhetri2x.f b/SRC/zhetri2x.f
index 481ff6f06..71e780eb6 100644
--- a/SRC/zhetri2x.f
+++ b/SRC/zhetri2x.f
@@ -273,11 +273,17 @@
 *
         CALL ZTRMM('L','U','C','U',NNB, NNB,
      $             CONE,A(CUT+1,CUT+1),LDA,WORK(U11+1,1),N+NB+1)
+*
+         DO I=1,NNB
+            DO J=I,NNB
+              A(CUT+I,CUT+J)=WORK(U11+I,J)
+            END DO
+         END DO
 *
 *          U01'invD*U01->A(CUT+I,CUT+J)
 *
          CALL ZGEMM('C','N',NNB,NNB,CUT,CONE,A(1,CUT+1),LDA,
-     $              WORK,N+NB+1, ZERO, A(CUT+1,CUT+1), LDA)
+     $              WORK,N+NB+1, ZERO, WORK(U11+1,1), N+NB+1)
 *
 *        U11 =  U11T*invD1*U11 + U01'invD*U01
 *
@@ -438,13 +444,19 @@
 *
         CALL ZTRMM('L',UPLO,'C','U',NNB, NNB,
      $             CONE,A(CUT+1,CUT+1),LDA,WORK(U11+1,1),N+NB+1)
-
+*
+         DO I=1,NNB
+            DO J=1,I
+              A(CUT+I,CUT+J)=WORK(U11+I,J)
+            END DO
+         END DO
+*
         IF ( (CUT+NNB) .LT. N ) THEN
 *
 *          L21T*invD2*L21->A(CUT+I,CUT+J)
 *
          CALL ZGEMM('C','N',NNB,NNB,N-NNB-CUT,CONE,A(CUT+NNB+1,CUT+1)
-     $             ,LDA,WORK,N+NB+1, ZERO, A(CUT+1,CUT+1), LDA)
+     $             ,LDA,WORK,N+NB+1, ZERO, WORK(U11+1,1), N+NB+1)
        
 *
 *        L11 =  L11T*invD1*L11 + U01'invD*U01
diff --git a/SRC/zsytri2x.f b/SRC/zsytri2x.f
index 93e23d54f..41cffcb2b 100644
--- a/SRC/zsytri2x.f
+++ b/SRC/zsytri2x.f
@@ -271,11 +271,17 @@
 *
         CALL ZTRMM('L','U','T','U',NNB, NNB,
      $             ONE,A(CUT+1,CUT+1),LDA,WORK(U11+1,1),N+NB+1)
+*
+         DO I=1,NNB
+            DO J=I,NNB
+              A(CUT+I,CUT+J)=WORK(U11+I,J)
+            END DO
+         END DO         
 *
 *          U01'invD*U01->A(CUT+I,CUT+J)
 *
          CALL ZGEMM('T','N',NNB,NNB,CUT,ONE,A(1,CUT+1),LDA,
-     $              WORK,N+NB+1, ZERO, A(CUT+1,CUT+1), LDA)
+     $              WORK,N+NB+1, ZERO, WORK(U11+1,1), N+NB+1)
 *
 *        U11 =  U11T*invD1*U11 + U01'invD*U01
 *
@@ -436,13 +442,20 @@
 *
         CALL ZTRMM('L',UPLO,'T','U',NNB, NNB,
      $             ONE,A(CUT+1,CUT+1),LDA,WORK(U11+1,1),N+NB+1)
+*
+         DO I=1,NNB
+            DO J=1,I
+              A(CUT+I,CUT+J)=WORK(U11+I,J)
+            END DO
+         END DO
+*
 
         IF ( (CUT+NNB) .LT. N ) THEN
 *
 *          L21T*invD2*L21->A(CUT+I,CUT+J)
 *
          CALL ZGEMM('T','N',NNB,NNB,N-NNB-CUT,ONE,A(CUT+NNB+1,CUT+1)
-     $             ,LDA,WORK,N+NB+1, ZERO, A(CUT+1,CUT+1), LDA)
+     $             ,LDA,WORK,N+NB+1, ZERO, WORK(U11+1,1), N+NB+1)
        
 *
 *        L11 =  L11T*invD1*L11 + U01'invD*U01