[BACK]Return to dlasyf_aa.f CVS log [TXT] [DIR] Up to [local] / rpl / lapack / lapack

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

version 1.2, 2017/06/17 11:06:27 version 1.3, 2018/05/29 06:55:20
Line 19 Line 19
 *  ===========  *  ===========
 *  *
 *       SUBROUTINE DLASYF_AA( UPLO, J1, M, NB, A, LDA, IPIV,  *       SUBROUTINE DLASYF_AA( UPLO, J1, M, NB, A, LDA, IPIV,
 *                             H, LDH, WORK, INFO )  *                             H, LDH, WORK )
 *  *
 *       .. Scalar Arguments ..  *       .. Scalar Arguments ..
 *       CHARACTER          UPLO  *       CHARACTER          UPLO
 *       INTEGER            J1, M, NB, LDA, LDH, INFO  *       INTEGER            J1, M, NB, LDA, LDH
 *       ..  *       ..
 *       .. Array Arguments ..  *       .. Array Arguments ..
 *       INTEGER            IPIV( * )  *       INTEGER            IPIV( * )
Line 99 Line 99
 *> \param[in] LDA  *> \param[in] LDA
 *> \verbatim  *> \verbatim
 *>          LDA is INTEGER  *>          LDA is INTEGER
 *>          The leading dimension of the array A.  LDA >= max(1,N).  *>          The leading dimension of the array A.  LDA >= max(1,M).
 *> \endverbatim  *> \endverbatim
 *>  *>
 *> \param[out] IPIV  *> \param[out] IPIV
 *> \verbatim  *> \verbatim
 *>          IPIV is INTEGER array, dimension (N)  *>          IPIV is INTEGER array, dimension (M)
 *>          Details of the row and column interchanges,  *>          Details of the row and column interchanges,
 *>          the row and column k were interchanged with the row and  *>          the row and column k were interchanged with the row and
 *>          column IPIV(k).  *>          column IPIV(k).
Line 127 Line 127
 *>          WORK is DOUBLE PRECISION workspace, dimension (M).  *>          WORK is DOUBLE PRECISION workspace, dimension (M).
 *> \endverbatim  *> \endverbatim
 *>  *>
 *> \param[out] INFO  
 *> \verbatim  
 *>          INFO is INTEGER  
 *>          = 0:  successful exit  
 *>          < 0:  if INFO = -i, the i-th argument had an illegal value  
 *>          > 0:  if INFO = i, D(i,i) is exactly zero.  The factorization  
 *>                has been completed, but the block diagonal matrix D is  
 *>                exactly singular, and division by zero will occur if it  
 *>                is used to solve a system of equations.  
 *> \endverbatim  
 *  *
 *  Authors:  *  Authors:
 *  ========  *  ========
Line 146 Line 136
 *> \author Univ. of Colorado Denver  *> \author Univ. of Colorado Denver
 *> \author NAG Ltd.  *> \author NAG Ltd.
 *  *
 *> \date December 2016  *> \date November 2017
 *  *
 *> \ingroup doubleSYcomputational  *> \ingroup doubleSYcomputational
 *  *
 *  =====================================================================  *  =====================================================================
       SUBROUTINE DLASYF_AA( UPLO, J1, M, NB, A, LDA, IPIV,        SUBROUTINE DLASYF_AA( UPLO, J1, M, NB, A, LDA, IPIV,
      $                      H, LDH, WORK, INFO )       $                      H, LDH, WORK )
 *  *
 *  -- LAPACK computational routine (version 3.7.0) --  *  -- LAPACK computational routine (version 3.8.0) --
 *  -- 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..--
 *     December 2016  *     November 2017
 *  *
       IMPLICIT NONE        IMPLICIT NONE
 *  *
 *     .. Scalar Arguments ..  *     .. Scalar Arguments ..
       CHARACTER          UPLO        CHARACTER          UPLO
       INTEGER            M, NB, J1, LDA, LDH, INFO        INTEGER            M, NB, J1, LDA, LDH
 *     ..  *     ..
 *     .. Array Arguments ..  *     .. Array Arguments ..
       INTEGER            IPIV( * )        INTEGER            IPIV( * )
