Diff for /rpl/lapack/lapack/zsytri2x.f between versions 1.2 and 1.3

version 1.2, 2010/12/21 13:53:56 version 1.3, 2011/07/22 07:38:20
Line 1 Line 1
       SUBROUTINE ZSYTRI2X( UPLO, N, A, LDA, IPIV, WORK, NB, INFO )        SUBROUTINE ZSYTRI2X( UPLO, N, A, LDA, IPIV, WORK, NB, INFO )
 *  *
 *  -- LAPACK routine (version 3.3.0) --  *  -- LAPACK routine (version 3.3.1) --
 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --  *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--  *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 *     November 2010  *  -- April 2011                                                      --
 *  *
 *  -- Written by Julie Langou of the Univ. of TN    --  *  -- Written by Julie Langou of the Univ. of TN    --
 *  *
Line 154 Line 154
   
       IF( UPPER ) THEN        IF( UPPER ) THEN
 *  *
 *        invA = P * inv(U')*inv(D)*inv(U)*P'.  *        invA = P * inv(U**T)*inv(D)*inv(U)*P**T.
 *  *
         CALL ZTRTRI( UPLO, 'U', N, A, LDA, INFO )          CALL ZTRTRI( UPLO, 'U', N, A, LDA, INFO )
 *  *
Line 182 Line 182
          END IF           END IF
         END DO          END DO
 *  *
 *       inv(U') = (inv(U))'  *       inv(U**T) = (inv(U))**T
 *  *
 *       inv(U')*inv(D)*inv(U)  *       inv(U**T)*inv(D)*inv(U)
 *  *
         CUT=N          CUT=N
         DO WHILE (CUT .GT. 0)          DO WHILE (CUT .GT. 0)
Line 267 Line 267
              END IF               END IF
            END DO             END DO
 *      *    
 *       U11T*invD1*U11->U11  *       U11**T*invD1*U11->U11
 *  *
         CALL ZTRMM('L','U','T','U',NNB, NNB,          CALL ZTRMM('L','U','T','U',NNB, NNB,
      $             ONE,A(CUT+1,CUT+1),LDA,WORK(U11+1,1),N+NB+1)       $             ONE,A(CUT+1,CUT+1),LDA,WORK(U11+1,1),N+NB+1)
 *  *
 *          U01'invD*U01->A(CUT+I,CUT+J)           DO I=1,NNB
               DO J=I,NNB
                 A(CUT+I,CUT+J)=WORK(U11+I,J)
               END DO
            END DO         
   *
   *          U01**T*invD*U01->A(CUT+I,CUT+J)
 *  *
          CALL ZGEMM('T','N',NNB,NNB,CUT,ONE,A(1,CUT+1),LDA,           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  *        U11 =  U11**T*invD1*U11 + U01**T*invD*U01
 *  *
          DO I=1,NNB           DO I=1,NNB
             DO J=I,NNB              DO J=I,NNB
Line 285 Line 291
             END DO              END DO
          END DO           END DO
 *  *
 *        U01 =  U00T*invD0*U01  *        U01 =  U00**T*invD0*U01
 *  *
          CALL ZTRMM('L',UPLO,'T','U',CUT, NNB,           CALL ZTRMM('L',UPLO,'T','U',CUT, NNB,
      $             ONE,A,LDA,WORK,N+NB+1)       $             ONE,A,LDA,WORK,N+NB+1)
Line 303 Line 309
 *  *
        END DO         END DO
 *  *
 *        Apply PERMUTATIONS P and P': P * inv(U')*inv(D)*inv(U) *P'  *        Apply PERMUTATIONS P and P**T: P * inv(U**T)*inv(D)*inv(U) *P**T
 *    *  
             I=1              I=1
             DO WHILE ( I .LE. N )              DO WHILE ( I .LE. N )
                IF( IPIV(I) .GT. 0 ) THEN                 IF( IPIV(I) .GT. 0 ) THEN
                   IP=IPIV(I)                    IP=IPIV(I)
                  IF (I .LT. IP) CALL ZSYSWAPR( UPLO, N, A, I ,IP )                   IF (I .LT. IP) CALL ZSYSWAPR( UPLO, N, A, LDA, I ,IP )
                  IF (I .GT. IP) CALL ZSYSWAPR( UPLO, N, A, IP ,I )                   IF (I .GT. IP) CALL ZSYSWAPR( UPLO, N, A, LDA, IP ,I )
                ELSE                 ELSE
                  IP=-IPIV(I)                   IP=-IPIV(I)
                  I=I+1                   I=I+1
                  IF ( (I-1) .LT. IP)                    IF ( (I-1) .LT. IP) 
      $                  CALL ZSYSWAPR( UPLO, N, A, I-1 ,IP )       $                  CALL ZSYSWAPR( UPLO, N, A, LDA, I-1 ,IP )
                  IF ( (I-1) .GT. IP)                    IF ( (I-1) .GT. IP) 
      $                  CALL ZSYSWAPR( UPLO, N, A, IP ,I-1 )       $                  CALL ZSYSWAPR( UPLO, N, A, LDA, IP ,I-1 )
               ENDIF                ENDIF
                I=I+1                 I=I+1
             END DO              END DO
Line 325 Line 331
 *  *
 *        LOWER...  *        LOWER...
 *  *
 *        invA = P * inv(U')*inv(D)*inv(U)*P'.  *        invA = P * inv(U**T)*inv(D)*inv(U)*P**T.
 *  *
          CALL ZTRTRI( UPLO, 'U', N, A, LDA, INFO )           CALL ZTRTRI( UPLO, 'U', N, A, LDA, INFO )
 *  *
Line 353 Line 359
          END IF           END IF
         END DO          END DO
 *  *
 *       inv(U') = (inv(U))'  *       inv(U**T) = (inv(U))**T
 *  *
 *       inv(U')*inv(D)*inv(U)  *       inv(U**T)*inv(D)*inv(U)
 *  *
         CUT=0          CUT=0
         DO WHILE (CUT .LT. N)          DO WHILE (CUT .LT. N)
Line 432 Line 438
              END IF               END IF
            END DO             END DO
 *      *    
 *       L11T*invD1*L11->L11  *       L11**T*invD1*L11->L11
 *  *
         CALL ZTRMM('L',UPLO,'T','U',NNB, NNB,          CALL ZTRMM('L',UPLO,'T','U',NNB, NNB,
      $             ONE,A(CUT+1,CUT+1),LDA,WORK(U11+1,1),N+NB+1)       $             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          IF ( (CUT+NNB) .LT. N ) THEN
 *  *
 *          L21T*invD2*L21->A(CUT+I,CUT+J)  *          L21**T*invD2*L21->A(CUT+I,CUT+J)
 *  *
          CALL ZGEMM('T','N',NNB,NNB,N-NNB-CUT,ONE,A(CUT+NNB+1,CUT+1)           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  *        L11 =  L11**T*invD1*L11 + U01**T*invD*U01
 *  *
          DO I=1,NNB           DO I=1,NNB
             DO J=1,I              DO J=1,I
Line 453 Line 466
             END DO              END DO
          END DO           END DO
 *  *
 *        U01 =  L22T*invD2*L21  *        U01 =  L22**T*invD2*L21
 *  *
          CALL ZTRMM('L',UPLO,'T','U', N-NNB-CUT, NNB,           CALL ZTRMM('L',UPLO,'T','U', N-NNB-CUT, NNB,
      $             ONE,A(CUT+NNB+1,CUT+NNB+1),LDA,WORK,N+NB+1)       $             ONE,A(CUT+NNB+1,CUT+NNB+1),LDA,WORK,N+NB+1)
Line 466 Line 479
          END DO           END DO
        ELSE         ELSE
 *  *
 *        L11 =  L11T*invD1*L11  *        L11 =  L11**T*invD1*L11
 *  *
          DO I=1,NNB           DO I=1,NNB
             DO J=1,I              DO J=1,I
Line 480 Line 493
            CUT=CUT+NNB             CUT=CUT+NNB
        END DO         END DO
 *  *
 *        Apply PERMUTATIONS P and P': P * inv(U')*inv(D)*inv(U) *P'  *        Apply PERMUTATIONS P and P**T: P * inv(U**T)*inv(D)*inv(U) *P**T
 *   * 
             I=N              I=N
             DO WHILE ( I .GE. 1 )              DO WHILE ( I .GE. 1 )
                IF( IPIV(I) .GT. 0 ) THEN                 IF( IPIV(I) .GT. 0 ) THEN
                   IP=IPIV(I)                    IP=IPIV(I)
                  IF (I .LT. IP) CALL ZSYSWAPR( UPLO, N, A, I ,IP  )                   IF (I .LT. IP) CALL ZSYSWAPR( UPLO, N, A, LDA, I ,IP  )
                  IF (I .GT. IP) CALL ZSYSWAPR( UPLO, N, A, IP ,I )                   IF (I .GT. IP) CALL ZSYSWAPR( UPLO, N, A, LDA, IP ,I )
                ELSE                 ELSE
                  IP=-IPIV(I)                   IP=-IPIV(I)
                  IF ( I .LT. IP) CALL ZSYSWAPR( UPLO, N, A, I ,IP )                   IF ( I .LT. IP) CALL ZSYSWAPR( UPLO, N, A, LDA, I ,IP )
                  IF ( I .GT. IP) CALL ZSYSWAPR(  UPLO, N, A, IP ,I )                   IF ( I .GT. IP) CALL ZSYSWAPR( UPLO, N, A, LDA, IP ,I )
                  I=I-1                   I=I-1
                ENDIF                 ENDIF
                I=I-1                 I=I-1

Removed from v.1.2  
changed lines
  Added in v.1.3


CVSweb interface <joel.bertrand@systella.fr>