diff options
author | igor175 <igor175@8a072113-8704-0410-8d35-dd094bca7971> | 2011-12-07 12:52:24 +0000 |
---|---|---|
committer | igor175 <igor175@8a072113-8704-0410-8d35-dd094bca7971> | 2011-12-07 12:52:24 +0000 |
commit | 3a4622a6c1cb8ba661b880e2598223e6a4494565 (patch) | |
tree | 8397c70e7339b9509772e7f142cdf2bfe1f12687 | |
parent | 1e32dac8817922ec7ae53406b4a1bd711d1c0c72 (diff) | |
download | lapack-3a4622a6c1cb8ba661b880e2598223e6a4494565.tar.gz lapack-3a4622a6c1cb8ba661b880e2598223e6a4494565.tar.bz2 lapack-3a4622a6c1cb8ba661b880e2598223e6a4494565.zip |
added standard precision (COMPLEX) routines for rook pivoting \\nalgorithm for symmetric indefinite matrices: slasyf_rook.f ssytri_rook.f ssytrs_rook.f ssytrf_rook.f ssysv_rook.f ssycon_rook.f ssytf2_rook.f
-rw-r--r-- | SRC/clasyf_rook.f | 894 | ||||
-rw-r--r-- | SRC/csycon_rook.f | 260 | ||||
-rw-r--r-- | SRC/csysv_rook.f | 293 | ||||
-rw-r--r-- | SRC/csytf2_rook.f | 824 | ||||
-rw-r--r-- | SRC/csytrf_rook.f | 393 | ||||
-rw-r--r-- | SRC/csytri_rook.f | 451 | ||||
-rw-r--r-- | SRC/csytrs_rook.f | 484 |
7 files changed, 3599 insertions, 0 deletions
diff --git a/SRC/clasyf_rook.f b/SRC/clasyf_rook.f new file mode 100644 index 00000000..b46fb987 --- /dev/null +++ b/SRC/clasyf_rook.f @@ -0,0 +1,894 @@ +*> \brief \b CLASYF_ROOK +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download CLASYF_ROOK + dependencies +*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/clasyf_rook.f"> +*> [TGZ]</a> +*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/clasyf_rook.f"> +*> [ZIP]</a> +*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/clasyf_rook.f"> +*> [TXT]</a> +*> \endhtmlonly +* +* Definition: +* =========== +* +* SUBROUTINE CLASYF_ROOK( UPLO, N, NB, KB, A, LDA, IPIV, W, LDW, INFO ) +* +* .. Scalar Arguments .. +* CHARACTER UPLO +* INTEGER INFO, KB, LDA, LDW, N, NB +* .. +* .. Array Arguments .. +* INTEGER IPIV( * ) +* COMPLEX A( LDA, * ), W( LDW, * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> CLASYF_ROOK computes a partial factorization of a complex symmetric +*> matrix A using the bounded Bunch-Kaufman ("rook") diagonal +*> pivoting method. The partial factorization has the form: +*> +*> A = ( I U12 ) ( A11 0 ) ( I 0 ) if UPLO = 'U', or: +*> ( 0 U22 ) ( 0 D ) ( U12**T U22**T ) +*> +*> A = ( L11 0 ) ( D 0 ) ( L11**T L21**T ) if UPLO = 'L' +*> ( L21 I ) ( 0 A22 ) ( 0 I ) +*> +*> where the order of D is at most NB. The actual order is returned in +*> the argument KB, and is either NB or NB-1, or N if N <= NB. +*> +*> CLASYF_ROOK is an auxiliary routine called by CSYTRF_ROOK. It uses +*> blocked code (calling Level 3 BLAS) to update the submatrix +*> A11 (if UPLO = 'U') or A22 (if UPLO = 'L'). +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] UPLO +*> \verbatim +*> UPLO is CHARACTER*1 +*> Specifies whether the upper or lower triangular part of the +*> symmetric matrix A is stored: +*> = 'U': Upper triangular +*> = 'L': Lower triangular +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The order of the matrix A. N >= 0. +*> \endverbatim +*> +*> \param[in] NB +*> \verbatim +*> NB is INTEGER +*> The maximum number of columns of the matrix A that should be +*> factored. NB should be at least 2 to allow for 2-by-2 pivot +*> blocks. +*> \endverbatim +*> +*> \param[out] KB +*> \verbatim +*> KB is INTEGER +*> The number of columns of A that were actually factored. +*> KB is either NB-1 or NB, or N if N <= NB. +*> \endverbatim +*> +*> \param[in,out] A +*> \verbatim +*> A is COMPLEX array, dimension (LDA,N) +*> On entry, the symmetric matrix A. If UPLO = 'U', the leading +*> n-by-n upper triangular part of A contains the upper +*> triangular part of the matrix A, and the strictly lower +*> triangular part of A is not referenced. If UPLO = 'L', the +*> leading n-by-n lower triangular part of A contains the lower +*> triangular part of the matrix A, and the strictly upper +*> triangular part of A is not referenced. +*> On exit, A contains details of the partial factorization. +*> \endverbatim +*> +*> \param[in] LDA +*> \verbatim +*> LDA is INTEGER +*> The leading dimension of the array A. LDA >= max(1,N). +*> \endverbatim +*> +*> \param[out] IPIV +*> \verbatim +*> IPIV is INTEGER array, dimension (N) +*> Details of the interchanges and the block structure of D. +*> +*> If UPLO = 'U': +*> Only the last KB elements of IPIV are set +*> +*> If IPIV(k) > 0, then rows and columns k and IPIV(k) +*> were interchanged and D(k,k) is a 1-by-1 diagonal block. +*> +*> If IPIV(k) < 0 and IPIV(k-1) < 0, then rows and +*> columns k and -IPIV(k) were interchanged and rows and +*> columns k-1 and -IPIV(k-1) were inerchaged, +*> D(k-1:k,k-1:k) is a 2-by-2 diagonal block. +*> +*> If UPLO = 'L': +*> Only the first KB elements are set. +*> +*> If IPIV(k) > 0, then rows and columns k and IPIV(k) +*> were interchanged and D(k,k) is a 1-by-1 diagonal block. +*> +*> If IPIV(k) < 0 and IPIV(k+1) < 0, then rows and +*> columns k and -IPIV(k) were interchanged and rows and +*> columns k+1 and -IPIV(k+1) were inerchaged, +*> D(k:k+1,k:k+1) is a 2-by-2 diagonal block. +*> \endverbatim +*> +*> \param[out] W +*> \verbatim +*> W is COMPLEX array, dimension (LDW,NB) +*> \endverbatim +*> +*> \param[in] LDW +*> \verbatim +*> LDW is INTEGER +*> The leading dimension of the array W. LDW >= max(1,N). +*> \endverbatim +*> +*> \param[out] INFO +*> \verbatim +*> INFO is INTEGER +*> = 0: successful exit +*> > 0: if INFO = k, D(k,k) is exactly zero. The factorization +*> has been completed, but the block diagonal matrix D is +*> exactly singular. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2011 +* +*> \ingroup complexSYcomputational +* +*> \par Contributors: +* ================== +*> +*> \verbatim +*> +*> November 2011, Igor Kozachenko, +*> Computer Science Division, +*> University of California, Berkeley +*> +*> September 2007, Sven Hammarling, Nicholas J. Higham, Craig Lucas, +*> School of Mathematics, +*> University of Manchester +*> +*> \endverbatim +* +* ===================================================================== + SUBROUTINE CLASYF_ROOK( UPLO, N, NB, KB, A, LDA, IPIV, W, LDW, + $ INFO ) +* +* -- LAPACK computational routine (version 3.3.1) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2011 +* +* .. Scalar Arguments .. + CHARACTER UPLO + INTEGER INFO, KB, LDA, LDW, N, NB +* .. +* .. Array Arguments .. + INTEGER IPIV( * ) + COMPLEX A( LDA, * ), W( LDW, * ) +* .. +* +* ===================================================================== +* +* .. Parameters .. + REAL ZERO, ONE + PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) + REAL EIGHT, SEVTEN + PARAMETER ( EIGHT = 8.0E+0, SEVTEN = 17.0E+0 ) + COMPLEX CONE, CZERO + PARAMETER ( CONE = ( 1.0E+0, 0.0E+0 ), + $ CZERO = ( 0.0E+0, 0.0E+0 ) ) +* .. +* .. Local Scalars .. + LOGICAL DONE + INTEGER IMAX, ITEMP, J, JB, JJ, JMAX, JP1, JP2, K, KK, + $ KW, KKW, KP, KSTEP, P, II + REAL ABSAKK, ALPHA, COLMAX, ROWMAX, STEMP, SFMIN + COMPLEX D11, D12, D21, D22, R1 + $ T, Z +* .. +* .. External Functions .. + LOGICAL LSAME + INTEGER ICAMAX + REAL SLAMCH + EXTERNAL LSAME, ICAMAX, SLAMCH +* .. +* .. External Subroutines .. + EXTERNAL CCOPY, CGEMM, CGEMV, CSCAL, CSWAP +* .. +* .. Intrinsic Functions .. + INTRINSIC ABS, MAX, MIN, SQRT, AIMAG, REAL +* .. +* .. Statement Functions .. + REAL CABS1 +* .. +* .. Statement Function definitions .. + CABS1( Z ) = ABS( REAL( Z ) ) + ABS( AIMAG( Z ) ) +* .. +* .. Executable Statements .. +* + INFO = 0 +* +* Initialize ALPHA for use in choosing pivot block size. +* + ALPHA = ( ONE+SQRT( SEVTEN ) ) / EIGHT +* +* Compute machine safe minimum +* + SFMIN = SLAMCH( 'S' ) +* + IF( LSAME( UPLO, 'U' ) ) THEN +* +* Factorize the trailing columns of A using the upper triangle +* of A and working backwards, and compute the matrix W = U12*D +* for use in updating A11 +* +* K is the main loop index, decreasing from N in steps of 1 or 2 +* + + K = N + 10 CONTINUE +* +* KW is the column of W which corresponds to column K of A +* + KW = NB + K - N +* +* Exit from loop +* + IF( ( K.LE.N-NB+1 .AND. NB.LT.N ) .OR. K.LT.1 ) + $ GO TO 30 + KSTEP = 1 + P = K +* +* Copy column K of A to column KW of W and update it +* + CALL CCOPY( K, A( 1, K ), 1, W( 1, KW ), 1 ) + IF( K.LT.N ) + $ CALL CGEMV( 'No transpose', K, N-K, -CONE, A( 1, K+1 ), LDA, + $ W( K, KW+1 ), LDW, CONE, W( 1, KW ), 1 ) +* +* Determine rows and columns to be interchanged and whether +* a 1-by-1 or 2-by-2 pivot block will be used +* + ABSAKK = CABS1( W( K, KW ) ) +* +* IMAX is the row-index of the largest off-diagonal element in +* column K, and COLMAX is its absolute value. +* Determine both COLMAX and IMAX. +* + IF( K.GT.1 ) THEN + IMAX = ICAMAX( K-1, W( 1, KW ), 1 ) + COLMAX = CABS1( W( IMAX, KW ) ) + ELSE + COLMAX = ZERO + END IF +* + IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN +* +* Column K is zero or underflow: set INFO and continue +* + IF( INFO.EQ.0 ) + $ INFO = K + KP = K + CALL CCOPY( K, W( 1, KW ), 1, A( 1, K ), 1 ) + ELSE +* ===================================================================== +* +* Test for interchange +* +* Equivalent to testing for (used to handle NaN and Inf) +* ABSAKK.GE.ALPHA*COLMAX +* + IF( .NOT.( ABSAKK.LT.ALPHA*COLMAX ) ) THEN +* +* no interchange, use 1-by-1 pivot block +* + KP = K +* + ELSE +* + DONE = .FALSE. +* +* Loop until pivot found +* + 12 CONTINUE +* +* Begin pivot search loop body +* +* +* Copy column IMAX to column KW-1 of W and update it +* + CALL CCOPY( IMAX, A( 1, IMAX ), 1, W( 1, KW-1 ), 1 ) + CALL CCOPY( K-IMAX, A( IMAX, IMAX+1 ), LDA, + $ W( IMAX+1, KW-1 ), 1 ) +* + IF( K.LT.N ) + $ CALL CGEMV( 'No transpose', K, N-K, -CONE, + $ A( 1, K+1 ), LDA, W( IMAX, KW+1 ), LDW, + $ CONE,W( 1, KW-1 ), 1 ) +* +* JMAX is the column-index of the largest off-diagonal +* element in row IMAX, and ROWMAX is its absolute value. +* Determine both ROWMAX and JMAX. +* + IF( IMAX.NE.K ) THEN + JMAX = IMAX + ICAMAX( K-IMAX, W( IMAX+1, KW-1 ), + $ 1 ) + ROWMAX = CABS1( W( JMAX, KW-1 ) ) + ELSE + ROWMAX = ZERO + END IF +* + IF( IMAX.GT.1 ) THEN + ITEMP = ICAMAX( IMAX-1, W( 1, KW-1 ), 1 ) + STEMP = CABS1( W( ITEMP, KW-1 ) ) + IF( STEMP.GT.ROWMAX ) THEN + ROWMAX = STEMP + JMAX = ITEMP + END IF + END IF +* +* Equivalent to testing for (used to handle NaN and Inf) +* CABS1( W( IMAX, KW-1 ) ).GE.ALPHA*ROWMAX +* + IF( .NOT.(CABS1( W( IMAX, KW-1 ) ).LT.ALPHA*ROWMAX ) ) + $ THEN +* +* interchange rows and columns K and IMAX, +* use 1-by-1 pivot block +* + KP = IMAX +* +* copy column KW-1 of W to column KW +* + CALL CCOPY( K, W( 1, KW-1 ), 1, W( 1, KW ), 1 ) +* + DONE = .TRUE. +* +* Equivalent to testing for ROWMAX .EQ. COLMAX, +* used to handle NaN and Inf +* + ELSE IF( ( P.EQ.JMAX ) .OR. ( ROWMAX.LE.COLMAX ) ) + $ THEN +* +* interchange rows and columns K-1 and IMAX, +* use 2-by-2 pivot block +* + KP = IMAX + KSTEP = 2 + DONE = .TRUE. + ELSE +* +* Pivot not found: set params and repeat +* + P = IMAX + COLMAX = ROWMAX + IMAX = JMAX +* +* Copy updated JMAXth (next IMAXth) column to Kth of W +* + CALL CCOPY( K, W( 1, KW-1 ), 1, W( 1, KW ), 1 ) +* + END IF +* +* End pivot search loop body +* + IF( .NOT. DONE ) GOTO 12 +* ===================================================================== + END IF +* + KK = K - KSTEP + 1 +* +* KKW is the column of W which corresponds to column KK of A +* + KKW = NB + KK - N +* + IF( ( KSTEP.EQ.2 ) .AND. ( P.NE.K ) ) THEN +* +* Copy non-updated column K to column P +* + CALL CCOPY( K-P, A( P+1, K ), 1, A( P, P+1 ), LDA ) + CALL CCOPY( P, A( 1, K ), 1, A( 1, P ), 1 ) +* +* Interchange rows K and P in last N-K+1 columns of A +* and last N-K+2 columns of W +* + CALL CSWAP( N-K+1, A( K, K ), LDA, A( P, K ), LDA ) + CALL CSWAP( N-KK+1, W( K, KKW ), LDW, W( P, KKW ), LDW ) + END IF +* +* Updated column KP is already stored in column KKW of W +* + IF( KP.NE.KK ) THEN +* +* Copy non-updated column KK to column KP +* + A( KP, K ) = A( KK, K ) + CALL CCOPY( K-1-KP, A( KP+1, KK ), 1, A( KP, KP+1 ), + $ LDA ) + CALL CCOPY( KP, A( 1, KK ), 1, A( 1, KP ), 1 ) +* +* Interchange rows KK and KP in last N-KK+1 columns +* of A and W +* + CALL CSWAP( N-KK+1, A( KK, KK ), LDA, A( KP, KK ), LDA ) + CALL CSWAP( N-KK+1, W( KK, KKW ), LDW, W( KP, KKW ), + $ LDW ) + END IF +* + IF( KSTEP.EQ.1 ) THEN +* +* 1-by-1 pivot block D(k): column KW of W now holds +* +* W(k) = U(k)*D(k) +* +* where U(k) is the k-th column of U +* +* Store U(k) in column k of A +* + CALL CCOPY( K, W( 1, KW ), 1, A( 1, K ), 1 ) + IF( K.GT.1 ) THEN + IF( CABS1( A( K, K ) ).GE.SFMIN ) THEN + R1 = CONE / A( K, K ) + CALL CSCAL( K-1, R1, A( 1, K ), 1 ) + ELSE IF( A( K, K ).NE.CZERO ) THEN + DO 14 II = 1, K - 1 + A( II, K ) = A( II, K ) / A( K, K ) + 14 CONTINUE + END IF + END IF +* + ELSE +* +* 2-by-2 pivot block D(k): columns KW and KW-1 of W now +* hold +* +* ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k) +* +* where U(k) and U(k-1) are the k-th and (k-1)-th columns +* of U +* + IF( K.GT.2 ) THEN +* +* Store U(k) and U(k-1) in columns k and k-1 of A +* + D12 = W( K-1, KW ) + D11 = W( K, KW ) / D12 + D22 = W( K-1, KW-1 ) / D12 + T = CONE / ( D11*D22-CONE ) + DO 20 J = 1, K - 2 + A( J, K-1 ) = T*( (D11*W( J, KW-1 )-W( J, KW ) ) / + $ D12 ) + A( J, K ) = T*( ( D22*W( J, KW )-W( J, KW-1 ) ) / + $ D12 ) + 20 CONTINUE + END IF +* +* Copy D(k) to A +* + A( K-1, K-1 ) = W( K-1, KW-1 ) + A( K-1, K ) = W( K-1, KW ) + A( K, K ) = W( K, KW ) + END IF + END IF +* +* Store details of the interchanges in IPIV +* + IF( KSTEP.EQ.1 ) THEN + IPIV( K ) = KP + ELSE + IPIV( K ) = -P + IPIV( K-1 ) = -KP + END IF +* +* Decrease K and return to the start of the main loop +* + K = K - KSTEP + GO TO 10 +* + 30 CONTINUE +* +* Update the upper triangle of A11 (= A(1:k,1:k)) as +* +* A11 := A11 - U12*D*U12**T = A11 - U12*W**T +* +* computing blocks of NB columns at a time +* + DO 50 J = ( ( K-1 ) / NB )*NB + 1, 1, -NB + JB = MIN( NB, K-J+1 ) +* +* Update the upper triangle of the diagonal block +* + DO 40 JJ = J, J + JB - 1 + CALL CGEMV( 'No transpose', JJ-J+1, N-K, -CONE, + $ A( J, K+1 ), LDA, W( JJ, KW+1 ), LDW, CONE, + $ A( J, JJ ), 1 ) + 40 CONTINUE +* +* Update the rectangular superdiagonal block +* + IF( J.GE.2 ) + $ CALL CGEMM( 'No transpose', 'Transpose', J-1, JB, + $ N-K, -CONE, A( 1, K+1 ), LDA, W( J, KW+1 ), LDW, + $ CONE, A( 1, J ), LDA ) + 50 CONTINUE +* +* Put U12 in standard form by partially undoing the interchanges +* in columns k+1:n +* + J = K + 1 + 60 CONTINUE +* + KSTEP = 1 + JP1 = 1 + JJ = J + JP2 = IPIV( J ) + IF( JP2.LT.0 ) THEN + JP2 = -JP2 + J = J + 1 + JP1 = -IPIV( J ) + KSTEP = 2 + END IF +* + J = J + 1 + IF( JP2.NE.JJ .AND. J.LE.N ) + $ CALL CSWAP( N-J+1, A( JP2, J ), LDA, A( JJ, J ), LDA ) + JJ = J - 1 + IF( JP1.NE.JJ .AND. KSTEP.EQ.2 ) + $ CALL CSWAP( N-J+1, A( JP1, J ), LDA, A( JJ, J ), LDA ) + IF( J.LE.N ) + $ GO TO 60 +* +* Set KB to the number of columns factorized +* + KB = N - K +* + ELSE +* +* Factorize the leading columns of A using the lower triangle +* of A and working forwards, and compute the matrix W = L21*D +* for use in updating A22 +* +* K is the main loop index, increasing from 1 in steps of 1 or 2 +* + K = 1 + 70 CONTINUE +* +* Exit from loop +* + IF( ( K.GE.NB .AND. NB.LT.N ) .OR. K.GT.N ) + $ GO TO 90 + KSTEP = 1 + P = K +* +* Copy column K of A to column K of W and update it +* + CALL CCOPY( N-K+1, A( K, K ), 1, W( K, K ), 1 ) + IF( K.GT.1 ) + $ CALL CGEMV( 'No transpose', N-K+1, K-1, -CONE, + $ A( K, 1 ), LDA, W( K, 1 ), LDW, CONE, W( K, K ), 1 ) + +* +* Determine rows and columns to be interchanged and whether +* a 1-by-1 or 2-by-2 pivot block will be used +* + ABSAKK = CABS1( W( K, K ) ) +* +* IMAX is the row-index of the largest off-diagonal element in +* column K, and COLMAX is its absolute value. +* Determine both COLMAX and IMAX. +* + IF( K.LT.N ) THEN + IMAX = K + ICAMAX( N-K, W( K+1, K ), 1 ) + COLMAX = CABS1( W( IMAX, K ) ) + ELSE + COLMAX = ZERO + END IF +* + IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN +* +* Column K is zero or underflow: set INFO and continue +* + IF( INFO.EQ.0 ) + $ INFO = K + KP = K + CALL CCOPY( N-K+1, W( K, K ), 1, A( K, K ), 1 ) +* + ELSE +* ===================================================================== +* +* Test for interchange +* +* Equivalent to testing for (used to handle NaN and Inf) +* ABSAKK.GE.ALPHA*COLMAX +* + IF( .NOT.( ABSAKK.LT.ALPHA*COLMAX ) ) THEN +* +* no interchange, use 1-by-1 pivot block +* + KP = K +* + ELSE +* + DONE = .FALSE. +* +* Loop until pivot found +* + 72 CONTINUE +* +* Begin pivot search loop body +* +* +* Copy column IMAX to column K+1 of W and update it +* + CALL CCOPY( IMAX-K, A( IMAX, K ), LDA, W( K, K+1 ), 1) + CALL CCOPY( N-IMAX+1, A( IMAX, IMAX ), 1, + $ W( IMAX, K+1 ), 1 ) + IF( K.GT.1 ) + $ CALL CGEMV( 'No transpose', N-K+1, K-1, -CONE, + $ A( K, 1 ), LDA, W( IMAX, 1 ), LDW, + $ CONE, W( K, K+1 ), 1 ) +* +* JMAX is the column-index of the largest off-diagonal +* element in row IMAX, and ROWMAX is its absolute value. +* Determine both ROWMAX and JMAX. +* + IF( IMAX.NE.K ) THEN + JMAX = K - 1 + ICAMAX( IMAX-K, W( K, K+1 ), 1 ) + ROWMAX = CABS1( W( JMAX, K+1 ) ) + ELSE + ROWMAX = ZERO + END IF +* + IF( IMAX.LT.N ) THEN + ITEMP = IMAX + ICAMAX( N-IMAX, W( IMAX+1, K+1 ), 1) + STEMP = CABS1( W( ITEMP, K+1 ) ) + IF( STEMP.GT.ROWMAX ) THEN + ROWMAX = STEMP + JMAX = ITEMP + END IF + END IF +* +* Equivalent to testing for (used to handle NaN and Inf) +* CABS1( W( IMAX, K+1 ) ).GE.ALPHA*ROWMAX +* + IF( .NOT.( CABS1( W( IMAX, K+1 ) ).LT.ALPHA*ROWMAX ) ) + $ THEN +* +* interchange rows and columns K and IMAX, +* use 1-by-1 pivot block +* + KP = IMAX +* +* copy column K+1 of W to column K +* + CALL CCOPY( N-K+1, W( K, K+1 ), 1, W( K, K ), 1 ) +* + DONE = .TRUE. +* +* Equivalent to testing for ROWMAX .EQ. COLMAX, +* used to handle NaN and Inf +* + ELSE IF( ( P.EQ.JMAX ) .OR. ( ROWMAX.LE.COLMAX ) ) + $ THEN +* +* interchange rows and columns K+1 and IMAX, +* use 2-by-2 pivot block +* + KP = IMAX + KSTEP = 2 + DONE = .TRUE. + ELSE +* +* Pivot not found: set params and repeat +* + P = IMAX + COLMAX = ROWMAX + IMAX = JMAX +* +* Copy updated JMAXth (next IMAXth) column to Kth of W +* + CALL CCOPY( N-K+1, W( K, K+1 ), 1, W( K, K ), 1 ) +* + END IF +* +* End pivot search loop body +* + IF( .NOT. DONE ) GOTO 72 +* ===================================================================== + END IF +* + KK = K + KSTEP - 1 +* + IF( ( KSTEP.EQ.2 ) .AND. ( P.NE.K ) ) THEN +* +* Copy non-updated column K to column P +* + CALL CCOPY( P-K, A( K, K ), 1, A( P, K ), LDA ) + CALL CCOPY( N-P+1, A( P, K ), 1, A( P, P ), 1 ) +* +* Interchange rows K and P in first K columns of A +* and first K+1 columns of W +* + CALL CSWAP( K, A( K, 1 ), LDA, A( P, 1 ), LDA ) + CALL CSWAP( KK, W( K, 1 ), LDW, W( P, 1 ), LDW ) + END IF +* +* Updated column KP is already stored in column KK of W +* + IF( KP.NE.KK ) THEN +* +* Copy non-updated column KK to column KP +* + A( KP, K ) = A( KK, K ) + CALL CCOPY( KP-K-1, A( K+1, KK ), 1, A( KP, K+1 ), LDA ) + CALL CCOPY( N-KP+1, A( KP, KK ), 1, A( KP, KP ), 1 ) +* +* Interchange rows KK and KP in first KK columns of A and W +* + CALL CSWAP( KK, A( KK, 1 ), LDA, A( KP, 1 ), LDA ) + CALL CSWAP( KK, W( KK, 1 ), LDW, W( KP, 1 ), LDW ) + END IF +* + IF( KSTEP.EQ.1 ) THEN +* +* 1-by-1 pivot block D(k): column k of W now holds +* +* W(k) = L(k)*D(k) +* +* where L(k) is the k-th column of L +* +* Store L(k) in column k of A +* + CALL CCOPY( N-K+1, W( K, K ), 1, A( K, K ), 1 ) + IF( K.LT.N ) THEN + IF( CABS1( A( K, K ) ).GE.SFMIN ) THEN + R1 = CONE / A( K, K ) + CALL CSCAL( N-K, R1, A( K+1, K ), 1 ) + ELSE IF( A( K, K ).NE.CZERO ) THEN + DO 74 II = K + 1, N + A( II, K ) = A( II, K ) / A( K, K ) + 74 CONTINUE + END IF + END IF +* + ELSE +* +* 2-by-2 pivot block D(k): columns k and k+1 of W now hold +* +* ( W(k) W(k+1) ) = ( L(k) L(k+1) )*D(k) +* +* where L(k) and L(k+1) are the k-th and (k+1)-th columns +* of L +* + IF( K.LT.N-1 ) THEN +* +* Store L(k) and L(k+1) in columns k and k+1 of A +* + D21 = W( K+1, K ) + D11 = W( K+1, K+1 ) / D21 + D22 = W( K, K ) / D21 + T = CONE / ( D11*D22-CONE ) + DO 80 J = K + 2, N + A( J, K ) = T*( ( D11*W( J, K )-W( J, K+1 ) ) / + $ D21 ) + A( J, K+1 ) = T*( ( D22*W( J, K+1 )-W( J, K ) ) / + $ D21 ) + 80 CONTINUE + END IF +* +* Copy D(k) to A +* + A( K, K ) = W( K, K ) + A( K+1, K ) = W( K+1, K ) + A( K+1, K+1 ) = W( K+1, K+1 ) + END IF + END IF +* +* Store details of the interchanges in IPIV +* + IF( KSTEP.EQ.1 ) THEN + IPIV( K ) = KP + ELSE + IPIV( K ) = -P + IPIV( K+1 ) = -KP + END IF +* +* Increase K and return to the start of the main loop +* + K = K + KSTEP + GO TO 70 +* + 90 CONTINUE +* +* Update the lower triangle of A22 (= A(k:n,k:n)) as +* +* A22 := A22 - L21*D*L21**T = A22 - L21*W**T +* +* computing blocks of NB columns at a time +* + DO 110 J = K, N, NB + JB = MIN( NB, N-J+1 ) +* +* Update the lower triangle of the diagonal block +* + DO 100 JJ = J, J + JB - 1 + CALL CGEMV( 'No transpose', J+JB-JJ, K-1, -CONE, + $ A( JJ, 1 ), LDA, W( JJ, 1 ), LDW, CONE, + $ A( JJ, JJ ), 1 ) + 100 CONTINUE +* +* Update the rectangular subdiagonal block +* + IF( J+JB.LE.N ) + $ CALL CGEMM( 'No transpose', 'Transpose', N-J-JB+1, JB, + $ K-1, -CONE, A( J+JB, 1 ), LDA, W( J, 1 ), LDW, + $ CONE, A( J+JB, J ), LDA ) + 110 CONTINUE +* +* Put L21 in standard form by partially undoing the interchanges +* in columns 1:k-1 +* + J = K - 1 + 120 CONTINUE +* + KSTEP = 1 + JP1 = 1 + JJ = J + JP2 = IPIV( J ) + IF( JP2.LT.0 ) THEN + JP2 = -JP2 + J = J - 1 + JP1 = -IPIV( J ) + KSTEP = 2 + END IF +* + J = J - 1 + IF( JP2.NE.JJ .AND. J.GE.1 ) + $ CALL CSWAP( J, A( JP2, 1 ), LDA, A( JJ, 1 ), LDA ) + JJ = J + 1 + IF( JP1.NE.JJ .AND. KSTEP.EQ.2 ) + $ CALL CSWAP( J, A( JP1, 1 ), LDA, A( JJ, 1 ), LDA ) + IF( J.GE.1 ) + $ GO TO 120 +* +* Set KB to the number of columns factorized +* + KB = K - 1 +* + END IF + RETURN +* +* End of CLASYF_ROOK +* + END diff --git a/SRC/csycon_rook.f b/SRC/csycon_rook.f new file mode 100644 index 00000000..50690dfd --- /dev/null +++ b/SRC/csycon_rook.f @@ -0,0 +1,260 @@ +*> \brief \b CSYCON_ROOK +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download CSYCON_ROOK + dependencies +*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/csycon_rook.f"> +*> [TGZ]</a> +*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/csycon_rook.f"> +*> [ZIP]</a> +*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/csycon_rook.f"> +*> [TXT]</a> +*> \endhtmlonly +* +* Definition: +* =========== +* +* SUBROUTINE CSYCON_ROOK( UPLO, N, A, LDA, IPIV, ANORM, RCOND, +* WORK, IWORK, INFO ) +* +* .. Scalar Arguments .. +* CHARACTER UPLO +* INTEGER INFO, LDA, N +* REAL ANORM, RCOND +* .. +* .. Array Arguments .. +* INTEGER IPIV( * ), IWORK( * ) +* COMPLEX A( LDA, * ), WORK( * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> CSYCON_ROOK estimates the reciprocal of the condition number (in the +*> 1-norm) of a complex symmetric matrix A using the factorization +*> A = U*D*U**T or A = L*D*L**T computed by CSYTRF_ROOK. +*> +*> An estimate is obtained for norm(inv(A)), and the reciprocal of the +*> condition number is computed as RCOND = 1 / (ANORM * norm(inv(A))). +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] UPLO +*> \verbatim +*> UPLO is CHARACTER*1 +*> Specifies whether the details of the factorization are stored +*> as an upper or lower triangular matrix. +*> = 'U': Upper triangular, form is A = U*D*U**T; +*> = 'L': Lower triangular, form is A = L*D*L**T. +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The order of the matrix A. N >= 0. +*> \endverbatim +*> +*> \param[in] A +*> \verbatim +*> A is COMPLEX array, dimension (LDA,N) +*> The block diagonal matrix D and the multipliers used to +*> obtain the factor U or L as computed by CSYTRF_ROOK. +*> \endverbatim +*> +*> \param[in] LDA +*> \verbatim +*> LDA is INTEGER +*> The leading dimension of the array A. LDA >= max(1,N). +*> \endverbatim +*> +*> \param[in] IPIV +*> \verbatim +*> IPIV is INTEGER array, dimension (N) +*> Details of the interchanges and the block structure of D +*> as determined by CSYTRF_ROOK. +*> \endverbatim +*> +*> \param[in] ANORM +*> \verbatim +*> ANORM is REAL +*> The 1-norm of the original matrix A. +*> \endverbatim +*> +*> \param[out] RCOND +*> \verbatim +*> RCOND is REAL +*> The reciprocal of the condition number of the matrix A, +*> computed as RCOND = 1/(ANORM * AINVNM), where AINVNM is an +*> estimate of the 1-norm of inv(A) computed in this routine. +*> \endverbatim +*> +*> \param[out] WORK +*> \verbatim +*> WORK is COMPLEX array, dimension (2*N) +*> \endverbatim +*> +*> \param[out] IWORK +*> \verbatim +*> IWORK is INTEGER array, dimension (N) +*> \endverbatim +*> +*> \param[out] INFO +*> \verbatim +*> INFO is INTEGER +*> = 0: successful exit +*> < 0: if INFO = -i, the i-th argument had an illegal value +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2011 +* +*> \ingroup complexSYcomputational +* +*> \par Contributors: +* ================== +*> \verbatim +*> +*> November 2011, Igor Kozachenko, +*> Computer Science Division, +*> University of California, Berkeley +*> +*> September 2007, Sven Hammarling, Nicholas J. Higham, Craig Lucas, +*> School of Mathematics, +*> University of Manchester +*> +*> \endverbatim +* +* ===================================================================== + SUBROUTINE CSYCON_ROOK( UPLO, N, A, LDA, IPIV, ANORM, RCOND, WORK, + $ IWORK, INFO ) +* +* -- LAPACK computational routine (version 3.4.0) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2011 +* +* .. Scalar Arguments .. + CHARACTER UPLO + INTEGER INFO, LDA, N + REAL ANORM, RCOND +* .. +* .. Array Arguments .. + INTEGER IPIV( * ), IWORK( * ) + COMPLEX A( LDA, * ), WORK( * ) +* .. +* +* ===================================================================== +* +* .. Parameters .. + REAL ONE, ZERO + PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) + COMPLEX CZERO + PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ) ) +* .. +* .. Local Scalars .. + LOGICAL UPPER + INTEGER I, KASE + REAL AINVNM +* .. +* .. Local Arrays .. + INTEGER ISAVE( 3 ) +* .. +* .. External Functions .. + LOGICAL LSAME + EXTERNAL LSAME +* .. +* .. External Subroutines .. + EXTERNAL CLACN2, CSYTRS_ROOK, XERBLA +* .. +* .. Intrinsic Functions .. + INTRINSIC MAX +* .. +* .. Executable Statements .. +* +* Test the input parameters. +* + INFO = 0 + UPPER = LSAME( UPLO, 'U' ) + IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN + INFO = -1 + ELSE IF( N.LT.0 ) THEN + INFO = -2 + ELSE IF( LDA.LT.MAX( 1, N ) ) THEN + INFO = -4 + ELSE IF( ANORM.LT.ZERO ) THEN + INFO = -6 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'CSYCON_ROOK', -INFO ) + RETURN + END IF +* +* Quick return if possible +* + RCOND = ZERO + IF( N.EQ.0 ) THEN + RCOND = ONE + RETURN + ELSE IF( ANORM.LE.ZERO ) THEN + RETURN + END IF +* +* Check that the diagonal matrix D is nonsingular. +* + IF( UPPER ) THEN +* +* Upper triangular storage: examine D from bottom to top +* + DO 10 I = N, 1, -1 + IF( IPIV( I ).GT.0 .AND. A( I, I ).EQ.CZERO ) + $ RETURN + 10 CONTINUE + ELSE +* +* Lower triangular storage: examine D from top to bottom. +* + DO 20 I = 1, N + IF( IPIV( I ).GT.0 .AND. A( I, I ).EQ.CZERO ) + $ RETURN + 20 CONTINUE + END IF +* +* Estimate the 1-norm of the inverse. +* + KASE = 0 + 30 CONTINUE + CALL CLACN2( N, WORK( N+1 ), WORK, IWORK, AINVNM, KASE, ISAVE ) + IF( KASE.NE.0 ) THEN +* +* Multiply by inv(L*D*L**T) or inv(U*D*U**T). +* + CALL CSYTRS_ROOK( UPLO, N, 1, A, LDA, IPIV, WORK, N, INFO ) + GO TO 30 + END IF +* +* Compute the estimate of the reciprocal condition number. +* + IF( AINVNM.NE.ZERO ) + $ RCOND = ( ONE / AINVNM ) / ANORM +* + RETURN +* +* End of CSYCON_ROOK +* + END diff --git a/SRC/csysv_rook.f b/SRC/csysv_rook.f new file mode 100644 index 00000000..efac6684 --- /dev/null +++ b/SRC/csysv_rook.f @@ -0,0 +1,293 @@ +*> \brief <b> CSYSV_ROOK computes the solution to system of linear equations A * X = B for SY matrices</b> +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download CSYSV_ROOK + dependencies +*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/csysv_rook.f"> +*> [TGZ]</a> +*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/csysv_rook.f"> +*> [ZIP]</a> +*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/csysv_rook.f"> +*> [TXT]</a> +*> \endhtmlonly +* +* Definition: +* =========== +* +* SUBROUTINE CSYSV_ROOK( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, WORK, +* LWORK, INFO ) +* +* .. Scalar Arguments .. +* CHARACTER UPLO +* INTEGER INFO, LDA, LDB, LWORK, N, NRHS +* .. +* .. Array Arguments .. +* INTEGER IPIV( * ) +* COMPLEX A( LDA, * ), B( LDB, * ), WORK( * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> CSYSV_ROOK computes the solution to a complex system of linear +*> equations +*> A * X = B, +*> where A is an N-by-N symmetric matrix and X and B are N-by-NRHS +*> matrices. +*> +*> The diagonal pivoting method is used to factor A as +*> A = U * D * U**T, if UPLO = 'U', or +*> A = L * D * L**T, if UPLO = 'L', +*> where U (or L) is a product of permutation and unit upper (lower) +*> triangular matrices, and D is symmetric and block diagonal with +*> 1-by-1 and 2-by-2 diagonal blocks. +*> +*> CSYTRF_ROOK is called to compute the factorization of a complex +*> symmetric matrix A using the bounded Bunch-Kaufman ("rook") diagonal +*> pivoting method. +*> +*> The factored form of A is then used to solve the system +*> of equations A * X = B by calling CSYTRS_ROOK. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] UPLO +*> \verbatim +*> UPLO is CHARACTER*1 +*> = 'U': Upper triangle of A is stored; +*> = 'L': Lower triangle of A is stored. +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The number of linear equations, i.e., the order of the +*> matrix A. N >= 0. +*> \endverbatim +*> +*> \param[in] NRHS +*> \verbatim +*> NRHS is INTEGER +*> The number of right hand sides, i.e., the number of columns +*> of the matrix B. NRHS >= 0. +*> \endverbatim +*> +*> \param[in,out] A +*> \verbatim +*> A is COMPLEX array, dimension (LDA,N) +*> On entry, the symmetric matrix A. If UPLO = 'U', the leading +*> N-by-N upper triangular part of A contains the upper +*> triangular part of the matrix A, and the strictly lower +*> triangular part of A is not referenced. If UPLO = 'L', the +*> leading N-by-N lower triangular part of A contains the lower +*> triangular part of the matrix A, and the strictly upper +*> triangular part of A is not referenced. +*> +*> On exit, if INFO = 0, the block diagonal matrix D and the +*> multipliers used to obtain the factor U or L from the +*> factorization A = U*D*U**T or A = L*D*L**T as computed by +*> CSYTRF_ROOK. +*> \endverbatim +*> +*> \param[in] LDA +*> \verbatim +*> LDA is INTEGER +*> The leading dimension of the array A. LDA >= max(1,N). +*> \endverbatim +*> +*> \param[out] IPIV +*> \verbatim +*> IPIV is INTEGER array, dimension (N) +*> Details of the interchanges and the block structure of D, +*> as determined by CSYTRF_ROOK. +*> +*> If UPLO = 'U': +*> If IPIV(k) > 0, then rows and columns k and IPIV(k) +*> were interchanged and D(k,k) is a 1-by-1 diagonal block. +*> +*> If IPIV(k) < 0 and IPIV(k-1) < 0, then rows and +*> columns k and -IPIV(k) were interchanged and rows and +*> columns k-1 and -IPIV(k-1) were inerchaged, +*> D(k-1:k,k-1:k) is a 2-by-2 diagonal block. +*> +*> If UPLO = 'L': +*> If IPIV(k) > 0, then rows and columns k and IPIV(k) +*> were interchanged and D(k,k) is a 1-by-1 diagonal block. +*> +*> If IPIV(k) < 0 and IPIV(k+1) < 0, then rows and +*> columns k and -IPIV(k) were interchanged and rows and +*> columns k+1 and -IPIV(k+1) were inerchaged, +*> D(k:k+1,k:k+1) is a 2-by-2 diagonal block. +*> \endverbatim +*> +*> \param[in,out] B +*> \verbatim +*> B is COMPLEX array, dimension (LDB,NRHS) +*> On entry, the N-by-NRHS right hand side matrix B. +*> On exit, if INFO = 0, the N-by-NRHS solution matrix X. +*> \endverbatim +*> +*> \param[in] LDB +*> \verbatim +*> LDB is INTEGER +*> The leading dimension of the array B. LDB >= max(1,N). +*> \endverbatim +*> +*> \param[out] WORK +*> \verbatim +*> WORK is COMPLEX array, dimension (MAX(1,LWORK)) +*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK. +*> \endverbatim +*> +*> \param[in] LWORK +*> \verbatim +*> LWORK is INTEGER +*> The length of WORK. LWORK >= 1, and for best performance +*> LWORK >= max(1,N*NB), where NB is the optimal blocksize for +*> CSYTRF_ROOK. +*> +*> TRS will be done with Level 2 BLAS +*> +*> If LWORK = -1, then a workspace query is assumed; the routine +*> only calculates the optimal size of the WORK array, returns +*> this value as the first entry of the WORK array, and no error +*> message related to LWORK is issued by XERBLA. +*> \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, so the solution could not be computed. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2011 +* +*> \ingroup complexSYsolve +* +*> \par Contributors: +* ================== +*> +*> \verbatim +*> +*> November 2011, Igor Kozachenko, +*> Computer Science Division, +*> University of California, Berkeley +*> +*> September 2007, Sven Hammarling, Nicholas J. Higham, Craig Lucas, +*> School of Mathematics, +*> University of Manchester +*> +*> \endverbatim +* +* ===================================================================== + SUBROUTINE CSYSV_ROOK( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, WORK, + $ LWORK, INFO ) +* +* -- LAPACK driver routine (version 3.4.0) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2011 +* +* .. Scalar Arguments .. + CHARACTER UPLO + INTEGER INFO, LDA, LDB, LWORK, N, NRHS +* .. +* .. Array Arguments .. + INTEGER IPIV( * ) + COMPLEX A( LDA, * ), B( LDB, * ), WORK( * ) +* .. +* +* ===================================================================== +* +* .. Local Scalars .. + LOGICAL LQUERY + INTEGER LWKOPT +* .. +* .. External Functions .. + LOGICAL LSAME + EXTERNAL LSAME +* .. +* .. External Subroutines .. + EXTERNAL XERBLA, CSYTRF_ROOK, CSYTRS_ROOK +* .. +* .. Intrinsic Functions .. + INTRINSIC MAX +* .. +* .. Executable Statements .. +* +* Test the input parameters. +* + INFO = 0 + LQUERY = ( LWORK.EQ.-1 ) + IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN + INFO = -1 + ELSE IF( N.LT.0 ) THEN + INFO = -2 + ELSE IF( NRHS.LT.0 ) THEN + INFO = -3 + ELSE IF( LDA.LT.MAX( 1, N ) ) THEN + INFO = -5 + ELSE IF( LDB.LT.MAX( 1, N ) ) THEN + INFO = -8 + ELSE IF( LWORK.LT.1 .AND. .NOT.LQUERY ) THEN + INFO = -10 + END IF +* + IF( INFO.EQ.0 ) THEN + IF( N.EQ.0 ) THEN + LWKOPT = 1 + ELSE + CALL DSYTRF( UPLO, N, A, LDA, IPIV, WORK, -1, INFO ) + LWKOPT = WORK(1) + END IF + WORK( 1 ) = LWKOPT + END IF +* + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'CSYSV_ROOK ', -INFO ) + RETURN + ELSE IF( LQUERY ) THEN + RETURN + END IF +* +* Compute the factorization A = U*D*U**T or A = L*D*L**T. +* + CALL CSYTRF_ROOK( UPLO, N, A, LDA, IPIV, WORK, LWORK, INFO ) + IF( INFO.EQ.0 ) THEN +* +* Solve the system A*X = B, overwriting B with X. +* +* Solve with TRS_ROOK ( Use Level 2 BLAS) +* + CALL CSYTRS_ROOK( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, INFO ) +* + END IF +* + WORK( 1 ) = LWKOPT +* + RETURN +* +* End of CSYSV_ROOK +* + END diff --git a/SRC/csytf2_rook.f b/SRC/csytf2_rook.f new file mode 100644 index 00000000..50b9c952 --- /dev/null +++ b/SRC/csytf2_rook.f @@ -0,0 +1,824 @@ +*> \brief \b CSYTF2_ROOK +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download CSYTF2_ROOK + dependencies +*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/csytf2_rook.f"> +*> [TGZ]</a> +*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/csytf2_rook.f"> +*> [ZIP]</a> +*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/csytf2_rook.f"> +*> [TXT]</a> +*> \endhtmlonly +* +* Definition: +* =========== +* +* SUBROUTINE CSYTF2_ROOK( UPLO, N, A, LDA, IPIV, INFO ) +* +* .. Scalar Arguments .. +* CHARACTER UPLO +* INTEGER INFO, LDA, N +* .. +* .. Array Arguments .. +* INTEGER IPIV( * ) +* COMPLEX A( LDA, * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> CSYTF2_ROOK computes the factorization of a complex symmetric matrix A +*> using the bounded Bunch-Kaufman ("rook") diagonal pivoting method: +*> +*> A = U*D*U**T or A = L*D*L**T +*> +*> where U (or L) is a product of permutation and unit upper (lower) +*> triangular matrices, U**T is the transpose of U, and D is symmetric and +*> block diagonal with 1-by-1 and 2-by-2 diagonal blocks. +*> +*> This is the unblocked version of the algorithm, calling Level 2 BLAS. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] UPLO +*> \verbatim +*> UPLO is CHARACTER*1 +*> Specifies whether the upper or lower triangular part of the +*> symmetric matrix A is stored: +*> = 'U': Upper triangular +*> = 'L': Lower triangular +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The order of the matrix A. N >= 0. +*> \endverbatim +*> +*> \param[in,out] A +*> \verbatim +*> A is COMPLEX array, dimension (LDA,N) +*> On entry, the symmetric matrix A. If UPLO = 'U', the leading +*> n-by-n upper triangular part of A contains the upper +*> triangular part of the matrix A, and the strictly lower +*> triangular part of A is not referenced. If UPLO = 'L', the +*> leading n-by-n lower triangular part of A contains the lower +*> triangular part of the matrix A, and the strictly upper +*> triangular part of A is not referenced. +*> +*> On exit, the block diagonal matrix D and the multipliers used +*> to obtain the factor U or L (see below for further details). +*> \endverbatim +*> +*> \param[in] LDA +*> \verbatim +*> LDA is INTEGER +*> The leading dimension of the array A. LDA >= max(1,N). +*> \endverbatim +*> +*> \param[out] IPIV +*> \verbatim +*> IPIV is INTEGER array, dimension (N) +*> Details of the interchanges and the block structure of D. +*> +*> If UPLO = 'U': +*> If IPIV(k) > 0, then rows and columns k and IPIV(k) +*> were interchanged and D(k,k) is a 1-by-1 diagonal block. +*> +*> If IPIV(k) < 0 and IPIV(k-1) < 0, then rows and +*> columns k and -IPIV(k) were interchanged and rows and +*> columns k-1 and -IPIV(k-1) were inerchaged, +*> D(k-1:k,k-1:k) is a 2-by-2 diagonal block. +*> +*> If UPLO = 'L': +*> If IPIV(k) > 0, then rows and columns k and IPIV(k) +*> were interchanged and D(k,k) is a 1-by-1 diagonal block. +*> +*> If IPIV(k) < 0 and IPIV(k+1) < 0, then rows and +*> columns k and -IPIV(k) were interchanged and rows and +*> columns k+1 and -IPIV(k+1) were inerchaged, +*> D(k:k+1,k:k+1) is a 2-by-2 diagonal block. +*> \endverbatim +*> +*> \param[out] INFO +*> \verbatim +*> INFO is INTEGER +*> = 0: successful exit +*> < 0: if INFO = -k, the k-th argument had an illegal value +*> > 0: if INFO = k, D(k,k) 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: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2011 +* +*> \ingroup complexSYcomputational +* +*> \par Further Details: +* ===================== +*> +*> \verbatim +*> +*> If UPLO = 'U', then A = U*D*U**T, where +*> U = P(n)*U(n)* ... *P(k)U(k)* ..., +*> i.e., U is a product of terms P(k)*U(k), where k decreases from n to +*> 1 in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1 +*> and 2-by-2 diagonal blocks D(k). P(k) is a permutation matrix as +*> defined by IPIV(k), and U(k) is a unit upper triangular matrix, such +*> that if the diagonal block D(k) is of order s (s = 1 or 2), then +*> +*> ( I v 0 ) k-s +*> U(k) = ( 0 I 0 ) s +*> ( 0 0 I ) n-k +*> k-s s n-k +*> +*> If s = 1, D(k) overwrites A(k,k), and v overwrites A(1:k-1,k). +*> If s = 2, the upper triangle of D(k) overwrites A(k-1,k-1), A(k-1,k), +*> and A(k,k), and v overwrites A(1:k-2,k-1:k). +*> +*> If UPLO = 'L', then A = L*D*L**T, where +*> L = P(1)*L(1)* ... *P(k)*L(k)* ..., +*> i.e., L is a product of terms P(k)*L(k), where k increases from 1 to +*> n in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1 +*> and 2-by-2 diagonal blocks D(k). P(k) is a permutation matrix as +*> defined by IPIV(k), and L(k) is a unit lower triangular matrix, such +*> that if the diagonal block D(k) is of order s (s = 1 or 2), then +*> +*> ( I 0 0 ) k-1 +*> L(k) = ( 0 I 0 ) s +*> ( 0 v I ) n-k-s+1 +*> k-1 s n-k-s+1 +*> +*> If s = 1, D(k) overwrites A(k,k), and v overwrites A(k+1:n,k). +*> If s = 2, the lower triangle of D(k) overwrites A(k,k), A(k+1,k), +*> and A(k+1,k+1), and v overwrites A(k+2:n,k:k+1). +*> \endverbatim +* +*> \par Contributors: +* ================== +*> +*> \verbatim +*> +*> November 2011, Igor Kozachenko, +*> Computer Science Division, +*> University of California, Berkeley +*> +*> September 2007, Sven Hammarling, Nicholas J. Higham, Craig Lucas, +*> School of Mathematics, +*> University of Manchester +*> +*> 01-01-96 - Based on modifications by +*> J. Lewis, Boeing Computer Services Company +*> A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA +*> 1-96 - Based on modifications by J. Lewis, Boeing Computer Services +*> Company +*> \endverbatim +* +* ===================================================================== + SUBROUTINE CSYTF2_ROOK( UPLO, N, A, LDA, IPIV, INFO ) +* +* -- LAPACK computational routine (version 3.3.1) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2011 +* +* .. Scalar Arguments .. + CHARACTER UPLO + INTEGER INFO, LDA, N +* .. +* .. Array Arguments .. + INTEGER IPIV( * ) + COMPLEX A( LDA, * ) +* .. +* +* ===================================================================== +* +* .. Parameters .. + REAL ZERO, ONE + PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) + REAL EIGHT, SEVTEN + PARAMETER ( EIGHT = 8.0E+0, SEVTEN = 17.0E+0 ) + COMPLEX CONE + PARAMETER ( CONE = ( 1.0E+0, 0.0E+0 ) ) +* .. +* .. Local Scalars .. + LOGICAL UPPER, DONE + INTEGER I, IMAX, J, JMAX, ITEMP, K, KK, KP, KSTEP, + $ P, II + REAL ABSAKK, ALPHA, COLMAX, ROWMAX, STEMP, SFMIN + COMPLEX D11, D12, D21, D22, + $ T, WK, WKM1, WKP1, Z +* .. +* .. External Functions .. + LOGICAL LSAME + INTEGER ICAMAX + REAL SLAMCH + EXTERNAL LSAME, ICAMAX, SLAMCH +* .. +* .. External Subroutines .. + EXTERNAL CSCAL, CSWAP, CSYR, XERBLA +* .. +* .. Intrinsic Functions .. + INTRINSIC ABS, MAX, SQRT, AIMAG, REAL +* .. +* .. Statement Functions .. + REAL CABS1 +* .. +* .. Statement Function definitions .. + CABS1( Z ) = ABS( REAL( Z ) ) + ABS( AIMAG( Z ) ) +* .. +* .. Executable Statements .. +* +* Test the input parameters. +* + INFO = 0 + UPPER = LSAME( UPLO, 'U' ) + IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN + INFO = -1 + ELSE IF( N.LT.0 ) THEN + INFO = -2 + ELSE IF( LDA.LT.MAX( 1, N ) ) THEN + INFO = -4 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'CSYTF2_ROOK', -INFO ) + RETURN + END IF +* +* Initialize ALPHA for use in choosing pivot block size. +* + ALPHA = ( ONE+SQRT( SEVTEN ) ) / EIGHT +* +* Compute machine safe minimum +* + SFMIN = SLAMCH( 'S' ) +* + IF( UPPER ) THEN +* +* Factorize A as U*D*U**T using the upper triangle of A +* +* K is the main loop index, decreasing from N to 1 in steps of +* 1 or 2 +* + K = N + 10 CONTINUE +* +* If K < 1, exit from loop +* + IF( K.LT.1 ) + $ GO TO 70 + KSTEP = 1 + P = K +* +* Determine rows and columns to be interchanged and whether +* a 1-by-1 or 2-by-2 pivot block will be used +* + ABSAKK = CABS1( A( K, K ) ) +* +* IMAX is the row-index of the largest off-diagonal element in +* column K, and COLMAX is its absolute value. +* Determine both COLMAX and IMAX. +* + IF( K.GT.1 ) THEN + IMAX = ICAMAX( K-1, A( 1, K ), 1 ) + COLMAX = CABS1( A( IMAX, K ) ) + ELSE + COLMAX = ZERO + END IF +* + IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) ) THEN +* +* Column K is zero: set INFO and continue +* + IF( INFO.EQ.0 ) + $ INFO = K + KP = K + ELSE +* +* Test for interchange +* +* Equivalent to testing for (used to handle NaN and Inf) +* ABSAKK.GE.ALPHA*COLMAX +* + IF( .NOT.( ABSAKK.LT.ALPHA*COLMAX ) ) THEN +* +* no interchange, +* use 1-by-1 pivot block +* + KP = K + ELSE +* + DONE = .FALSE. +* +* Loop until pivot found +* + 12 CONTINUE +* +* Begin pivot search loop body +* +* JMAX is the column-index of the largest off-diagonal +* element in row IMAX, and ROWMAX is its absolute value. +* Determine both ROWMAX and JMAX. +* + IF( IMAX.NE.K ) THEN + JMAX = IMAX + ICAMAX( K-IMAX, A( IMAX, IMAX+1 ), + $ LDA ) + ROWMAX = CABS1( A( IMAX, JMAX ) ) + ELSE + ROWMAX = ZERO + END IF +* + IF( IMAX.GT.1 ) THEN + ITEMP = ICAMAX( IMAX-1, A( 1, IMAX ), 1 ) + STEMP = CABS1( A( ITEMP, IMAX ) ) + IF( STEMP.GT.ROWMAX ) THEN + ROWMAX = STEMP + JMAX = ITEMP + END IF + END IF +* +* Equivalent to testing for (used to handle NaN and Inf) +* CABS1( A( IMAX, IMAX ) ).GE.ALPHA*ROWMAX +* + IF( .NOT.( CABS1(A( IMAX, IMAX )).LT.ALPHA*ROWMAX ) ) + $ THEN +* +* interchange rows and columns K and IMAX, +* use 1-by-1 pivot block +* + KP = IMAX + DONE = .TRUE. +* +* Equivalent to testing for ROWMAX .EQ. COLMAX, +* used to handle NaN and Inf +* + ELSE IF( ( P.EQ.JMAX ).OR.( ROWMAX.LE.COLMAX ) ) THEN +* +* interchange rows and columns K+1 and IMAX, +* use 2-by-2 pivot block +* + KP = IMAX + KSTEP = 2 + DONE = .TRUE. + ELSE +* +* Pivot NOT found, set variables and repeat +* + P = IMAX + COLMAX = ROWMAX + IMAX = JMAX + END IF +* +* End pivot search loop body +* + IF( .NOT. DONE ) GOTO 12 +* + END IF +* +* Swap TWO rows and TWO columns +* +* First swap +* + IF( ( KSTEP.EQ.2 ) .AND. ( P.NE.K ) ) THEN +* +* Interchange rows and column K and P in the leading +* submatrix A(1:k,1:k) if we have a 2-by-2 pivot +* + IF( P.GT.1 ) + $ CALL CSWAP( P-1, A( 1, K ), 1, A( 1, P ), 1 ) + IF( P.LT.(K-1) ) + $ CALL CSWAP( K-P-1, A( P+1, K ), 1, A( P, P+1 ), + $ LDA ) + T = A( K, K ) + A( K, K ) = A( P, P ) + A( P, P ) = T + END IF +* +* Second swap +* + KK = K - KSTEP + 1 + IF( KP.NE.KK ) THEN +* +* Interchange rows and columns KK and KP in the leading +* submatrix A(1:k,1:k) +* + IF( KP.GT.1 ) + $ CALL CSWAP( KP-1, A( 1, KK ), 1, A( 1, KP ), 1 ) + IF( ( KK.GT.1 ) .AND. ( KP.LT.(KK-1) ) ) + $ CALL CSWAP( KK-KP-1, A( KP+1, KK ), 1, A( KP, KP+1 ), + $ LDA ) + T = A( KK, KK ) + A( KK, KK ) = A( KP, KP ) + A( KP, KP ) = T + IF( KSTEP.EQ.2 ) THEN + T = A( K-1, K ) + A( K-1, K ) = A( KP, K ) + A( KP, K ) = T + END IF + END IF +* +* Update the leading submatrix +* + IF( KSTEP.EQ.1 ) THEN +* +* 1-by-1 pivot block D(k): column k now holds +* +* W(k) = U(k)*D(k) +* +* where U(k) is the k-th column of U +* + IF( K.GT.1 ) THEN +* +* Perform a rank-1 update of A(1:k-1,1:k-1) and +* store U(k) in column k +* + IF( CABS1( A( K, K ) ).GE.SFMIN ) THEN +* +* Perform a rank-1 update of A(1:k-1,1:k-1) as +* A := A - U(k)*D(k)*U(k)**T +* = A - W(k)*1/D(k)*W(k)**T +* + D11 = CONE / A( K, K ) + CALL CSYR( UPLO, K-1, -D11, A( 1, K ), 1, A, LDA ) +* +* Store U(k) in column k +* + CALL CSCAL( K-1, D11, A( 1, K ), 1 ) + ELSE +* +* Store L(k) in column K +* + D11 = A( K, K ) + DO 16 II = 1, K - 1 + A( II, K ) = A( II, K ) / D11 + 16 CONTINUE +* +* Perform a rank-1 update of A(k+1:n,k+1:n) as +* A := A - U(k)*D(k)*U(k)**T +* = A - W(k)*(1/D(k))*W(k)**T +* = A - (W(k)/D(k))*(D(k))*(W(k)/D(K))**T +* + CALL CSYR( UPLO, K-1, -D11, A( 1, K ), 1, A, LDA ) + END IF + END IF +* + ELSE +* +* 2-by-2 pivot block D(k): columns k and k-1 now hold +* +* ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k) +* +* where U(k) and U(k-1) are the k-th and (k-1)-th columns +* of U +* +* Perform a rank-2 update of A(1:k-2,1:k-2) as +* +* A := A - ( U(k-1) U(k) )*D(k)*( U(k-1) U(k) )**T +* = A - ( ( A(k-1)A(k) )*inv(D(k)) ) * ( A(k-1)A(k) )**T +* +* and store L(k) and L(k+1) in columns k and k+1 +* + IF( K.GT.2 ) THEN +* + D12 = A( K-1, K ) + D22 = A( K-1, K-1 ) / D12 + D11 = A( K, K ) / D12 + T = CONE / ( D11*D22-CONE ) +* + DO 30 J = K - 2, 1, -1 +* + WKM1 = T*( D11*A( J, K-1 )-A( J, K ) ) + WK = T*( D22*A( J, K )-A( J, K-1 ) ) +* + DO 20 I = J, 1, -1 + A( I, J ) = A( I, J ) - (A( I, K ) / D12 )*WK - + $ ( A( I, K-1 ) / D12 )*WKM1 + 20 CONTINUE +* +* Store U(k) and U(k-1) in cols k and k-1 for row J +* + A( J, K ) = WK / D12 + A( J, K-1 ) = WKM1 / D12 +* + 30 CONTINUE +* + END IF +* + END IF + END IF +* +* Store details of the interchanges in IPIV +* + IF( KSTEP.EQ.1 ) THEN + IPIV( K ) = KP + ELSE + IPIV( K ) = -P + IPIV( K-1 ) = -KP + END IF +* +* Decrease K and return to the start of the main loop +* + K = K - KSTEP + GO TO 10 +* + ELSE +* +* Factorize A as L*D*L**T using the lower triangle of A +* +* K is the main loop index, increasing from 1 to N in steps of +* 1 or 2 +* + K = 1 + 40 CONTINUE +* +* If K > N, exit from loop +* + IF( K.GT.N ) + $ GO TO 70 + KSTEP = 1 + P = K +* +* Determine rows and columns to be interchanged and whether +* a 1-by-1 or 2-by-2 pivot block will be used +* + ABSAKK = CABS1( A( K, K ) ) +* +* IMAX is the row-index of the largest off-diagonal element in +* column K, and COLMAX is its absolute value. +* Determine both COLMAX and IMAX. +* + IF( K.LT.N ) THEN + IMAX = K + ICAMAX( N-K, A( K+1, K ), 1 ) + COLMAX = CABS1( A( IMAX, K ) ) + ELSE + COLMAX = ZERO + END IF +* + IF( ( MAX( ABSAKK, COLMAX ).EQ.ZERO ) ) THEN +* +* Column K is zero: set INFO and continue +* + IF( INFO.EQ.0 ) + $ INFO = K + KP = K + ELSE +* +* Test for interchange +* +* Equivalent to testing for (used to handle NaN and Inf) +* ABSAKK.GE.ALPHA*COLMAX +* + IF( .NOT.( ABSAKK.LT.ALPHA*COLMAX ) ) THEN +* +* no interchange, use 1-by-1 pivot block +* + KP = K + ELSE +* + DONE = .FALSE. +* +* Loop until pivot found +* + 42 CONTINUE +* +* Begin pivot search loop body +* +* JMAX is the column-index of the largest off-diagonal +* element in row IMAX, and ROWMAX is its absolute value. +* Determine both ROWMAX and JMAX. +* + IF( IMAX.NE.K ) THEN + JMAX = K - 1 + ICAMAX( IMAX-K, A( IMAX, K ), LDA ) + ROWMAX = CABS1( A( IMAX, JMAX ) ) + ELSE + ROWMAX = ZERO + END IF +* + IF( IMAX.LT.N ) THEN + ITEMP = IMAX + ICAMAX( N-IMAX, A( IMAX+1, IMAX ), + $ 1 ) + STEMP = CABS1( A( ITEMP, IMAX ) ) + IF( STEMP.GT.ROWMAX ) THEN + ROWMAX = STEMP + JMAX = ITEMP + END IF + END IF +* +* Equivalent to testing for (used to handle NaN and Inf) +* CABS1( A( IMAX, IMAX ) ).GE.ALPHA*ROWMAX +* + IF( .NOT.( CABS1(A( IMAX, IMAX )).LT.ALPHA*ROWMAX ) ) + $ THEN +* +* interchange rows and columns K and IMAX, +* use 1-by-1 pivot block +* + KP = IMAX + DONE = .TRUE. +* +* Equivalent to testing for ROWMAX .EQ. COLMAX, +* used to handle NaN and Inf +* + ELSE IF( ( P.EQ.JMAX ).OR.( ROWMAX.LE.COLMAX ) ) THEN +* +* interchange rows and columns K+1 and IMAX, +* use 2-by-2 pivot block +* + KP = IMAX + KSTEP = 2 + DONE = .TRUE. + ELSE +* +* Pivot NOT found, set variables and repeat +* + P = IMAX + COLMAX = ROWMAX + IMAX = JMAX + END IF +* +* End pivot search loop body +* + IF( .NOT. DONE ) GOTO 42 +* + END IF +* +* Swap TWO rows and TWO columns +* +* First swap +* + IF( ( KSTEP.EQ.2 ) .AND. ( P.NE.K ) ) THEN +* +* Interchange rows and column K and P in the trailing +* submatrix A(k:n,k:n) if we have a 2-by-2 pivot +* + IF( P.LT.N ) + $ CALL CSWAP( N-P, A( P+1, K ), 1, A( P+1, P ), 1 ) + IF( P.GT.(K+1) ) + $ CALL CSWAP( P-K-1, A( K+1, K ), 1, A( P, K+1 ), LDA ) + T = A( K, K ) + A( K, K ) = A( P, P ) + A( P, P ) = T + END IF +* +* Second swap +* + KK = K + KSTEP - 1 + IF( KP.NE.KK ) THEN +* +* Interchange rows and columns KK and KP in the trailing +* submatrix A(k:n,k:n) +* + IF( KP.LT.N ) + $ CALL CSWAP( N-KP, A( KP+1, KK ), 1, A( KP+1, KP ), 1 ) + IF( ( KK.LT.N ) .AND. ( KP.GT.(KK+1) ) ) + $ CALL CSWAP( KP-KK-1, A( KK+1, KK ), 1, A( KP, KK+1 ), + $ LDA ) + T = A( KK, KK ) + A( KK, KK ) = A( KP, KP ) + A( KP, KP ) = T + IF( KSTEP.EQ.2 ) THEN + T = A( K+1, K ) + A( K+1, K ) = A( KP, K ) + A( KP, K ) = T + END IF + END IF +* +* Update the trailing submatrix +* + IF( KSTEP.EQ.1 ) THEN +* +* 1-by-1 pivot block D(k): column k now holds +* +* W(k) = L(k)*D(k) +* +* where L(k) is the k-th column of L +* + IF( K.LT.N ) THEN +* +* Perform a rank-1 update of A(k+1:n,k+1:n) and +* store L(k) in column k +* + IF( CABS1( A( K, K ) ).GE.SFMIN ) THEN +* +* Perform a rank-1 update of A(k+1:n,k+1:n) as +* A := A - L(k)*D(k)*L(k)**T +* = A - W(k)*(1/D(k))*W(k)**T +* + D11 = CONE / A( K, K ) + CALL CSYR( UPLO, N-K, -D11, A( K+1, K ), 1, + $ A( K+1, K+1 ), LDA ) +* +* Store L(k) in column k +* + CALL CSCAL( N-K, D11, A( K+1, K ), 1 ) + ELSE +* +* Store L(k) in column k +* + D11 = A( K, K ) + DO 46 II = K + 1, N + A( II, K ) = A( II, K ) / D11 + 46 CONTINUE +* +* Perform a rank-1 update of A(k+1:n,k+1:n) as +* A := A - L(k)*D(k)*L(k)**T +* = A - W(k)*(1/D(k))*W(k)**T +* = A - (W(k)/D(k))*(D(k))*(W(k)/D(K))**T +* + CALL CSYR( UPLO, N-K, -D11, A( K+1, K ), 1, + $ A( K+1, K+1 ), LDA ) + END IF + END IF +* + ELSE +* +* 2-by-2 pivot block D(k): columns k and k+1 now hold +* +* ( W(k) W(k+1) ) = ( L(k) L(k+1) )*D(k) +* +* where L(k) and L(k+1) are the k-th and (k+1)-th columns +* of L +* +* +* Perform a rank-2 update of A(k+2:n,k+2:n) as +* +* A := A - ( L(k) L(k+1) ) * D(k) * ( L(k) L(k+1) )**T +* = A - ( ( A(k)A(k+1) )*inv(D(k) ) * ( A(k)A(k+1) )**T +* +* and store L(k) and L(k+1) in columns k and k+1 +* + IF( K.LT.N-1 ) THEN +* + D21 = A( K+1, K ) + D11 = A( K+1, K+1 ) / D21 + D22 = A( K, K ) / D21 + T = CONE / ( D11*D22-CONE ) +* + DO 60 J = K + 2, N +* +* Compute D21 * ( W(k)W(k+1) ) * inv(D(k)) for row J +* + WK = T*( D11*A( J, K )-A( J, K+1 ) ) + WKP1 = T*( D22*A( J, K+1 )-A( J, K ) ) +* +* Perform a rank-2 update of A(k+2:n,k+2:n) +* + DO 50 I = J, N + A( I, J ) = A( I, J ) - ( A( I, K ) / D21 )*WK - + $ ( A( I, K+1 ) / D21 )*WKP1 + 50 CONTINUE +* +* Store L(k) and L(k+1) in cols k and k+1 for row J +* + A( J, K ) = WK / D21 + A( J, K+1 ) = WKP1 / D21 +* + 60 CONTINUE +* + END IF +* + END IF + END IF +* +* Store details of the interchanges in IPIV +* + IF( KSTEP.EQ.1 ) THEN + IPIV( K ) = KP + ELSE + IPIV( K ) = -P + IPIV( K+1 ) = -KP + END IF +* +* Increase K and return to the start of the main loop +* + K = K + KSTEP + GO TO 40 +* + END IF +* + 70 CONTINUE +* + RETURN +* +* End of CSYTF2_ROOK +* + END diff --git a/SRC/csytrf_rook.f b/SRC/csytrf_rook.f new file mode 100644 index 00000000..d4169d74 --- /dev/null +++ b/SRC/csytrf_rook.f @@ -0,0 +1,393 @@ +*> \brief \b CSYTRF_ROOK +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download CSYTRF_ROOK + dependencies +*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/csytrf_rook.f"> +*> [TGZ]</a> +*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/csytrf_rook.f"> +*> [ZIP]</a> +*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/csytrf_rook.f"> +*> [TXT]</a> +*> \endhtmlonly +* +* Definition: +* =========== +* +* SUBROUTINE CSYTRF_ROOK( UPLO, N, A, LDA, IPIV, WORK, LWORK, INFO ) +* +* .. Scalar Arguments .. +* CHARACTER UPLO +* INTEGER INFO, LDA, LWORK, N +* .. +* .. Array Arguments .. +* INTEGER IPIV( * ) +* COMPLEX A( LDA, * ), WORK( * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> CSYTRF_ROOK computes the factorization of a complex symmetric matrix A +*> using the bounded Bunch-Kaufman ("rook") diagonal pivoting method. +*> The form of the factorization is +*> +*> A = U*D*U**T or A = L*D*L**T +*> +*> where U (or L) is a product of permutation and unit upper (lower) +*> triangular matrices, and D is symmetric and block diagonal with +*> 1-by-1 and 2-by-2 diagonal blocks. +*> +*> This is the blocked version of the algorithm, calling Level 3 BLAS. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] UPLO +*> \verbatim +*> UPLO is CHARACTER*1 +*> = 'U': Upper triangle of A is stored; +*> = 'L': Lower triangle of A is stored. +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The order of the matrix A. N >= 0. +*> \endverbatim +*> +*> \param[in,out] A +*> \verbatim +*> A is COMPLEX array, dimension (LDA,N) +*> On entry, the symmetric matrix A. If UPLO = 'U', the leading +*> N-by-N upper triangular part of A contains the upper +*> triangular part of the matrix A, and the strictly lower +*> triangular part of A is not referenced. If UPLO = 'L', the +*> leading N-by-N lower triangular part of A contains the lower +*> triangular part of the matrix A, and the strictly upper +*> triangular part of A is not referenced. +*> +*> On exit, the block diagonal matrix D and the multipliers used +*> to obtain the factor U or L (see below for further details). +*> \endverbatim +*> +*> \param[in] LDA +*> \verbatim +*> LDA is INTEGER +*> The leading dimension of the array A. LDA >= max(1,N). +*> \endverbatim +*> +*> \param[out] IPIV +*> \verbatim +*> IPIV is INTEGER array, dimension (N) +*> Details of the interchanges and the block structure of D. +*> +*> If UPLO = 'U': +*> If IPIV(k) > 0, then rows and columns k and IPIV(k) +*> were interchanged and D(k,k) is a 1-by-1 diagonal block. +*> +*> If IPIV(k) < 0 and IPIV(k-1) < 0, then rows and +*> columns k and -IPIV(k) were interchanged and rows and +*> columns k-1 and -IPIV(k-1) were inerchaged, +*> D(k-1:k,k-1:k) is a 2-by-2 diagonal block. +*> +*> If UPLO = 'L': +*> If IPIV(k) > 0, then rows and columns k and IPIV(k) +*> were interchanged and D(k,k) is a 1-by-1 diagonal block. +*> +*> If IPIV(k) < 0 and IPIV(k+1) < 0, then rows and +*> columns k and -IPIV(k) were interchanged and rows and +*> columns k+1 and -IPIV(k+1) were inerchaged, +*> D(k:k+1,k:k+1) is a 2-by-2 diagonal block. +*> \endverbatim +*> +*> \param[out] WORK +*> \verbatim +*> WORK is COMPLEX array, dimension (MAX(1,LWORK)). +*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK. +*> \endverbatim +*> +*> \param[in] LWORK +*> \verbatim +*> LWORK is INTEGER +*> The length of WORK. LWORK >=1. For best performance +*> LWORK >= N*NB, where NB is the block size returned by ILAENV. +*> +*> If LWORK = -1, then a workspace query is assumed; the routine +*> only calculates the optimal size of the WORK array, returns +*> this value as the first entry of the WORK array, and no error +*> message related to LWORK is issued by XERBLA. +*> \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: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2011 +* +*> \ingroup complexSYcomputational +* +*> \par Further Details: +* ===================== +*> +*> \verbatim +*> +*> If UPLO = 'U', then A = U*D*U**T, where +*> U = P(n)*U(n)* ... *P(k)U(k)* ..., +*> i.e., U is a product of terms P(k)*U(k), where k decreases from n to +*> 1 in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1 +*> and 2-by-2 diagonal blocks D(k). P(k) is a permutation matrix as +*> defined by IPIV(k), and U(k) is a unit upper triangular matrix, such +*> that if the diagonal block D(k) is of order s (s = 1 or 2), then +*> +*> ( I v 0 ) k-s +*> U(k) = ( 0 I 0 ) s +*> ( 0 0 I ) n-k +*> k-s s n-k +*> +*> If s = 1, D(k) overwrites A(k,k), and v overwrites A(1:k-1,k). +*> If s = 2, the upper triangle of D(k) overwrites A(k-1,k-1), A(k-1,k), +*> and A(k,k), and v overwrites A(1:k-2,k-1:k). +*> +*> If UPLO = 'L', then A = L*D*L**T, where +*> L = P(1)*L(1)* ... *P(k)*L(k)* ..., +*> i.e., L is a product of terms P(k)*L(k), where k increases from 1 to +*> n in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1 +*> and 2-by-2 diagonal blocks D(k). P(k) is a permutation matrix as +*> defined by IPIV(k), and L(k) is a unit lower triangular matrix, such +*> that if the diagonal block D(k) is of order s (s = 1 or 2), then +*> +*> ( I 0 0 ) k-1 +*> L(k) = ( 0 I 0 ) s +*> ( 0 v I ) n-k-s+1 +*> k-1 s n-k-s+1 +*> +*> If s = 1, D(k) overwrites A(k,k), and v overwrites A(k+1:n,k). +*> If s = 2, the lower triangle of D(k) overwrites A(k,k), A(k+1,k), +*> and A(k+1,k+1), and v overwrites A(k+2:n,k:k+1). +*> \endverbatim +* +*> \par Contributors: +* ================== +*> +*> \verbatim +*> +*> November 2011, Igor Kozachenko, +*> Computer Science Division, +*> University of California, Berkeley +*> +*> September 2007, Sven Hammarling, Nicholas J. Higham, Craig Lucas, +*> School of Mathematics, +*> University of Manchester +*> +*> \endverbatim +* +* ===================================================================== + SUBROUTINE CSYTRF_ROOK( UPLO, N, A, LDA, IPIV, WORK, LWORK, INFO ) +* +* -- LAPACK computational routine (version 3.4.0) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2011 +* +* .. Scalar Arguments .. + CHARACTER UPLO + INTEGER INFO, LDA, LWORK, N +* .. +* .. Array Arguments .. + INTEGER IPIV( * ) + COMPLEX A( LDA, * ), WORK( * ) +* .. +* +* ===================================================================== +* +* .. Local Scalars .. + LOGICAL LQUERY, UPPER + INTEGER IINFO, IWS, J, K, KB, LDWORK, LWKOPT, NB, NBMIN +* .. +* .. External Functions .. + LOGICAL LSAME + INTEGER ILAENV + EXTERNAL LSAME, ILAENV +* .. +* .. External Subroutines .. + EXTERNAL CLASYF_ROOK, CSYTF2_ROOK, XERBLA +* .. +* .. Intrinsic Functions .. + INTRINSIC MAX +* .. +* .. Executable Statements .. +* +* Test the input parameters. +* + INFO = 0 + UPPER = LSAME( UPLO, 'U' ) + LQUERY = ( LWORK.EQ.-1 ) + IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN + INFO = -1 + ELSE IF( N.LT.0 ) THEN + INFO = -2 + ELSE IF( LDA.LT.MAX( 1, N ) ) THEN + INFO = -4 + ELSE IF( LWORK.LT.1 .AND. .NOT.LQUERY ) THEN + INFO = -7 + END IF +* + IF( INFO.EQ.0 ) THEN +* +* Determine the block size +* + NB = ILAENV( 1, 'CSYTRF_ROOK', UPLO, N, -1, -1, -1 ) + LWKOPT = N*NB + WORK( 1 ) = LWKOPT + END IF +* + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'CSYTRF_ROOK', -INFO ) + RETURN + ELSE IF( LQUERY ) THEN + RETURN + END IF +* + NBMIN = 2 + LDWORK = N + IF( NB.GT.1 .AND. NB.LT.N ) THEN + IWS = LDWORK*NB + IF( LWORK.LT.IWS ) THEN + NB = MAX( LWORK / LDWORK, 1 ) + NBMIN = MAX( 2, ILAENV( 2, 'CSYTRF_ROOK', + $ UPLO, N, -1, -1, -1 ) ) + END IF + ELSE + IWS = 1 + END IF + IF( NB.LT.NBMIN ) + $ NB = N +* + IF( UPPER ) THEN +* +* Factorize A as U*D*U**T using the upper triangle of A +* +* K is the main loop index, decreasing from N to 1 in steps of +* KB, where KB is the number of columns factorized by CLASYF_ROOK; +* KB is either NB or NB-1, or K for the last block +* + K = N + 10 CONTINUE +* +* If K < 1, exit from loop +* + IF( K.LT.1 ) + $ GO TO 40 +* + IF( K.GT.NB ) THEN +* +* Factorize columns k-kb+1:k of A and use blocked code to +* update columns 1:k-kb +* + CALL CLASYF_ROOK( UPLO, K, NB, KB, A, LDA, + $ IPIV, WORK, LDWORK, IINFO ) + ELSE +* +* Use unblocked code to factorize columns 1:k of A +* + CALL CSYTF2_ROOK( UPLO, K, A, LDA, IPIV, IINFO ) + KB = K + END IF +* +* Set INFO on the first occurrence of a zero pivot +* + IF( INFO.EQ.0 .AND. IINFO.GT.0 ) + $ INFO = IINFO +* +* No need to adjust IPIV +* +* Decrease K and return to the start of the main loop +* + K = K - KB + GO TO 10 +* + ELSE +* +* Factorize A as L*D*L**T using the lower triangle of A +* +* K is the main loop index, increasing from 1 to N in steps of +* KB, where KB is the number of columns factorized by CLASYF_ROOK; +* KB is either NB or NB-1, or N-K+1 for the last block +* + K = 1 + 20 CONTINUE +* +* If K > N, exit from loop +* + IF( K.GT.N ) + $ GO TO 40 +* + IF( K.LE.N-NB ) THEN +* +* Factorize columns k:k+kb-1 of A and use blocked code to +* update columns k+kb:n +* + CALL CLASYF_ROOK( UPLO, N-K+1, NB, KB, A( K, K ), LDA, + $ IPIV( K ), WORK, LDWORK, IINFO ) + ELSE +* +* Use unblocked code to factorize columns k:n of A +* + CALL CSYTF2_ROOK( UPLO, N-K+1, A( K, K ), LDA, IPIV( K ), + $ IINFO ) + KB = N - K + 1 + END IF +* +* Set INFO on the first occurrence of a zero pivot +* + IF( INFO.EQ.0 .AND. IINFO.GT.0 ) + $ INFO = IINFO + K - 1 +* +* Adjust IPIV +* + DO 30 J = K, K + KB - 1 + IF( IPIV( J ).GT.0 ) THEN + IPIV( J ) = IPIV( J ) + K - 1 + ELSE + IPIV( J ) = IPIV( J ) - K + 1 + END IF + 30 CONTINUE +* +* Increase K and return to the start of the main loop +* + K = K + KB + GO TO 20 +* + END IF +* + 40 CONTINUE + WORK( 1 ) = LWKOPT + RETURN +* +* End of CSYTRF_ROOK +* + END diff --git a/SRC/csytri_rook.f b/SRC/csytri_rook.f new file mode 100644 index 00000000..0b57a713 --- /dev/null +++ b/SRC/csytri_rook.f @@ -0,0 +1,451 @@ +*> \brief \b CSYTRI_ROOK +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download CSYTRI_ROOK + dependencies +*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/csytri_rook.f"> +*> [TGZ]</a> +*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/csytri_rook.f"> +*> [ZIP]</a> +*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/csytri_rook.f"> +*> [TXT]</a> +*> \endhtmlonly +* +* Definition: +* =========== +* +* SUBROUTINE CSYTRI_ROOK( UPLO, N, A, LDA, IPIV, WORK, INFO ) +* +* .. Scalar Arguments .. +* CHARACTER UPLO +* INTEGER INFO, LDA, N +* .. +* .. Array Arguments .. +* INTEGER IPIV( * ) +* COMPLEX A( LDA, * ), WORK( * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> CSYTRI_ROOK computes the inverse of a complex symmetric +*> matrix A using the factorization A = U*D*U**T or A = L*D*L**T +*> computed by CSYTRF_ROOK. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] UPLO +*> \verbatim +*> UPLO is CHARACTER*1 +*> Specifies whether the details of the factorization are stored +*> as an upper or lower triangular matrix. +*> = 'U': Upper triangular, form is A = U*D*U**T; +*> = 'L': Lower triangular, form is A = L*D*L**T. +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The order of the matrix A. N >= 0. +*> \endverbatim +*> +*> \param[in,out] A +*> \verbatim +*> A is COMPLEX array, dimension (LDA,N) +*> On entry, the block diagonal matrix D and the multipliers +*> used to obtain the factor U or L as computed by CSYTRF_ROOK. +*> +*> On exit, if INFO = 0, the (symmetric) inverse of the original +*> matrix. 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. +*> \endverbatim +*> +*> \param[in] LDA +*> \verbatim +*> LDA is INTEGER +*> The leading dimension of the array A. LDA >= max(1,N). +*> \endverbatim +*> +*> \param[in] IPIV +*> \verbatim +*> IPIV is INTEGER array, dimension (N) +*> Details of the interchanges and the block structure of D +*> as determined by CSYTRF_ROOK. +*> \endverbatim +*> +*> \param[out] WORK +*> \verbatim +*> WORK is COMPLEX array, dimension (N) +*> \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) = 0; the matrix is singular and its +*> inverse could not be computed. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2011 +* +*> \ingroup complexSYcomputational +* +*> \par Contributors: +* ================== +*> +*> \verbatim +*> +*> November 2011, Igor Kozachenko, +*> Computer Science Division, +*> University of California, Berkeley +*> +*> September 2007, Sven Hammarling, Nicholas J. Higham, Craig Lucas, +*> School of Mathematics, +*> University of Manchester +*> +*> \endverbatim +* +* ===================================================================== + SUBROUTINE CSYTRI_ROOK( UPLO, N, A, LDA, IPIV, WORK, INFO ) +* +* -- LAPACK computational routine (version 3.4.0) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2011 +* +* .. Scalar Arguments .. + CHARACTER UPLO + INTEGER INFO, LDA, N +* .. +* .. Array Arguments .. + INTEGER IPIV( * ) + COMPLEX A( LDA, * ), WORK( * ) +* .. +* +* ===================================================================== +* +* .. Parameters .. + COMPLEX CONE, CZERO + PARAMETER ( CONE = ( 1.0E+0, 0.0E+0 ), + $ CZERO = ( 0.0E+0, 0.0E+0 ) ) +* .. +* .. Local Scalars .. + LOGICAL UPPER + INTEGER K, KP, KSTEP + COMPLEX AK, AKKP1, AKP1, D, T, TEMP +* .. +* .. External Functions .. + LOGICAL LSAME + COMPLEX CDOTU + EXTERNAL LSAME, CDOTU +* .. +* .. External Subroutines .. + EXTERNAL CCOPY, CSWAP, CSYMV, XERBLA +* .. +* .. Intrinsic Functions .. + INTRINSIC MAX +* .. +* .. Executable Statements .. +* +* Test the input parameters. +* + INFO = 0 + UPPER = LSAME( UPLO, 'U' ) + IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN + INFO = -1 + ELSE IF( N.LT.0 ) THEN + INFO = -2 + ELSE IF( LDA.LT.MAX( 1, N ) ) THEN + INFO = -4 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'CSYTRI_ROOK', -INFO ) + RETURN + END IF +* +* Quick return if possible +* + IF( N.EQ.0 ) + $ RETURN +* +* Check that the diagonal matrix D is nonsingular. +* + IF( UPPER ) THEN +* +* Upper triangular storage: examine D from bottom to top +* + DO 10 INFO = N, 1, -1 + IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.CZERO ) + $ RETURN + 10 CONTINUE + ELSE +* +* Lower triangular storage: examine D from top to bottom. +* + DO 20 INFO = 1, N + IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.CZERO ) + $ RETURN + 20 CONTINUE + END IF + INFO = 0 +* + IF( UPPER ) THEN +* +* Compute inv(A) from the factorization A = U*D*U**T. +* +* K is the main loop index, increasing from 1 to N in steps of +* 1 or 2, depending on the size of the diagonal blocks. +* + K = 1 + 30 CONTINUE +* +* If K > N, exit from loop. +* + IF( K.GT.N ) + $ GO TO 40 +* + IF( IPIV( K ).GT.0 ) THEN +* +* 1 x 1 diagonal block +* +* Invert the diagonal block. +* + A( K, K ) = CONE / A( K, K ) +* +* Compute column K of the inverse. +* + IF( K.GT.1 ) THEN + CALL CCOPY( K-1, A( 1, K ), 1, WORK, 1 ) + CALL CSYMV( UPLO, K-1, -CONE, A, LDA, WORK, 1, CZERO, + $ A( 1, K ), 1 ) + A( K, K ) = A( K, K ) - CDOTU( K-1, WORK, 1, A( 1, K ), + $ 1 ) + END IF + KSTEP = 1 + ELSE +* +* 2 x 2 diagonal block +* +* Invert the diagonal block. +* + T = A( K, K+1 ) + AK = A( K, K ) / T + AKP1 = A( K+1, K+1 ) / T + AKKP1 = A( K, K+1 ) / T + D = T*( AK*AKP1-CONE ) + A( K, K ) = AKP1 / D + A( K+1, K+1 ) = AK / D + A( K, K+1 ) = -AKKP1 / D +* +* Compute columns K and K+1 of the inverse. +* + IF( K.GT.1 ) THEN + CALL CCOPY( K-1, A( 1, K ), 1, WORK, 1 ) + CALL CSYMV( UPLO, K-1, -CONE, A, LDA, WORK, 1, CZERO, + $ A( 1, K ), 1 ) + A( K, K ) = A( K, K ) - CDOTU( K-1, WORK, 1, A( 1, K ), + $ 1 ) + A( K, K+1 ) = A( K, K+1 ) - + $ CDOTU( K-1, A( 1, K ), 1, A( 1, K+1 ), 1 ) + CALL CCOPY( K-1, A( 1, K+1 ), 1, WORK, 1 ) + CALL CSYMV( UPLO, K-1, -CONE, A, LDA, WORK, 1, CZERO, + $ A( 1, K+1 ), 1 ) + A( K+1, K+1 ) = A( K+1, K+1 ) - + $ CDOTU( K-1, WORK, 1, A( 1, K+1 ), 1 ) + END IF + KSTEP = 2 + END IF +* + IF( KSTEP.EQ.1 ) THEN +* +* Interchange rows and columns K and IPIV(K) in the leading +* submatrix A(1:k+1,1:k+1) +* + KP = IPIV( K ) + IF( KP.NE.K ) THEN + IF( KP.GT.1 ) + $ CALL CSWAP( KP-1, A( 1, K ), 1, A( 1, KP ), 1 ) + CALL CSWAP( K-KP-1, A( KP+1, K ), 1, A( KP, KP+1 ), LDA ) + TEMP = A( K, K ) + A( K, K ) = A( KP, KP ) + A( KP, KP ) = TEMP + END IF + ELSE +* +* Interchange rows and columns K and K+1 with -IPIV(K) and +* -IPIV(K+1)in the leading submatrix A(1:k+1,1:k+1) +* + KP = -IPIV( K ) + IF( KP.NE.K ) THEN + IF( KP.GT.1 ) + $ CALL CSWAP( KP-1, A( 1, K ), 1, A( 1, KP ), 1 ) + CALL CSWAP( K-KP-1, A( KP+1, K ), 1, A( KP, KP+1 ), LDA ) +* + TEMP = A( K, K ) + A( K, K ) = A( KP, KP ) + A( KP, KP ) = TEMP + TEMP = A( K, K+1 ) + A( K, K+1 ) = A( KP, K+1 ) + A( KP, K+1 ) = TEMP + END IF +* + K = K + 1 + KP = -IPIV( K ) + IF( KP.NE.K ) THEN + IF( KP.GT.1 ) + $ CALL CSWAP( KP-1, A( 1, K ), 1, A( 1, KP ), 1 ) + CALL CSWAP( K-KP-1, A( KP+1, K ), 1, A( KP, KP+1 ), LDA ) + TEMP = A( K, K ) + A( K, K ) = A( KP, KP ) + A( KP, KP ) = TEMP + END IF + END IF +* + K = K + 1 + GO TO 30 + 40 CONTINUE +* + ELSE +* +* Compute inv(A) from the factorization A = L*D*L**T. +* +* K is the main loop index, increasing from 1 to N in steps of +* 1 or 2, depending on the size of the diagonal blocks. +* + K = N + 50 CONTINUE +* +* If K < 1, exit from loop. +* + IF( K.LT.1 ) + $ GO TO 60 +* + IF( IPIV( K ).GT.0 ) THEN +* +* 1 x 1 diagonal block +* +* Invert the diagonal block. +* + A( K, K ) = CONE / A( K, K ) +* +* Compute column K of the inverse. +* + IF( K.LT.N ) THEN + CALL CCOPY( N-K, A( K+1, K ), 1, WORK, 1 ) + CALL CSYMV( UPLO, N-K,-CONE, A( K+1, K+1 ), LDA, WORK, 1, + $ CZERO, A( K+1, K ), 1 ) + A( K, K ) = A( K, K ) - CDOTU( N-K, WORK, 1, A( K+1, K ), + $ 1 ) + END IF + KSTEP = 1 + ELSE +* +* 2 x 2 diagonal block +* +* Invert the diagonal block. +* + T = A( K, K-1 ) + AK = A( K-1, K-1 ) / T + AKP1 = A( K, K ) / T + AKKP1 = A( K, K-1 ) / T + D = T*( AK*AKP1-CONE ) + A( K-1, K-1 ) = AKP1 / D + A( K, K ) = AK / D + A( K, K-1 ) = -AKKP1 / D +* +* Compute columns K-1 and K of the inverse. +* + IF( K.LT.N ) THEN + CALL CCOPY( N-K, A( K+1, K ), 1, WORK, 1 ) + CALL CSYMV( UPLO, N-K,-CONE, A( K+1, K+1 ), LDA, WORK, 1, + $ CZERO, A( K+1, K ), 1 ) + A( K, K ) = A( K, K ) - CDOTU( N-K, WORK, 1, A( K+1, K ), + $ 1 ) + A( K, K-1 ) = A( K, K-1 ) - + $ CDOTU( N-K, A( K+1, K ), 1, A( K+1, K-1 ), + $ 1 ) + CALL CCOPY( N-K, A( K+1, K-1 ), 1, WORK, 1 ) + CALL CSYMV( UPLO, N-K,-CONE, A( K+1, K+1 ), LDA, WORK, 1, + $ CZERO, A( K+1, K-1 ), 1 ) + A( K-1, K-1 ) = A( K-1, K-1 ) - + $ CDOTU( N-K, WORK, 1, A( K+1, K-1 ), 1 ) + END IF + KSTEP = 2 + END IF +* + IF( KSTEP.EQ.1 ) THEN +* +* Interchange rows and columns K and IPIV(K) in the trailing +* submatrix A(k-1:n,k-1:n) +* + KP = IPIV( K ) + IF( KP.NE.K ) THEN + IF( KP.LT.N ) + $ CALL CSWAP( N-KP, A( KP+1, K ), 1, A( KP+1, KP ), 1 ) + CALL CSWAP( KP-K-1, A( K+1, K ), 1, A( KP, K+1 ), LDA ) + TEMP = A( K, K ) + A( K, K ) = A( KP, KP ) + A( KP, KP ) = TEMP + END IF + ELSE +* +* Interchange rows and columns K and K-1 with -IPIV(K) and +* -IPIV(K-1) in the trailing submatrix A(k-1:n,k-1:n) +* + KP = -IPIV( K ) + IF( KP.NE.K ) THEN + IF( KP.LT.N ) + $ CALL CSWAP( N-KP, A( KP+1, K ), 1, A( KP+1, KP ), 1 ) + CALL CSWAP( KP-K-1, A( K+1, K ), 1, A( KP, K+1 ), LDA ) +* + TEMP = A( K, K ) + A( K, K ) = A( KP, KP ) + A( KP, KP ) = TEMP + TEMP = A( K, K-1 ) + A( K, K-1 ) = A( KP, K-1 ) + A( KP, K-1 ) = TEMP + END IF +* + K = K - 1 + KP = -IPIV( K ) + IF( KP.NE.K ) THEN + IF( KP.LT.N ) + $ CALL CSWAP( N-KP, A( KP+1, K ), 1, A( KP+1, KP ), 1 ) + CALL CSWAP( KP-K-1, A( K+1, K ), 1, A( KP, K+1 ), LDA ) + TEMP = A( K, K ) + A( K, K ) = A( KP, KP ) + A( KP, KP ) = TEMP + END IF + END IF +* + K = K - 1 + GO TO 50 + 60 CONTINUE + END IF +* + RETURN +* +* End of CSYTRI_ROOK +* + END diff --git a/SRC/csytrs_rook.f b/SRC/csytrs_rook.f new file mode 100644 index 00000000..44727212 --- /dev/null +++ b/SRC/csytrs_rook.f @@ -0,0 +1,484 @@ +*> \brief \b CSYTRS_ROOK +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download CSYTRS_ROOK + dependencies +*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/csytrs_rook.f"> +*> [TGZ]</a> +*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/csytrs_rook.f"> +*> [ZIP]</a> +*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/csytrs_rook.f"> +*> [TXT]</a> +*> \endhtmlonly +* +* Definition: +* =========== +* +* SUBROUTINE CSYTRS_ROOK( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, INFO ) +* +* .. Scalar Arguments .. +* CHARACTER UPLO +* INTEGER INFO, LDA, LDB, N, NRHS +* .. +* .. Array Arguments .. +* INTEGER IPIV( * ) +* COMPLEX A( LDA, * ), B( LDB, * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> CSYTRS_ROOK solves a system of linear equations A*X = B with +*> a complex symmetric matrix A using the factorization A = U*D*U**T or +*> A = L*D*L**T computed by CSYTRF_ROOK. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] UPLO +*> \verbatim +*> UPLO is CHARACTER*1 +*> Specifies whether the details of the factorization are stored +*> as an upper or lower triangular matrix. +*> = 'U': Upper triangular, form is A = U*D*U**T; +*> = 'L': Lower triangular, form is A = L*D*L**T. +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The order of the matrix A. N >= 0. +*> \endverbatim +*> +*> \param[in] NRHS +*> \verbatim +*> NRHS is INTEGER +*> The number of right hand sides, i.e., the number of columns +*> of the matrix B. NRHS >= 0. +*> \endverbatim +*> +*> \param[in] A +*> \verbatim +*> A is COMPLEX array, dimension (LDA,N) +*> The block diagonal matrix D and the multipliers used to +*> obtain the factor U or L as computed by CSYTRF_ROOK. +*> \endverbatim +*> +*> \param[in] LDA +*> \verbatim +*> LDA is INTEGER +*> The leading dimension of the array A. LDA >= max(1,N). +*> \endverbatim +*> +*> \param[in] IPIV +*> \verbatim +*> IPIV is INTEGER array, dimension (N) +*> Details of the interchanges and the block structure of D +*> as determined by CSYTRF_ROOK. +*> \endverbatim +*> +*> \param[in,out] B +*> \verbatim +*> B is COMPLEX array, dimension (LDB,NRHS) +*> On entry, the right hand side matrix B. +*> On exit, the solution matrix X. +*> \endverbatim +*> +*> \param[in] LDB +*> \verbatim +*> LDB is INTEGER +*> The leading dimension of the array B. LDB >= max(1,N). +*> \endverbatim +*> +*> \param[out] INFO +*> \verbatim +*> INFO is INTEGER +*> = 0: successful exit +*> < 0: if INFO = -i, the i-th argument had an illegal value +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2011 +* +*> \ingroup complexSYcomputational +* +*> \par Contributors: +* ================== +*> +*> \verbatim +*> +*> November 2011, Igor Kozachenko, +*> Computer Science Division, +*> University of California, Berkeley +*> +*> September 2007, Sven Hammarling, Nicholas J. Higham, Craig Lucas, +*> School of Mathematics, +*> University of Manchester +*> +*> \endverbatim +* +* ===================================================================== + SUBROUTINE CSYTRS_ROOK( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, + $ INFO ) +* +* -- LAPACK computational routine (version 3.4.0) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2011 +* +* .. Scalar Arguments .. + CHARACTER UPLO + INTEGER INFO, LDA, LDB, N, NRHS +* .. +* .. Array Arguments .. + INTEGER IPIV( * ) + COMPLEX A( LDA, * ), B( LDB, * ) +* .. +* +* ===================================================================== +* +* .. Parameters .. + COMPLEX CONE + PARAMETER ( CONE = ( 1.0E+0, 0.0E+0 ) ) +* .. +* .. Local Scalars .. + LOGICAL UPPER + INTEGER J, K, KP + COMPLEX AK, AKM1, AKM1K, BK, BKM1, DENOM +* .. +* .. External Functions .. + LOGICAL LSAME + EXTERNAL LSAME +* .. +* .. External Subroutines .. + EXTERNAL CGEMV, CGERU, CSCAL, CSWAP, XERBLA +* .. +* .. Intrinsic Functions .. + INTRINSIC MAX +* .. +* .. Executable Statements .. +* + INFO = 0 + UPPER = LSAME( UPLO, 'U' ) + IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN + INFO = -1 + ELSE IF( N.LT.0 ) THEN + INFO = -2 + ELSE IF( NRHS.LT.0 ) THEN + INFO = -3 + ELSE IF( LDA.LT.MAX( 1, N ) ) THEN + INFO = -5 + ELSE IF( LDB.LT.MAX( 1, N ) ) THEN + INFO = -8 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'CSYTRS_ROOK', -INFO ) + RETURN + END IF +* +* Quick return if possible +* + IF( N.EQ.0 .OR. NRHS.EQ.0 ) + $ RETURN +* + IF( UPPER ) THEN +* +* Solve A*X = B, where A = U*D*U**T. +* +* First solve U*D*X = B, overwriting B with X. +* +* K is the main loop index, decreasing from N to 1 in steps of +* 1 or 2, depending on the size of the diagonal blocks. +* + K = N + 10 CONTINUE +* +* If K < 1, exit from loop. +* + IF( K.LT.1 ) + $ GO TO 30 +* + IF( IPIV( K ).GT.0 ) THEN +* +* 1 x 1 diagonal block +* +* Interchange rows K and IPIV(K). +* + KP = IPIV( K ) + IF( KP.NE.K ) + $ CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) +* +* Multiply by inv(U(K)), where U(K) is the transformation +* stored in column K of A. +* + CALL CGERU( K-1, NRHS, -CONE, A( 1, K ), 1, B( K, 1 ), LDB, + $ B( 1, 1 ), LDB ) +* +* Multiply by the inverse of the diagonal block. +* + CALL CSCAL( NRHS, CONE / A( K, K ), B( K, 1 ), LDB ) + K = K - 1 + ELSE +* +* 2 x 2 diagonal block +* +* Interchange rows K and -IPIV(K) THEN K-1 and -IPIV(K-1) +* + KP = -IPIV( K ) + IF( KP.NE.K ) + $ CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) +* + KP = -IPIV( K-1 ) + IF( KP.NE.K-1 ) + $ CALL CSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ), LDB ) +* +* Multiply by inv(U(K)), where U(K) is the transformation +* stored in columns K-1 and K of A. +* + IF( K.GT.2 ) THEN + CALL CGERU( K-2, NRHS,-CONE, A( 1, K ), 1, B( K, 1 ), + $ LDB, B( 1, 1 ), LDB ) + CALL CGERU( K-2, NRHS,-CONE, A( 1, K-1 ), 1, B( K-1, 1 ), + $ LDB, B( 1, 1 ), LDB ) + END IF +* +* Multiply by the inverse of the diagonal block. +* + AKM1K = A( K-1, K ) + AKM1 = A( K-1, K-1 ) / AKM1K + AK = A( K, K ) / AKM1K + DENOM = AKM1*AK - CONE + DO 20 J = 1, NRHS + BKM1 = B( K-1, J ) / AKM1K + BK = B( K, J ) / AKM1K + B( K-1, J ) = ( AK*BKM1-BK ) / DENOM + B( K, J ) = ( AKM1*BK-BKM1 ) / DENOM + 20 CONTINUE + K = K - 2 + END IF +* + GO TO 10 + 30 CONTINUE +* +* Next solve U**T *X = B, overwriting B with X. +* +* K is the main loop index, increasing from 1 to N in steps of +* 1 or 2, depending on the size of the diagonal blocks. +* + K = 1 + 40 CONTINUE +* +* If K > N, exit from loop. +* + IF( K.GT.N ) + $ GO TO 50 +* + IF( IPIV( K ).GT.0 ) THEN +* +* 1 x 1 diagonal block +* +* Multiply by inv(U**T(K)), where U(K) is the transformation +* stored in column K of A. +* + IF( K.GT.1 ) + $ CALL CGEMV( 'Transpose', K-1, NRHS, -CONE, B, + $ LDB, A( 1, K ), 1, CONE, B( K, 1 ), LDB ) +* +* Interchange rows K and IPIV(K). +* + KP = IPIV( K ) + IF( KP.NE.K ) + $ CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) + K = K + 1 + ELSE +* +* 2 x 2 diagonal block +* +* Multiply by inv(U**T(K+1)), where U(K+1) is the transformation +* stored in columns K and K+1 of A. +* + IF( K.GT.1 ) THEN + CALL CGEMV( 'Transpose', K-1, NRHS, -CONE, B, + $ LDB, A( 1, K ), 1, CONE, B( K, 1 ), LDB ) + CALL CGEMV( 'Transpose', K-1, NRHS, -CONE, B, + $ LDB, A( 1, K+1 ), 1, CONE, B( K+1, 1 ), LDB ) + END IF +* +* Interchange rows K and -IPIV(K) THEN K+1 and -IPIV(K+1). +* + KP = -IPIV( K ) + IF( KP.NE.K ) + $ CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) +* + KP = -IPIV( K+1 ) + IF( KP.NE.K+1 ) + $ CALL CSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ), LDB ) +* + K = K + 2 + END IF +* + GO TO 40 + 50 CONTINUE +* + ELSE +* +* Solve A*X = B, where A = L*D*L**T. +* +* First solve L*D*X = B, overwriting B with X. +* +* K is the main loop index, increasing from 1 to N in steps of +* 1 or 2, depending on the size of the diagonal blocks. +* + K = 1 + 60 CONTINUE +* +* If K > N, exit from loop. +* + IF( K.GT.N ) + $ GO TO 80 +* + IF( IPIV( K ).GT.0 ) THEN +* +* 1 x 1 diagonal block +* +* Interchange rows K and IPIV(K). +* + KP = IPIV( K ) + IF( KP.NE.K ) + $ CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) +* +* Multiply by inv(L(K)), where L(K) is the transformation +* stored in column K of A. +* + IF( K.LT.N ) + $ CALL CGERU( N-K, NRHS, -CONE, A( K+1, K ), 1, B( K, 1 ), + $ LDB, B( K+1, 1 ), LDB ) +* +* Multiply by the inverse of the diagonal block. +* + CALL CSCAL( NRHS, CONE / A( K, K ), B( K, 1 ), LDB ) + K = K + 1 + ELSE +* +* 2 x 2 diagonal block +* +* Interchange rows K and -IPIV(K) THEN K+1 and -IPIV(K+1) +* + KP = -IPIV( K ) + IF( KP.NE.K ) + $ CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) +* + KP = -IPIV( K+1 ) + IF( KP.NE.K+1 ) + $ CALL CSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ), LDB ) +* +* Multiply by inv(L(K)), where L(K) is the transformation +* stored in columns K and K+1 of A. +* + IF( K.LT.N-1 ) THEN + CALL CGERU( N-K-1, NRHS,-CONE, A( K+2, K ), 1, B( K, 1 ), + $ LDB, B( K+2, 1 ), LDB ) + CALL CGERU( N-K-1, NRHS,-CONE, A( K+2, K+1 ), 1, + $ B( K+1, 1 ), LDB, B( K+2, 1 ), LDB ) + END IF +* +* Multiply by the inverse of the diagonal block. +* + AKM1K = A( K+1, K ) + AKM1 = A( K, K ) / AKM1K + AK = A( K+1, K+1 ) / AKM1K + DENOM = AKM1*AK - CONE + DO 70 J = 1, NRHS + BKM1 = B( K, J ) / AKM1K + BK = B( K+1, J ) / AKM1K + B( K, J ) = ( AK*BKM1-BK ) / DENOM + B( K+1, J ) = ( AKM1*BK-BKM1 ) / DENOM + 70 CONTINUE + K = K + 2 + END IF +* + GO TO 60 + 80 CONTINUE +* +* Next solve L**T *X = B, overwriting B with X. +* +* K is the main loop index, decreasing from N to 1 in steps of +* 1 or 2, depending on the size of the diagonal blocks. +* + K = N + 90 CONTINUE +* +* If K < 1, exit from loop. +* + IF( K.LT.1 ) + $ GO TO 100 +* + IF( IPIV( K ).GT.0 ) THEN +* +* 1 x 1 diagonal block +* +* Multiply by inv(L**T(K)), where L(K) is the transformation +* stored in column K of A. +* + IF( K.LT.N ) + $ CALL CGEMV( 'Transpose', N-K, NRHS, -CONE, B( K+1, 1 ), + $ LDB, A( K+1, K ), 1, CONE, B( K, 1 ), LDB ) +* +* Interchange rows K and IPIV(K). +* + KP = IPIV( K ) + IF( KP.NE.K ) + $ CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) + K = K - 1 + ELSE +* +* 2 x 2 diagonal block +* +* Multiply by inv(L**T(K-1)), where L(K-1) is the transformation +* stored in columns K-1 and K of A. +* + IF( K.LT.N ) THEN + CALL CGEMV( 'Transpose', N-K, NRHS, -CONE, B( K+1, 1 ), + $ LDB, A( K+1, K ), 1, CONE, B( K, 1 ), LDB ) + CALL CGEMV( 'Transpose', N-K, NRHS, -CONE, B( K+1, 1 ), + $ LDB, A( K+1, K-1 ), 1, CONE, B( K-1, 1 ), + $ LDB ) + END IF +* +* Interchange rows K and -IPIV(K) THEN K-1 and -IPIV(K-1) +* + KP = -IPIV( K ) + IF( KP.NE.K ) + $ CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) +* + KP = -IPIV( K-1 ) + IF( KP.NE.K-1 ) + $ CALL CSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ), LDB ) +* + K = K - 2 + END IF +* + GO TO 90 + 100 CONTINUE + END IF +* + RETURN +* +* End of CSYTRS_ROOK +* + END |