Line 176 Line 166
       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )        PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
 *  *
 *     .. Local Scalars ..  *     .. Local Scalars ..
       INTEGER            J, K, K1, I1, I2        INTEGER            J, K, K1, I1, I2, MJ
       DOUBLE PRECISION   PIV, ALPHA        DOUBLE PRECISION   PIV, ALPHA
 *     ..  *     ..
 *     .. External Functions ..  *     .. External Functions ..
Line 185 Line 175
       EXTERNAL           LSAME, ILAENV, IDAMAX        EXTERNAL           LSAME, ILAENV, IDAMAX
 *     ..  *     ..
 *     .. External Subroutines ..  *     .. External Subroutines ..
       EXTERNAL           XERBLA        EXTERNAL           DGEMV, DAXPY, DCOPY, DSWAP, DSCAL, DLASET,
        $                   XERBLA
 *     ..  *     ..
 *     .. Intrinsic Functions ..  *     .. Intrinsic Functions ..
       INTRINSIC          MAX        INTRINSIC          MAX
 *     ..  *     ..
 *     .. Executable Statements ..  *     .. Executable Statements ..
 *  *
       INFO = 0  
       J = 1        J = 1
 *  *
 *     K1 is the first column of the panel to be factorized  *     K1 is the first column of the panel to be factorized
Line 216 Line 206
 *         > for the rest of the columns, J1 is 2, and J1+J-1 is J+1,  *         > for the rest of the columns, J1 is 2, and J1+J-1 is J+1,
 *  *
          K = J1+J-1           K = J1+J-1
            IF( J.EQ.M ) THEN
   *
   *            Only need to compute T(J, J)
   *
                MJ = 1
            ELSE
                MJ = M-J+1
            END IF
 *  *
 *        H(J:N, J) := A(J, J:N) - H(J:N, 1:(J-1)) * L(J1:(J-1), J),  *        H(J:M, J) := A(J, J:M) - H(J:M, 1:(J-1)) * L(J1:(J-1), J),
 *         where H(J:N, J) has been initialized to be A(J, J:N)  *         where H(J:M, J) has been initialized to be A(J, J:M)
 *  *
          IF( K.GT.2 ) THEN           IF( K.GT.2 ) THEN
 *  *
Line 228 Line 226
 *         > for the rest of the columns, K is J+1, skipping only the  *         > for the rest of the columns, K is J+1, skipping only the
 *           first column  *           first column
 *  *
             CALL DGEMV( 'No transpose', M-J+1, J-K1,              CALL DGEMV( 'No transpose', MJ, J-K1,
      $                 -ONE, H( J, K1 ), LDH,       $                 -ONE, H( J, K1 ), LDH,
      $                       A( 1, J ), 1,       $                       A( 1, J ), 1,
      $                  ONE, H( J, J ), 1 )       $                  ONE, H( J, J ), 1 )
          END IF           END IF
 *  *
 *        Copy H(i:n, i) into WORK  *        Copy H(i:M, i) into WORK
 *  *
          CALL DCOPY( M-J+1, H( J, J ), 1, WORK( 1 ), 1 )           CALL DCOPY( MJ, H( J, J ), 1, WORK( 1 ), 1 )
 *  *
          IF( J.GT.K1 ) THEN           IF( J.GT.K1 ) THEN
 *  *
 *           Compute WORK := WORK - L(J-1, J:N) * T(J-1,J),  *           Compute WORK := WORK - L(J-1, J:M) * T(J-1,J),
 *            where A(J-1, J) stores T(J-1, J) and A(J-2, J:N) stores U(J-1, J:N)  *            where A(J-1, J) stores T(J-1, J) and A(J-2, J:M) stores U(J-1, J:M)
 *  *
             ALPHA = -A( K-1, J )              ALPHA = -A( K-1, J )
             CALL DAXPY( M-J+1, ALPHA, A( K-2, J ), LDA, WORK( 1 ), 1 )              CALL DAXPY( MJ, ALPHA, A( K-2, J ), LDA, WORK( 1 ), 1 )
          END IF           END IF
 *  *
 *        Set A(J, J) = T(J, J)  *        Set A(J, J) = T(J, J)
