--- /dev/null
+*> \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
--- /dev/null
+*> \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
--- /dev/null
+*> \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
--- /dev/null
+*> \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
--- /dev/null
+*> \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
--- /dev/null
+*> \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
--- /dev/null
+*> \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
projects / platform / upstream / lapack.git / commitdiff