Line 253 Line 251
 *  *
          IF( J.LT.M ) THEN           IF( J.LT.M ) THEN
 *  *
 *           Compute WORK(2:N) = T(J, J) L(J, (J+1):N)  *           Compute WORK(2:M) = T(J, J) L(J, (J+1):M)
 *            where A(J, J) stores T(J, J) and A(J-1, (J+1):N) stores U(J, (J+1):N)  *            where A(J, J) stores T(J, J) and A(J-1, (J+1):M) stores U(J, (J+1):M)
 *  *
             IF( K.GT.1 ) THEN              IF( K.GT.1 ) THEN
                ALPHA = -A( K, J )                 ALPHA = -A( K, J )
Line 262 Line 260
      $                                 WORK( 2 ), 1 )       $                                 WORK( 2 ), 1 )
             ENDIF              ENDIF
 *  *
 *           Find max(|WORK(2:n)|)  *           Find max(|WORK(2:M)|)
 *  *
             I2 = IDAMAX( M-J, WORK( 2 ), 1 ) + 1              I2 = IDAMAX( M-J, WORK( 2 ), 1 ) + 1
             PIV = WORK( I2 )              PIV = WORK( I2 )
Line 277 Line 275
                WORK( I2 ) = WORK( I1 )                 WORK( I2 ) = WORK( I1 )
                WORK( I1 ) = PIV                 WORK( I1 ) = PIV
 *  *
 *              Swap A(I1, I1+1:N) with A(I1+1:N, I2)  *              Swap A(I1, I1+1:M) with A(I1+1:M, I2)
 *  *
                I1 = I1+J-1                 I1 = I1+J-1
                I2 = I2+J-1                 I2 = I2+J-1
                CALL DSWAP( I2-I1-1, A( J1+I1-1, I1+1 ), LDA,                 CALL DSWAP( I2-I1-1, A( J1+I1-1, I1+1 ), LDA,
      $                              A( J1+I1, I2 ), 1 )       $                              A( J1+I1, I2 ), 1 )
 *  *
 *              Swap A(I1, I2+1:N) with A(I2, I2+1:N)  *              Swap A(I1, I2+1:M) with A(I2, I2+1:M)
 *  *
                CALL DSWAP( M-I2, A( J1+I1-1, I2+1 ), LDA,                 CALL DSWAP( M-I2, A( J1+I1-1, I2+1 ), LDA,
      $                           A( J1+I2-1, I2+1 ), LDA )       $                           A( J1+I2-1, I2+1 ), LDA )
Line 315 Line 313
 *           Set A(J, J+1) = T(J, J+1)  *           Set A(J, J+1) = T(J, J+1)
 *  *
             A( K, J+1 ) = WORK( 2 )              A( K, J+1 ) = WORK( 2 )
             IF( (A( K, J ).EQ.ZERO ) .AND.  
      $        ( (J.EQ.M) .OR. (A( K, J+1 ).EQ.ZERO))) THEN  
                 IF(INFO .EQ. 0) THEN  
                     INFO = J  
                 ENDIF  
             END IF  
 *  *
             IF( J.LT.NB ) THEN              IF( J.LT.NB ) THEN
 *  *
 *              Copy A(J+1:N, J+1) into H(J:N, J),  *              Copy A(J+1:M, J+1) into H(J:M, J),
 *  *
                CALL DCOPY( M-J, A( K+1, J+1 ), LDA,                 CALL DCOPY( M-J, A( K+1, J+1 ), LDA,
      $                          H( J+1, J+1 ), 1 )       $                          H( J+1, J+1 ), 1 )
             END IF              END IF
 *  *
 *           Compute L(J+2, J+1) = WORK( 3:N ) / T(J, J+1),  *           Compute L(J+2, J+1) = WORK( 3:M ) / T(J, J+1),
 *            where A(J, J+1) = T(J, J+1) and A(J+2:N, J) = L(J+2:N, J+1)  *            where A(J, J+1) = T(J, J+1) and A(J+2:M, J) = L(J+2:M, J+1)
 *  *
             IF( A( K, J+1 ).NE.ZERO ) THEN              IF( A( K, J+1 ).NE.ZERO ) THEN
                ALPHA = ONE / A( K, J+1 )                 ALPHA = ONE / A( K, J+1 )
Line 341 Line 333
                CALL DLASET( 'Full', 1, M-J-1, ZERO, ZERO,                 CALL DLASET( 'Full', 1, M-J-1, ZERO, ZERO,
      $                      A( K, J+2 ), LDA)       $                      A( K, J+2 ), LDA)
             END IF              END IF
          ELSE  
             IF( (A( K, J ).EQ.ZERO) .AND. (INFO.EQ.0) ) THEN  
                INFO = J  
             END IF  
          END IF           END IF
          J = J + 1           J = J + 1
          GO TO 10           GO TO 10
Line 366 Line 354
 *         > for the rest of the columns, J1 is 2, and J1+J-1 is J+1,  *         > for the rest of the columns, J1 is 2, and J1+J-1 is J+1,
 *  *
          K = J1+J-1           K = J1+J-1
            IF( J.EQ.M ) THEN
   *
   *            Only need to compute T(J, J)
   *
                MJ = 1
            ELSE
                MJ = M-J+1
            END IF
 *  *
 *        H(J:N, J) := A(J:N, J) - H(J:N, 1:(J-1)) * L(J, J1:(J-1))^T,  *        H(J:M, J) := A(J:M, J) - H(J:M, 1:(J-1)) * L(J, J1:(J-1))^T,
 *         where H(J:N, J) has been initialized to be A(J:N, J)  *         where H(J:M, J) has been initialized to be A(J:M, J)
 *  *
          IF( K.GT.2 ) THEN           IF( K.GT.2 ) THEN
 *  *
Line 378 Line 374
 *         > for the rest of the columns, K is J+1, skipping only the  *         > for the rest of the columns, K is J+1, skipping only the
 *           first column  *           first column
 *  *
             CALL DGEMV( 'No transpose', M-J+1, J-K1,              CALL DGEMV( 'No transpose', MJ, J-K1,
      $                 -ONE, H( J, K1 ), LDH,       $                 -ONE, H( J, K1 ), LDH,
      $                       A( J, 1 ), LDA,       $                       A( J, 1 ), LDA,
      $                  ONE, H( J, J ), 1 )       $                  ONE, H( J, J ), 1 )
          END IF           END IF
 *  *
 *        Copy H(J:N, J) into WORK  *        Copy H(J:M, J) into WORK
 *  *
          CALL DCOPY( M-J+1, H( J, J ), 1, WORK( 1 ), 1 )           CALL DCOPY( MJ, H( J, J ), 1, WORK( 1 ), 1 )
 *  *
          IF( J.GT.K1 ) THEN           IF( J.GT.K1 ) THEN
 *  *
 *           Compute WORK := WORK - L(J:N, J-1) * T(J-1,J),  *           Compute WORK := WORK - L(J:M, J-1) * T(J-1,J),
 *            where A(J-1, J) = T(J-1, J) and A(J, J-2) = L(J, J-1)  *            where A(J-1, J) = T(J-1, J) and A(J, J-2) = L(J, J-1)
 *  *
             ALPHA = -A( J, K-1 )              ALPHA = -A( J, K-1 )
             CALL DAXPY( M-J+1, ALPHA, A( J, K-2 ), 1, WORK( 1 ), 1 )              CALL DAXPY( MJ, ALPHA, A( J, K-2 ), 1, WORK( 1 ), 1 )
          END IF           END IF
 *  *
 *        Set A(J, J) = T(J, J)  *        Set A(J, J) = T(J, J)
Line 403 Line 399
 *  *
          IF( J.LT.M ) THEN           IF( J.LT.M ) THEN
 *  *
 *           Compute WORK(2:N) = T(J, J) L((J+1):N, J)  *           Compute WORK(2:M) = T(J, J) L((J+1):M, J)
 *            where A(J, J) = T(J, J) and A((J+1):N, J-1) = L((J+1):N, J)  *            where A(J, J) = T(J, J) and A((J+1):M, J-1) = L((J+1):M, J)
 *  *
             IF( K.GT.1 ) THEN              IF( K.GT.1 ) THEN
                ALPHA = -A( J, K )                 ALPHA = -A( J, K )
Line 412 Line 408
      $                                 WORK( 2 ), 1 )       $                                 WORK( 2 ), 1 )
             ENDIF              ENDIF
 *  *
 *           Find max(|WORK(2:n)|)  *           Find max(|WORK(2:M)|)
 *  *
             I2 = IDAMAX( M-J, WORK( 2 ), 1 ) + 1              I2 = IDAMAX( M-J, WORK( 2 ), 1 ) + 1
             PIV = WORK( I2 )              PIV = WORK( I2 )
Line 427 Line 423
                WORK( I2 ) = WORK( I1 )                 WORK( I2 ) = WORK( I1 )
                WORK( I1 ) = PIV                 WORK( I1 ) = PIV
 *  *
 *              Swap A(I1+1:N, I1) with A(I2, I1+1:N)  *              Swap A(I1+1:M, I1) with A(I2, I1+1:M)
 *  *
                I1 = I1+J-1                 I1 = I1+J-1
                I2 = I2+J-1                 I2 = I2+J-1
                CALL DSWAP( I2-I1-1, A( I1+1, J1+I1-1 ), 1,                 CALL DSWAP( I2-I1-1, A( I1+1, J1+I1-1 ), 1,
      $                              A( I2, J1+I1 ), LDA )       $                              A( I2, J1+I1 ), LDA )
 *  *
 *              Swap A(I2+1:N, I1) with A(I2+1:N, I2)  *              Swap A(I2+1:M, I1) with A(I2+1:M, I2)
 *  *
                CALL DSWAP( M-I2, A( I2+1, J1+I1-1 ), 1,                 CALL DSWAP( M-I2, A( I2+1, J1+I1-1 ), 1,
      $                           A( I2+1, J1+I2-1 ), 1 )       $                           A( I2+1, J1+I2-1 ), 1 )
Line 465 Line 461
 *           Set A(J+1, J) = T(J+1, J)  *           Set A(J+1, J) = T(J+1, J)
 *  *
             A( J+1, K ) = WORK( 2 )              A( J+1, K ) = WORK( 2 )
             IF( (A( J, K ).EQ.ZERO) .AND.  
      $        ( (J.EQ.M) .OR. (A( J+1, K ).EQ.ZERO)) ) THEN  
                 IF (INFO .EQ. 0)  
      $              INFO = J  
             END IF  
 *  *
             IF( J.LT.NB ) THEN              IF( J.LT.NB ) THEN
 *  *
 *              Copy A(J+1:N, J+1) into H(J+1:N, J),  *              Copy A(J+1:M, J+1) into H(J+1:M, J),
 *  *
                CALL DCOPY( M-J, A( J+1, K+1 ), 1,                 CALL DCOPY( M-J, A( J+1, K+1 ), 1,
      $                          H( J+1, J+1 ), 1 )       $                          H( J+1, J+1 ), 1 )
             END IF              END IF
 *  *
 *           Compute L(J+2, J+1) = WORK( 3:N ) / T(J, J+1),  *           Compute L(J+2, J+1) = WORK( 3:M ) / T(J, J+1),
 *            where A(J, J+1) = T(J, J+1) and A(J+2:N, J) = L(J+2:N, J+1)  *            where A(J, J+1) = T(J, J+1) and A(J+2:M, J) = L(J+2:M, J+1)
 *  *
             IF( A( J+1, K ).NE.ZERO ) THEN              IF( A( J+1, K ).NE.ZERO ) THEN
                ALPHA = ONE / A( J+1, K )                 ALPHA = ONE / A( J+1, K )
Line 490 Line 481
                CALL DLASET( 'Full', M-J-1, 1, ZERO, ZERO,                 CALL DLASET( 'Full', M-J-1, 1, ZERO, ZERO,
      $                      A( J+2, K ), LDA )       $                      A( J+2, K ), LDA )
             END IF              END IF
          ELSE  
             IF( (A( J, K ).EQ.ZERO) .AND. (INFO.EQ.0) ) THEN  
                INFO = J  
             END IF  
          END IF           END IF
          J = J + 1           J = J + 1
          GO TO 30           GO TO 30
Removed from v.1.2  
changed lines
  Added in v.1.3

CVSweb interface <joel.bertrand@systella.fr>