--- /dev/null
+*> \brief \b CCHKHE_ROOK
+*
+* =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+* http://www.netlib.org/lapack/explore-html/
+*
+* Definition:
+* ===========
+*
+* SUBROUTINE CCHKHE_ROOK( DOTYPE, NN, NVAL, NNB, NBVAL, NNS, NSVAL,
+* THRESH, TSTERR, NMAX, A, AFAC, AINV, B, X,
+* XACT, WORK, RWORK, IWORK, NOUT )
+*
+* .. Scalar Arguments ..
+* LOGICAL TSTERR
+* INTEGER NMAX, NN, NNB, NNS, NOUT
+* REAL THRESH
+* ..
+* .. Array Arguments ..
+* LOGICAL DOTYPE( * )
+* INTEGER IWORK( * ), NBVAL( * ), NSVAL( * ), NVAL( * )
+* REAL RWORK( * )
+* COMPLEX A( * ), AFAC( * ), AINV( * ), B( * ),
+* $ WORK( * ), X( * ), XACT( * )
+* ..
+*
+*
+*> \par Purpose:
+* =============
+*>
+*> \verbatim
+*>
+*> CCHKHE_ROOK tests CHETRF_ROOK, -TRI_ROOK, -TRS_ROOK,
+*> and -CON_ROOK.
+*> \endverbatim
+*
+* Arguments:
+* ==========
+*
+*> \param[in] DOTYPE
+*> \verbatim
+*> DOTYPE is LOGICAL array, dimension (NTYPES)
+*> The matrix types to be used for testing. Matrices of type j
+*> (for 1 <= j <= NTYPES) are used for testing if DOTYPE(j) =
+*> .TRUE.; if DOTYPE(j) = .FALSE., then type j is not used.
+*> \endverbatim
+*>
+*> \param[in] NN
+*> \verbatim
+*> NN is INTEGER
+*> The number of values of N contained in the vector NVAL.
+*> \endverbatim
+*>
+*> \param[in] NVAL
+*> \verbatim
+*> NVAL is INTEGER array, dimension (NN)
+*> The values of the matrix dimension N.
+*> \endverbatim
+*>
+*> \param[in] NNB
+*> \verbatim
+*> NNB is INTEGER
+*> The number of values of NB contained in the vector NBVAL.
+*> \endverbatim
+*>
+*> \param[in] NBVAL
+*> \verbatim
+*> NBVAL is INTEGER array, dimension (NBVAL)
+*> The values of the blocksize NB.
+*> \endverbatim
+*>
+*> \param[in] NNS
+*> \verbatim
+*> NNS is INTEGER
+*> The number of values of NRHS contained in the vector NSVAL.
+*> \endverbatim
+*>
+*> \param[in] NSVAL
+*> \verbatim
+*> NSVAL is INTEGER array, dimension (NNS)
+*> The values of the number of right hand sides NRHS.
+*> \endverbatim
+*>
+*> \param[in] THRESH
+*> \verbatim
+*> THRESH is REAL
+*> The threshold value for the test ratios. A result is
+*> included in the output file if RESULT >= THRESH. To have
+*> every test ratio printed, use THRESH = 0.
+*> \endverbatim
+*>
+*> \param[in] TSTERR
+*> \verbatim
+*> TSTERR is LOGICAL
+*> Flag that indicates whether error exits are to be tested.
+*> \endverbatim
+*>
+*> \param[in] NMAX
+*> \verbatim
+*> NMAX is INTEGER
+*> The maximum value permitted for N, used in dimensioning the
+*> work arrays.
+*> \endverbatim
+*>
+*> \param[out] A
+*> \verbatim
+*> A is COMPLEX array, dimension (NMAX*NMAX)
+*> \endverbatim
+*>
+*> \param[out] AFAC
+*> \verbatim
+*> AFAC is COMPLEX array, dimension (NMAX*NMAX)
+*> \endverbatim
+*>
+*> \param[out] AINV
+*> \verbatim
+*> AINV is COMPLEX array, dimension (NMAX*NMAX)
+*> \endverbatim
+*>
+*> \param[out] B
+*> \verbatim
+*> B is COMPLEX array, dimension (NMAX*NSMAX)
+*> where NSMAX is the largest entry in NSVAL.
+*> \endverbatim
+*>
+*> \param[out] X
+*> \verbatim
+*> X is COMPLEX array, dimension (NMAX*NSMAX)
+*> \endverbatim
+*>
+*> \param[out] XACT
+*> \verbatim
+*> XACT is COMPLEX array, dimension (NMAX*NSMAX)
+*> \endverbatim
+*>
+*> \param[out] WORK
+*> \verbatim
+*> WORK is COMPLEX array, dimension
+*> (NMAX*max(3,NSMAX))
+*> \endverbatim
+*>
+*> \param[out] RWORK
+*> \verbatim
+*> RWORK is REAL array, dimension
+*> (max(NMAX,2*NSMAX))
+*> \endverbatim
+*>
+*> \param[out] IWORK
+*> \verbatim
+*> IWORK is INTEGER array, dimension (2*NMAX)
+*> \endverbatim
+*>
+*> \param[in] NOUT
+*> \verbatim
+*> NOUT is INTEGER
+*> The unit number for output.
+*> \endverbatim
+*
+* Authors:
+* ========
+*
+*> \author Univ. of Tennessee
+*> \author Univ. of California Berkeley
+*> \author Univ. of Colorado Denver
+*> \author NAG Ltd.
+*
+*> \date April 2013
+*
+*> \ingroup complex_lin
+*
+* =====================================================================
+ SUBROUTINE CCHKHE_ROOK( DOTYPE, NN, NVAL, NNB, NBVAL, NNS, NSVAL,
+ $ THRESH, TSTERR, NMAX, A, AFAC, AINV, B, X,
+ $ XACT, WORK, RWORK, IWORK, NOUT )
+*
+* -- LAPACK test routine (version 3.4.1) --
+* -- LAPACK is a software package provided by Univ. of Tennessee, --
+* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+* April 2012
+*
+* .. Scalar Arguments ..
+ LOGICAL TSTERR
+ INTEGER NMAX, NN, NNB, NNS, NOUT
+ REAL THRESH
+* ..
+* .. Array Arguments ..
+ LOGICAL DOTYPE( * )
+ INTEGER IWORK( * ), NBVAL( * ), NSVAL( * ), NVAL( * )
+ REAL RWORK( * )
+ COMPLEX A( * ), AFAC( * ), AINV( * ), B( * ),
+ $ WORK( * ), X( * ), XACT( * )
+* ..
+*
+* =====================================================================
+*
+* .. Parameters ..
+ REAL ZERO, ONE
+ PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 )
+ REAL ONEHALF
+ PARAMETER ( ONEHALF = 0.5E+0 )
+ REAL EIGHT, SEVTEN
+ PARAMETER ( EIGHT = 8.0E+0, SEVTEN = 17.0E+0 )
+ COMPLEX CZERO
+ PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ) )
+ INTEGER NTYPES
+ PARAMETER ( NTYPES = 10 )
+ INTEGER NTESTS
+ PARAMETER ( NTESTS = 7 )
+* ..
+* .. Local Scalars ..
+ LOGICAL TRFCON, ZEROT
+ CHARACTER DIST, TYPE, UPLO, XTYPE
+ CHARACTER*3 PATH, MATPATH
+ INTEGER I, I1, I2, IMAT, IN, INB, INFO, IOFF, IRHS,
+ $ ITEMP, ITEMP2, IUPLO, IZERO, J, K, KL, KU, LDA,
+ $ LWORK, MODE, N, NB, NERRS, NFAIL, NIMAT, NRHS,
+ $ NRUN, NT
+ REAL ALPHA, ANORM, CNDNUM, CONST, LAM_MAX, LAM_MIN,
+ $ RCOND, RCONDC, STEMP
+* ..
+* .. Local Arrays ..
+ CHARACTER UPLOS( 2 )
+ INTEGER ISEED( 4 ), ISEEDY( 4 )
+ REAL RESULT( NTESTS )
+ COMPLEX BLOCK( 2, 2 ), CDUMMY( 1 )
+* ..
+* .. External Functions ..
+ REAL CLANGE, CLANHE, SGET06
+ EXTERNAL CLANGE, CLANHE, SGET06
+* ..
+* .. External Subroutines ..
+ EXTERNAL ALAERH, ALAHD, ALASUM, CERRHE, CHEEVX, CGET04,
+ $ CLACPY, CLARHS, CLATB4, CLATMS, CPOT02,
+ $ CPOT03, CHECON_ROOK, CHET01_ROOK, CHETRF_ROOK,
+ $ CHETRI_ROOK, CHETRS_ROOK, XLAENV
+* ..
+* .. Intrinsic Functions ..
+ INTRINSIC ABS, MAX, MIN, SQRT
+* ..
+* .. Scalars in Common ..
+ LOGICAL LERR, OK
+ CHARACTER*32 SRNAMT
+ INTEGER INFOT, NUNIT
+* ..
+* .. Common blocks ..
+ COMMON / INFOC / INFOT, NUNIT, OK, LERR
+ COMMON / SRNAMC / SRNAMT
+* ..
+* .. Data statements ..
+ DATA ISEEDY / 1988, 1989, 1990, 1991 /
+ DATA UPLOS / 'U', 'L' /
+* ..
+* .. Executable Statements ..
+*
+* Initialize constants and the random number seed.
+*
+ ALPHA = ( ONE+SQRT( SEVTEN ) ) / EIGHT
+*
+* Test path
+*
+ PATH( 1: 1 ) = 'Complex precision'
+ PATH( 2: 3 ) = 'HR'
+*
+* Path to generate matrices
+*
+ MATPATH( 1: 1 ) = 'Complex precision'
+ MATPATH( 2: 3 ) = 'HE'
+*
+ NRUN = 0
+ NFAIL = 0
+ NERRS = 0
+ DO 10 I = 1, 4
+ ISEED( I ) = ISEEDY( I )
+ 10 CONTINUE
+*
+* Test the error exits
+*
+ IF( TSTERR )
+ $ CALL CERRHE( PATH, NOUT )
+ INFOT = 0
+*
+* Set the minimum block size for which the block routine should
+* be used, which will be later returned by ILAENV
+*
+ CALL XLAENV( 2, 2 )
+*
+* Do for each value of N in NVAL
+*
+ DO 270 IN = 1, NN
+ N = NVAL( IN )
+ LDA = MAX( N, 1 )
+ XTYPE = 'N'
+ NIMAT = NTYPES
+ IF( N.LE.0 )
+ $ NIMAT = 1
+*
+ IZERO = 0
+*
+* Do for each value of matrix type IMAT
+*
+ DO 260 IMAT = 1, NIMAT
+*
+* Do the tests only if DOTYPE( IMAT ) is true.
+*
+ IF( .NOT.DOTYPE( IMAT ) )
+ $ GO TO 260
+*
+* Skip types 3, 4, 5, or 6 if the matrix size is too small.
+*
+ ZEROT = IMAT.GE.3 .AND. IMAT.LE.6
+ IF( ZEROT .AND. N.LT.IMAT-2 )
+ $ GO TO 260
+*
+* Do first for UPLO = 'U', then for UPLO = 'L'
+*
+ DO 250 IUPLO = 1, 2
+ UPLO = UPLOS( IUPLO )
+*
+* Begin generate the test matrix A.
+*
+* Set up parameters with CLATB4 for the matrix generator
+* based on the type of matrix to be generated.
+*
+ CALL CLATB4( MATPATH, IMAT, N, N, TYPE, KL, KU, ANORM,
+ $ MODE, CNDNUM, DIST )
+*
+* Generate a matrix with CLATMS.
+*
+ SRNAMT = 'CLATMS'
+ CALL CLATMS( N, N, DIST, ISEED, TYPE, RWORK, MODE,
+ $ CNDNUM, ANORM, KL, KU, UPLO, A, LDA,
+ $ WORK, INFO )
+*
+* Check error code from CLATMS and handle error.
+*
+ IF( INFO.NE.0 ) THEN
+ CALL ALAERH( PATH, 'CLATMS', INFO, 0, UPLO, N, N,
+ $ -1, -1, -1, IMAT, NFAIL, NERRS, NOUT )
+*
+* Skip all tests for this generated matrix
+*
+ GO TO 250
+ END IF
+*
+* For matrix types 3-6, zero one or more rows and
+* columns of the matrix to test that INFO is returned
+* correctly.
+*
+ IF( ZEROT ) THEN
+ IF( IMAT.EQ.3 ) THEN
+ IZERO = 1
+ ELSE IF( IMAT.EQ.4 ) THEN
+ IZERO = N
+ ELSE
+ IZERO = N / 2 + 1
+ END IF
+*
+ IF( IMAT.LT.6 ) THEN
+*
+* Set row and column IZERO to zero.
+*
+ IF( IUPLO.EQ.1 ) THEN
+ IOFF = ( IZERO-1 )*LDA
+ DO 20 I = 1, IZERO - 1
+ A( IOFF+I ) = CZERO
+ 20 CONTINUE
+ IOFF = IOFF + IZERO
+ DO 30 I = IZERO, N
+ A( IOFF ) = CZERO
+ IOFF = IOFF + LDA
+ 30 CONTINUE
+ ELSE
+ IOFF = IZERO
+ DO 40 I = 1, IZERO - 1
+ A( IOFF ) = CZERO
+ IOFF = IOFF + LDA
+ 40 CONTINUE
+ IOFF = IOFF - IZERO
+ DO 50 I = IZERO, N
+ A( IOFF+I ) = CZERO
+ 50 CONTINUE
+ END IF
+ ELSE
+ IOFF = 0
+ IF( IUPLO.EQ.1 ) THEN
+*
+* Set the first IZERO rows and columns to zero.
+*
+ DO 70 J = 1, N
+ I2 = MIN( J, IZERO )
+ DO 60 I = 1, I2
+ A( IOFF+I ) = CZERO
+ 60 CONTINUE
+ IOFF = IOFF + LDA
+ 70 CONTINUE
+ ELSE
+*
+* Set the last IZERO rows and columns to zero.
+*
+ DO 90 J = 1, N
+ I1 = MAX( J, IZERO )
+ DO 80 I = I1, N
+ A( IOFF+I ) = CZERO
+ 80 CONTINUE
+ IOFF = IOFF + LDA
+ 90 CONTINUE
+ END IF
+ END IF
+ ELSE
+ IZERO = 0
+ END IF
+*
+* End generate the test matrix A.
+*
+* Do for each value of NB in NBVAL
+*
+ DO 240 INB = 1, NNB
+*
+* Set the optimal blocksize, which will be later
+* returned by ILAENV.
+*
+ NB = NBVAL( INB )
+ CALL XLAENV( 1, NB )
+*
+* Copy the test matrix A into matrix AFAC which
+* will be factorized in place. This is needed to
+* preserve the test matrix A for subsequent tests.
+*
+ CALL CLACPY( UPLO, N, N, A, LDA, AFAC, LDA )
+*
+* Compute the L*D*L**T or U*D*U**T factorization of the
+* matrix. IWORK stores details of the interchanges and
+* the block structure of D. AINV is a work array for
+* block factorization, LWORK is the length of AINV.
+*
+ LWORK = MAX( 2, NB )*LDA
+ SRNAMT = 'CHETRF_ROOK'
+ CALL CHETRF_ROOK( UPLO, N, AFAC, LDA, IWORK, AINV,
+ $ LWORK, INFO )
+*
+* Adjust the expected value of INFO to account for
+* pivoting.
+*
+ K = IZERO
+ IF( K.GT.0 ) THEN
+ 100 CONTINUE
+ IF( IWORK( K ).LT.0 ) THEN
+ IF( IWORK( K ).NE.-K ) THEN
+ K = -IWORK( K )
+ GO TO 100
+ END IF
+ ELSE IF( IWORK( K ).NE.K ) THEN
+ K = IWORK( K )
+ GO TO 100
+ END IF
+ END IF
+*
+* Check error code from CHETRF_ROOK and handle error.
+*
+ IF( INFO.NE.K)
+ $ CALL ALAERH( PATH, 'CHETRF_ROOK', INFO, K,
+ $ UPLO, N, N, -1, -1, NB, IMAT,
+ $ NFAIL, NERRS, NOUT )
+*
+* Set the condition estimate flag if the INFO is not 0.
+*
+ IF( INFO.NE.0 ) THEN
+ TRFCON = .TRUE.
+ ELSE
+ TRFCON = .FALSE.
+ END IF
+*
+*+ TEST 1
+* Reconstruct matrix from factors and compute residual.
+*
+ CALL CHET01_ROOK( UPLO, N, A, LDA, AFAC, LDA, IWORK,
+ $ AINV, LDA, RWORK, RESULT( 1 ) )
+ NT = 1
+*
+*+ TEST 2
+* Form the inverse and compute the residual,
+* if the factorization was competed without INFO > 0
+* (i.e. there is no zero rows and columns).
+* Do it only for the first block size.
+*
+ IF( INB.EQ.1 .AND. .NOT.TRFCON ) THEN
+ CALL CLACPY( UPLO, N, N, AFAC, LDA, AINV, LDA )
+ SRNAMT = 'CHETRI_ROOK'
+ CALL CHETRI_ROOK( UPLO, N, AINV, LDA, IWORK, WORK,
+ $ INFO )
+*
+* Check error code from CHETRI_ROOK and handle error.
+*
+ IF( INFO.NE.0 )
+ $ CALL ALAERH( PATH, 'CHETRI_ROOK', INFO, -1,
+ $ UPLO, N, N, -1, -1, -1, IMAT,
+ $ NFAIL, NERRS, NOUT )
+*
+* Compute the residual for a Hermitian matrix times
+* its inverse.
+*
+ CALL CPOT03( UPLO, N, A, LDA, AINV, LDA, WORK, LDA,
+ $ RWORK, RCONDC, RESULT( 2 ) )
+ NT = 2
+ END IF
+*
+* Print information about the tests that did not pass
+* the threshold.
+*
+ DO 110 K = 1, NT
+ IF( RESULT( K ).GE.THRESH ) THEN
+ IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
+ $ CALL ALAHD( NOUT, PATH )
+ WRITE( NOUT, FMT = 9999 )UPLO, N, NB, IMAT, K,
+ $ RESULT( K )
+ NFAIL = NFAIL + 1
+ END IF
+ 110 CONTINUE
+ NRUN = NRUN + NT
+*
+*+ TEST 3
+* Compute largest element in U or L
+*
+ RESULT( 3 ) = ZERO
+ STEMP = ZERO
+*
+ CONST = ( ( ALPHA**2-ONE ) / ( ALPHA**2-ONEHALF ) ) /
+ $ ( ONE-ALPHA )
+*
+ IF( IUPLO.EQ.1 ) THEN
+*
+* Compute largest element in U
+*
+ K = N
+ 120 CONTINUE
+ IF( K.LE.1 )
+ $ GO TO 130
+*
+ IF( IWORK( K ).GT.ZERO ) THEN
+*
+* Get max absolute value from elements
+* in column k in in U
+*
+ STEMP = CLANGE( 'M', K-1, 1,
+ $ AFAC( ( K-1 )*LDA+1 ), LDA, RWORK )
+ ELSE
+*
+* Get max absolute value from elements
+* in columns k and k-1 in U
+*
+ STEMP = CLANGE( 'M', K-2, 2,
+ $ AFAC( ( K-2 )*LDA+1 ), LDA, RWORK )
+ K = K - 1
+*
+ END IF
+*
+* STEMP should be bounded by CONST
+*
+ STEMP = STEMP - CONST + THRESH
+ IF( STEMP.GT.RESULT( 3 ) )
+ $ RESULT( 3 ) = STEMP
+*
+ K = K - 1
+*
+ GO TO 120
+ 130 CONTINUE
+*
+ ELSE
+*
+* Compute largest element in L
+*
+ K = 1
+ 140 CONTINUE
+ IF( K.GE.N )
+ $ GO TO 150
+*
+ IF( IWORK( K ).GT.ZERO ) THEN
+*
+* Get max absolute value from elements
+* in column k in L
+*
+ STEMP = CLANGE( 'M', N-K, 1,
+ $ AFAC( ( K-1 )*LDA+K+1 ), LDA, RWORK )
+ ELSE
+*
+* Get max absolute value from elements
+* in columns k and k+1 in L
+*
+ STEMP = CLANGE( 'M', N-K-1, 2,
+ $ AFAC( ( K-1 )*LDA+K+2 ), LDA, RWORK )
+ K = K + 1
+*
+ END IF
+*
+* STEMP should be bounded by CONST
+*
+ STEMP = STEMP - CONST + THRESH
+ IF( STEMP.GT.RESULT( 3 ) )
+ $ RESULT( 3 ) = STEMP
+*
+ K = K + 1
+*
+ GO TO 140
+ 150 CONTINUE
+ END IF
+*
+*
+*+ TEST 4
+* Compute largest 2-Norm of 2-by-2 diag blocks
+*
+ RESULT( 4 ) = ZERO
+ STEMP = ZERO
+*
+ CONST = ( ( ALPHA**2-ONE ) / ( ALPHA**2-ONEHALF ) )*
+ $ ( ( ONE + ALPHA ) / ( ONE - ALPHA ) )
+*
+ IF( IUPLO.EQ.1 ) THEN
+*
+* Loop backward for UPLO = 'U'
+*
+ K = N
+ 160 CONTINUE
+ IF( K.LE.1 )
+ $ GO TO 170
+*
+ IF( IWORK( K ).LT.ZERO ) THEN
+*
+* Get the two eigenvalues of a 2-by-2 block,
+* store them in WORK array
+*
+ BLOCK( 1, 1 ) = AFAC( ( K-2 )*LDA+K-1 )
+ BLOCK( 2, 1 ) = AFAC( ( K-2 )*LDA+K )
+ BLOCK( 1, 2 ) = BLOCK( 2, 1 )
+ BLOCK( 2, 2 ) = AFAC( (K-1)*LDA+K )
+*
+ CALL CHEEVX( 'N', 'N', 'N', 'N', 2, BLOCK,
+ $ 2, WORK, CDUMMY, 1, CDUMMY, 1,
+ $ ITEMP, ITEMP2, RWORK, STEMP,
+ $ RWORK( 3 ), RWORK( 5 ), WORK( 3 ),
+ $ 4, RWORK( 7 ), INFO )
+*
+ LAM_MAX = MAX( ABS( WORK( 1 ) ),
+ $ ABS( WORK( 2 ) ) )
+ LAM_MIN = MIN( ABS( WORK( 1 ) ),
+ $ ABS( WORK( 2 ) ) )
+*
+ STEMP = LAM_MAX / LAM_MIN
+*
+* STEMP should be bounded by CONST
+*
+ STEMP = ABS( STEMP ) - CONST + THRESH
+ IF( STEMP.GT.RESULT( 4 ) )
+ $ RESULT( 4 ) = STEMP
+ K = K - 1
+*
+ END IF
+*
+ K = K - 1
+*
+ GO TO 160
+ 170 CONTINUE
+*
+ ELSE
+*
+* Loop forward for UPLO = 'L'
+*
+ K = 1
+ 180 CONTINUE
+ IF( K.GE.N )
+ $ GO TO 190
+*
+ IF( IWORK( K ).LT.ZERO ) THEN
+*
+* Get the two eigenvalues of a 2-by-2 block,
+* store them in WORK array
+*
+ BLOCK( 1, 1 ) = AFAC( ( K-1 )*LDA+K )
+ BLOCK( 2, 1 ) = AFAC( ( K-1 )*LDA+K+1 )
+ BLOCK( 1, 2 ) = BLOCK( 2, 1 )
+ BLOCK( 2, 2 ) = AFAC( K*LDA+K+1 )
+*
+ CALL CHEEVX( 'N', 'N', 'N', 'N', 2, BLOCK,
+ $ 2, WORK, CDUMMY, 1, CDUMMY, 1,
+ $ ITEMP, ITEMP2, RWORK, STEMP,
+ $ RWORK( 3 ), RWORK( 5 ), WORK( 3 ),
+ $ 4, RWORK( 7 ), INFO )
+*
+ LAM_MAX = MAX( ABS( WORK( 1 ) ),
+ $ ABS( WORK( 2 ) ) )
+ LAM_MIN = MIN( ABS( WORK( 1 ) ),
+ $ ABS( WORK( 2 ) ) )
+*
+ STEMP = LAM_MAX / LAM_MIN
+*
+* STEMP should be bounded by CONST
+*
+ STEMP = ABS( STEMP ) - CONST + THRESH
+ IF( STEMP.GT.RESULT( 4 ) )
+ $ RESULT( 4 ) = STEMP
+ K = K + 1
+*
+ END IF
+*
+ K = K + 1
+*
+ GO TO 180
+ 190 CONTINUE
+ END IF
+*
+* Print information about the tests that did not pass
+* the threshold.
+*
+ DO 200 K = 3, 4
+ IF( RESULT( K ).GE.THRESH ) THEN
+ IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
+ $ CALL ALAHD( NOUT, PATH )
+ WRITE( NOUT, FMT = 9999 )UPLO, N, NB, IMAT, K,
+ $ RESULT( K )
+ NFAIL = NFAIL + 1
+ END IF
+ 200 CONTINUE
+ NRUN = NRUN + NT
+*
+* Skip the other tests if this is not the first block
+* size.
+*
+ IF( INB.GT.1 )
+ $ GO TO 240
+*
+* Do only the condition estimate if INFO is not 0.
+*
+ IF( TRFCON ) THEN
+ RCONDC = ZERO
+ GO TO 230
+ END IF
+*
+ DO 220 IRHS = 1, NNS
+ NRHS = NSVAL( IRHS )
+*
+* Begin loop over NRHS values
+*
+*
+*+ TEST 5 ( Using TRS_ROOK)
+* Solve and compute residual for A * X = B.
+*
+* Choose a set of NRHS random solution vectors
+* stored in XACT and set up the right hand side B
+*
+ SRNAMT = 'CLARHS'
+ CALL CLARHS( MATPATH, XTYPE, UPLO, ' ', N, N,
+ $ KL, KU, NRHS, A, LDA, XACT, LDA,
+ $ B, LDA, ISEED, INFO )
+ CALL CLACPY( 'Full', N, NRHS, B, LDA, X, LDA )
+*
+ SRNAMT = 'CHETRS_ROOK'
+ CALL CHETRS_ROOK( UPLO, N, NRHS, AFAC, LDA, IWORK,
+ $ X, LDA, INFO )
+*
+* Check error code from CHETRS_ROOK and handle error.
+*
+ IF( INFO.NE.0 )
+ $ CALL ALAERH( PATH, 'CHETRS_ROOK', INFO, 0,
+ $ UPLO, N, N, -1, -1, NRHS, IMAT,
+ $ NFAIL, NERRS, NOUT )
+*
+ CALL CLACPY( 'Full', N, NRHS, B, LDA, WORK, LDA )
+*
+* Compute the residual for the solution
+*
+ CALL CPOT02( UPLO, N, NRHS, A, LDA, X, LDA, WORK,
+ $ LDA, RWORK, RESULT( 5 ) )
+*
+*+ TEST 6
+* Check solution from generated exact solution.
+*
+ CALL CGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC,
+ $ RESULT( 6 ) )
+*
+* Print information about the tests that did not pass
+* the threshold.
+*
+ DO 210 K = 5, 6
+ IF( RESULT( K ).GE.THRESH ) THEN
+ IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
+ $ CALL ALAHD( NOUT, PATH )
+ WRITE( NOUT, FMT = 9998 )UPLO, N, NRHS,
+ $ IMAT, K, RESULT( K )
+ NFAIL = NFAIL + 1
+ END IF
+ 210 CONTINUE
+ NRUN = NRUN + 2
+*
+* End loop over NRHS values
+*
+ 220 CONTINUE
+*
+*+ TEST 7
+* Get an estimate of RCOND = 1/CNDNUM.
+*
+ 230 CONTINUE
+ ANORM = CLANHE( '1', UPLO, N, A, LDA, RWORK )
+ SRNAMT = 'CHECON_ROOK'
+ CALL CHECON_ROOK( UPLO, N, AFAC, LDA, IWORK, ANORM,
+ $ RCOND, WORK, INFO )
+*
+* Check error code from CHECON_ROOK and handle error.
+*
+ IF( INFO.NE.0 )
+ $ CALL ALAERH( PATH, 'CHECON_ROOK', INFO, 0,
+ $ UPLO, N, N, -1, -1, -1, IMAT,
+ $ NFAIL, NERRS, NOUT )
+*
+* Compute the test ratio to compare to values of RCOND
+*
+ RESULT( 7 ) = SGET06( RCOND, RCONDC )
+*
+* Print information about the tests that did not pass
+* the threshold.
+*
+ IF( RESULT( 7 ).GE.THRESH ) THEN
+ IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
+ $ CALL ALAHD( NOUT, PATH )
+ WRITE( NOUT, FMT = 9997 )UPLO, N, IMAT, 7,
+ $ RESULT( 7 )
+ NFAIL = NFAIL + 1
+ END IF
+ NRUN = NRUN + 1
+ 240 CONTINUE
+*
+ 250 CONTINUE
+ 260 CONTINUE
+ 270 CONTINUE
+*
+* Print a summary of the results.
+*
+ CALL ALASUM( PATH, NOUT, NFAIL, NRUN, NERRS )
+*
+ 9999 FORMAT( ' UPLO = ''', A1, ''', N =', I5, ', NB =', I4, ', type ',
+ $ I2, ', test ', I2, ', ratio =', G12.5 )
+ 9998 FORMAT( ' UPLO = ''', A1, ''', N =', I5, ', NRHS=', I3, ', type ',
+ $ I2, ', test(', I2, ') =', G12.5 )
+ 9997 FORMAT( ' UPLO = ''', A1, ''', N =', I5, ',', 10X, ' type ', I2,
+ $ ', test(', I2, ') =', G12.5 )
+ RETURN
+*
+* End of CCHKHE_ROOK
+*
+ END
--- /dev/null
+*> \brief \b CHET01_ROOK
+*
+* =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+* http://www.netlib.org/lapack/explore-html/
+*
+* Definition:
+* ===========
+*
+* SUBROUTINE CHET01_ROOK( UPLO, N, A, LDA, AFAC, LDAFAC, IPIV, C, LDC,
+* RWORK, RESID )
+*
+* .. Scalar Arguments ..
+* CHARACTER UPLO
+* INTEGER LDA, LDAFAC, LDC, N
+* REAL RESID
+* ..
+* .. Array Arguments ..
+* INTEGER IPIV( * )
+* REAL RWORK( * )
+* COMPLEX A( LDA, * ), AFAC( LDAFAC, * ), C( LDC, * )
+* ..
+*
+*
+*> \par Purpose:
+* =============
+*>
+*> \verbatim
+*>
+*> CHET01_ROOK reconstructs a complex Hermitian indefinite matrix A from its
+*> block L*D*L' or U*D*U' factorization and computes the residual
+*> norm( C - A ) / ( N * norm(A) * EPS ),
+*> where C is the reconstructed matrix, EPS is the machine epsilon,
+*> L' is the transpose of L, and U' is the transpose of U.
+*> \endverbatim
+*
+* Arguments:
+* ==========
+*
+*> \param[in] UPLO
+*> \verbatim
+*> UPLO is CHARACTER*1
+*> Specifies whether the upper or lower triangular part of the
+*> complex Hermitian matrix A is stored:
+*> = 'U': Upper triangular
+*> = 'L': Lower triangular
+*> \endverbatim
+*>
+*> \param[in] N
+*> \verbatim
+*> N is INTEGER
+*> The number of rows and columns of the matrix A. N >= 0.
+*> \endverbatim
+*>
+*> \param[in] A
+*> \verbatim
+*> A is COMPLEX array, dimension (LDA,N)
+*> The original complex Hermitian matrix A.
+*> \endverbatim
+*>
+*> \param[in] LDA
+*> \verbatim
+*> LDA is INTEGER
+*> The leading dimension of the array A. LDA >= max(1,N)
+*> \endverbatim
+*>
+*> \param[in] AFAC
+*> \verbatim
+*> AFAC is COMPLEX array, dimension (LDAFAC,N)
+*> The factored form of the matrix A. AFAC contains the block
+*> diagonal matrix D and the multipliers used to obtain the
+*> factor L or U from the block L*D*L' or U*D*U' factorization
+*> as computed by CSYTRF_ROOK.
+*> \endverbatim
+*>
+*> \param[in] LDAFAC
+*> \verbatim
+*> LDAFAC is INTEGER
+*> The leading dimension of the array AFAC. LDAFAC >= max(1,N).
+*> \endverbatim
+*>
+*> \param[in] IPIV
+*> \verbatim
+*> IPIV is INTEGER array, dimension (N)
+*> The pivot indices from CSYTRF_ROOK.
+*> \endverbatim
+*>
+*> \param[out] C
+*> \verbatim
+*> C is COMPLEX array, dimension (LDC,N)
+*> \endverbatim
+*>
+*> \param[in] LDC
+*> \verbatim
+*> LDC is INTEGER
+*> The leading dimension of the array C. LDC >= max(1,N).
+*> \endverbatim
+*>
+*> \param[out] RWORK
+*> \verbatim
+*> RWORK is REAL array, dimension (N)
+*> \endverbatim
+*>
+*> \param[out] RESID
+*> \verbatim
+*> RESID is REAL
+*> If UPLO = 'L', norm(L*D*L' - A) / ( N * norm(A) * EPS )
+*> If UPLO = 'U', norm(U*D*U' - A) / ( N * norm(A) * EPS )
+*> \endverbatim
+*
+* Authors:
+* ========
+*
+*> \author Univ. of Tennessee
+*> \author Univ. of California Berkeley
+*> \author Univ. of Colorado Denver
+*> \author NAG Ltd.
+*
+*> \date April 2013
+*
+*> \ingroup complex_lin
+*
+* =====================================================================
+ SUBROUTINE CHET01_ROOK( UPLO, N, A, LDA, AFAC, LDAFAC, IPIV, C,
+ $ LDC, RWORK, RESID )
+*
+* -- LAPACK test 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 LDA, LDAFAC, LDC, N
+ REAL RESID
+* ..
+* .. Array Arguments ..
+ INTEGER IPIV( * )
+ REAL RWORK( * )
+ COMPLEX A( LDA, * ), AFAC( LDAFAC, * ), C( LDC, * )
+* ..
+*
+* =====================================================================
+*
+* .. Parameters ..
+ REAL ZERO, ONE
+ PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 )
+ COMPLEX CZERO, CONE
+ PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ),
+ $ CONE = ( 1.0E+0, 0.0E+0 ) )
+* ..
+* .. Local Scalars ..
+ INTEGER I, INFO, J
+ REAL ANORM, EPS
+* ..
+* .. External Functions ..
+ LOGICAL LSAME
+ REAL CLANHE, SLAMCH
+ EXTERNAL LSAME, CLANHE, SLAMCH
+* ..
+* .. External Subroutines ..
+ EXTERNAL CLASET, CLAVHE_ROOK
+* ..
+* .. Intrinsic Functions ..
+ INTRINSIC AIMAG, REAL
+* ..
+* .. Executable Statements ..
+*
+* Quick exit if N = 0.
+*
+ IF( N.LE.0 ) THEN
+ RESID = ZERO
+ RETURN
+ END IF
+*
+* Determine EPS and the norm of A.
+*
+ EPS = SLAMCH( 'Epsilon' )
+ ANORM = CLANHE( '1', UPLO, N, A, LDA, RWORK )
+*
+* Check the imaginary parts of the diagonal elements and return with
+* an error code if any are nonzero.
+*
+ DO 10 J = 1, N
+ IF( AIMAG( AFAC( J, J ) ).NE.ZERO ) THEN
+ RESID = ONE / EPS
+ RETURN
+ END IF
+ 10 CONTINUE
+*
+* Initialize C to the identity matrix.
+*
+ CALL CLASET( 'Full', N, N, CZERO, CONE, C, LDC )
+*
+* Call CLAVHE_ROOK to form the product D * U' (or D * L' ).
+*
+ CALL CLAVHE_ROOK( UPLO, 'Conjugate', 'Non-unit', N, N, AFAC,
+ $ LDAFAC, IPIV, C, LDC, INFO )
+*
+* Call CLAVHE_ROOK again to multiply by U (or L ).
+*
+ CALL CLAVHE_ROOK( UPLO, 'No transpose', 'Unit', N, N, AFAC,
+ $ LDAFAC, IPIV, C, LDC, INFO )
+*
+* Compute the difference C - A .
+*
+ IF( LSAME( UPLO, 'U' ) ) THEN
+ DO 30 J = 1, N
+ DO 20 I = 1, J - 1
+ C( I, J ) = C( I, J ) - A( I, J )
+ 20 CONTINUE
+ C( J, J ) = C( J, J ) - REAL( A( J, J ) )
+ 30 CONTINUE
+ ELSE
+ DO 50 J = 1, N
+ C( J, J ) = C( J, J ) - REAL( A( J, J ) )
+ DO 40 I = J + 1, N
+ C( I, J ) = C( I, J ) - A( I, J )
+ 40 CONTINUE
+ 50 CONTINUE
+ END IF
+*
+* Compute norm( C - A ) / ( N * norm(A) * EPS )
+*
+ RESID = CLANHE( '1', UPLO, N, C, LDC, RWORK )
+*
+ IF( ANORM.LE.ZERO ) THEN
+ IF( RESID.NE.ZERO )
+ $ RESID = ONE / EPS
+ ELSE
+ RESID = ( ( RESID/REAL( N ) )/ANORM ) / EPS
+ END IF
+*
+ RETURN
+*
+* End of CHET01_ROOK
+*
+ END
--- /dev/null
+*> \brief \b CLAVHE_ROOK
+*
+* =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+* http://www.netlib.org/lapack/explore-html/
+*
+* Definition:
+* ===========
+*
+* SUBROUTINE CLAVHE_ROOK( UPLO, TRANS, DIAG, N, NRHS, A, LDA, IPIV, B,
+* LDB, INFO )
+*
+* .. Scalar Arguments ..
+* CHARACTER DIAG, TRANS, UPLO
+* INTEGER INFO, LDA, LDB, N, NRHS
+* ..
+* .. Array Arguments ..
+* INTEGER IPIV( * )
+* COMPLEX A( LDA, * ), B( LDB, * )
+* ..
+*
+*
+*> \par Purpose:
+* =============
+*>
+*> \verbatim
+*>
+*> CLAVHE_ROOK performs one of the matrix-vector operations
+*> x := A*x or x := A^H*x,
+*> where x is an N element vector and A is one of the factors
+*> from the Hermitian factorization computed by CHETRF_ROOK.
+*>
+*> CHETRF_ROOK produces a factorization of the form
+*> U * D * U^H or L * D * L^H,
+*> where U (or L) is a product of permutation and unit upper (lower)
+*> triangular matrices, U^H (or L^H) is the conjugate transpose of
+*> U (or L), and D is Hermitian and block diagonal with 1 x 1 and
+*> 2 x 2 diagonal blocks. The multipliers for the transformations
+*> and the upper or lower triangular parts of the diagonal blocks
+*> are stored in the leading upper or lower triangle of the 2-D
+*> array A.
+*>
+*> If TRANS = 'N', multiplies by U or U * D (or L or L * D)
+*> If TRANS = 'T', multiplies by U' or D * U' (or L' or D * L')
+*> If TRANS = 'C', multiplies by U' or D * U' (or L' or D * L')
+*> \endverbatim
+*
+* Arguments:
+* ==========
+*
+*> \param[in] UPLO
+*> \verbatim
+*> UPLO is CHARACTER*1
+*> Specifies whether the factor stored in A is upper or lower
+*> triangular.
+*> = 'U': Upper triangular
+*> = 'L': Lower triangular
+*> \endverbatim
+*>
+*> \param[in] TRANS
+*> \verbatim
+*> TRANS is CHARACTER*1
+*> Specifies the operation to be performed:
+*> = 'N': x := A*x
+*> = 'T': x := A^H*x
+*> = 'C': x := A^H*x
+*> \endverbatim
+*>
+*> \param[in] DIAG
+*> \verbatim
+*> DIAG is CHARACTER*1
+*> Specifies whether or not the diagonal blocks are unit
+*> matrices. If the diagonal blocks are assumed to be unit,
+*> then A = U or A = L, otherwise A = U*D or A = L*D.
+*> = 'U': Diagonal blocks are assumed to be unit matrices.
+*> = 'N': Diagonal blocks are assumed to be non-unit matrices.
+*> \endverbatim
+*>
+*> \param[in] N
+*> \verbatim
+*> N is INTEGER
+*> The number of rows and columns 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 vectors
+*> x to be multiplied by A. 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 CHETRF_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 CHETRF_ROOK.
+*> 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 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: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, B contains NRHS vectors of length N.
+*> On exit, B is overwritten with the product A * B.
+*> \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 = -k, the k-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 April 2013
+*
+*> \ingroup complex_lin
+*
+* =====================================================================
+ SUBROUTINE CLAVHE_ROOK( UPLO, TRANS, DIAG, N, NRHS, A, LDA, IPIV,
+ $ B, LDB, INFO )
+*
+* -- LAPACK test 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 DIAG, TRANS, 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 NOUNIT
+ INTEGER J, K, KP
+ COMPLEX D11, D12, D21, D22, T1, T2
+* ..
+* .. External Functions ..
+ LOGICAL LSAME
+ EXTERNAL LSAME
+* ..
+* .. External Subroutines ..
+ EXTERNAL CGEMV, CGERU, CLACGV, CSCAL, CSWAP, XERBLA
+* ..
+* .. Intrinsic Functions ..
+ INTRINSIC ABS, CONJG, MAX
+* ..
+* .. Executable Statements ..
+*
+* Test the input parameters.
+*
+ INFO = 0
+ IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
+ INFO = -1
+ ELSE IF( .NOT.LSAME( TRANS, 'N' ) .AND. .NOT.LSAME( TRANS, 'C' ) )
+ $ THEN
+ INFO = -2
+ ELSE IF( .NOT.LSAME( DIAG, 'U' ) .AND. .NOT.LSAME( DIAG, 'N' ) )
+ $ THEN
+ INFO = -3
+ ELSE IF( N.LT.0 ) THEN
+ INFO = -4
+ ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
+ INFO = -6
+ ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
+ INFO = -9
+ END IF
+ IF( INFO.NE.0 ) THEN
+ CALL XERBLA( 'CLAVHE_ROOK ', -INFO )
+ RETURN
+ END IF
+*
+* Quick return if possible.
+*
+ IF( N.EQ.0 )
+ $ RETURN
+*
+ NOUNIT = LSAME( DIAG, 'N' )
+*------------------------------------------
+*
+* Compute B := A * B (No transpose)
+*
+*------------------------------------------
+ IF( LSAME( TRANS, 'N' ) ) THEN
+*
+* Compute B := U*B
+* where U = P(m)*inv(U(m))* ... *P(1)*inv(U(1))
+*
+ IF( LSAME( UPLO, 'U' ) ) THEN
+*
+* Loop forward applying the transformations.
+*
+ K = 1
+ 10 CONTINUE
+ IF( K.GT.N )
+ $ GO TO 30
+ IF( IPIV( K ).GT.0 ) THEN
+*
+* 1 x 1 pivot block
+*
+* Multiply by the diagonal element if forming U * D.
+*
+ IF( NOUNIT )
+ $ CALL CSCAL( NRHS, A( K, K ), B( K, 1 ), LDB )
+*
+* Multiply by P(K) * inv(U(K)) if K > 1.
+*
+ IF( K.GT.1 ) THEN
+*
+* Apply the transformation.
+*
+ CALL CGERU( K-1, NRHS, CONE, A( 1, K ), 1, B( K, 1 ),
+ $ LDB, B( 1, 1 ), LDB )
+*
+* Interchange if P(K) != I.
+*
+ KP = IPIV( K )
+ IF( KP.NE.K )
+ $ CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+ END IF
+ K = K + 1
+ ELSE
+*
+* 2 x 2 pivot block
+*
+* Multiply by the diagonal block if forming U * D.
+*
+ IF( NOUNIT ) THEN
+ D11 = A( K, K )
+ D22 = A( K+1, K+1 )
+ D12 = A( K, K+1 )
+ D21 = CONJG( D12 )
+ DO 20 J = 1, NRHS
+ T1 = B( K, J )
+ T2 = B( K+1, J )
+ B( K, J ) = D11*T1 + D12*T2
+ B( K+1, J ) = D21*T1 + D22*T2
+ 20 CONTINUE
+ END IF
+*
+* Multiply by P(K) * inv(U(K)) if K > 1.
+*
+ IF( K.GT.1 ) THEN
+*
+* Apply the transformations.
+*
+ CALL CGERU( K-1, NRHS, CONE, A( 1, K ), 1, B( K, 1 ),
+ $ LDB, B( 1, 1 ), LDB )
+ CALL CGERU( K-1, NRHS, CONE, A( 1, K+1 ), 1,
+ $ B( K+1, 1 ), LDB, B( 1, 1 ), LDB )
+*
+* Interchange if a permutation was applied at the
+* K-th step of the factorization.
+*
+* Swap the first of pair with IMAXth
+*
+ KP = ABS( IPIV( K ) )
+ IF( KP.NE.K )
+ $ CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+*
+* NOW swap the first of pair with Pth
+*
+ KP = ABS( IPIV( K+1 ) )
+ IF( KP.NE.K+1 )
+ $ CALL CSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ),
+ $ LDB )
+ END IF
+ K = K + 2
+ END IF
+ GO TO 10
+ 30 CONTINUE
+*
+* Compute B := L*B
+* where L = P(1)*inv(L(1))* ... *P(m)*inv(L(m)) .
+*
+ ELSE
+*
+* Loop backward applying the transformations to B.
+*
+ K = N
+ 40 CONTINUE
+ IF( K.LT.1 )
+ $ GO TO 60
+*
+* Test the pivot index. If greater than zero, a 1 x 1
+* pivot was used, otherwise a 2 x 2 pivot was used.
+*
+ IF( IPIV( K ).GT.0 ) THEN
+*
+* 1 x 1 pivot block:
+*
+* Multiply by the diagonal element if forming L * D.
+*
+ IF( NOUNIT )
+ $ CALL CSCAL( NRHS, A( K, K ), B( K, 1 ), LDB )
+*
+* Multiply by P(K) * inv(L(K)) if K < N.
+*
+ IF( K.NE.N ) THEN
+ KP = IPIV( K )
+*
+* Apply the transformation.
+*
+ CALL CGERU( N-K, NRHS, CONE, A( K+1, K ), 1,
+ $ B( K, 1 ), LDB, B( K+1, 1 ), LDB )
+*
+* Interchange if a permutation was applied at the
+* K-th step of the factorization.
+*
+ IF( KP.NE.K )
+ $ CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+ END IF
+ K = K - 1
+*
+ ELSE
+*
+* 2 x 2 pivot block:
+*
+* Multiply by the diagonal block if forming L * D.
+*
+ IF( NOUNIT ) THEN
+ D11 = A( K-1, K-1 )
+ D22 = A( K, K )
+ D21 = A( K, K-1 )
+ D12 = CONJG( D21 )
+ DO 50 J = 1, NRHS
+ T1 = B( K-1, J )
+ T2 = B( K, J )
+ B( K-1, J ) = D11*T1 + D12*T2
+ B( K, J ) = D21*T1 + D22*T2
+ 50 CONTINUE
+ END IF
+*
+* Multiply by P(K) * inv(L(K)) if K < N.
+*
+ IF( K.NE.N ) THEN
+*
+* Apply the transformation.
+*
+ CALL CGERU( N-K, NRHS, CONE, A( K+1, K ), 1,
+ $ B( K, 1 ), LDB, B( K+1, 1 ), LDB )
+ CALL CGERU( N-K, NRHS, CONE, A( K+1, K-1 ), 1,
+ $ B( K-1, 1 ), LDB, B( K+1, 1 ), LDB )
+*
+* Interchange if a permutation was applied at the
+* K-th step of the factorization.
+*
+*
+* Swap the second of pair with IMAXth
+*
+ KP = ABS( IPIV( K ) )
+ IF( KP.NE.K )
+ $ CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+*
+* NOW swap the first of pair with Pth
+*
+ KP = ABS( IPIV( K-1 ) )
+ IF( KP.NE.K-1 )
+ $ CALL CSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ),
+ $ LDB )
+*
+ END IF
+ K = K - 2
+ END IF
+ GO TO 40
+ 60 CONTINUE
+ END IF
+*--------------------------------------------------
+*
+* Compute B := A^H * B (conjugate transpose)
+*
+*--------------------------------------------------
+ ELSE
+*
+* Form B := U^H*B
+* where U = P(m)*inv(U(m))* ... *P(1)*inv(U(1))
+* and U^H = inv(U^H(1))*P(1)* ... *inv(U^H(m))*P(m)
+*
+ IF( LSAME( UPLO, 'U' ) ) THEN
+*
+* Loop backward applying the transformations.
+*
+ K = N
+ 70 IF( K.LT.1 )
+ $ GO TO 90
+*
+* 1 x 1 pivot block.
+*
+ IF( IPIV( K ).GT.0 ) THEN
+ IF( K.GT.1 ) THEN
+*
+* Interchange if P(K) != I.
+*
+ KP = IPIV( K )
+ IF( KP.NE.K )
+ $ CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+*
+* Apply the transformation
+* y = y - B' conjg(x),
+* where x is a column of A and y is a row of B.
+*
+ CALL CLACGV( NRHS, B( K, 1 ), LDB )
+ CALL CGEMV( 'Conjugate', K-1, NRHS, CONE, B, LDB,
+ $ A( 1, K ), 1, CONE, B( K, 1 ), LDB )
+ CALL CLACGV( NRHS, B( K, 1 ), LDB )
+ END IF
+ IF( NOUNIT )
+ $ CALL CSCAL( NRHS, A( K, K ), B( K, 1 ), LDB )
+ K = K - 1
+*
+* 2 x 2 pivot block.
+*
+ ELSE
+ IF( K.GT.2 ) THEN
+*
+* Swap the second of pair with Pth
+*
+ KP = ABS( IPIV( K ) )
+ IF( KP.NE.K )
+ $ CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+*
+* Now swap the first of pair with IMAX(r)th
+*
+ KP = ABS( IPIV( K-1 ) )
+ IF( KP.NE.K-1 )
+ $ CALL CSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ),
+ $ LDB )
+*
+* Apply the transformations
+* y = y - B' conjg(x),
+* where x is a block column of A and y is a block
+* row of B.
+*
+ CALL CLACGV( NRHS, B( K, 1 ), LDB )
+ CALL CGEMV( 'Conjugate', K-2, NRHS, CONE, B, LDB,
+ $ A( 1, K ), 1, CONE, B( K, 1 ), LDB )
+ CALL CLACGV( NRHS, B( K, 1 ), LDB )
+*
+ CALL CLACGV( NRHS, B( K-1, 1 ), LDB )
+ CALL CGEMV( 'Conjugate', K-2, NRHS, CONE, B, LDB,
+ $ A( 1, K-1 ), 1, CONE, B( K-1, 1 ), LDB )
+ CALL CLACGV( NRHS, B( K-1, 1 ), LDB )
+ END IF
+*
+* Multiply by the diagonal block if non-unit.
+*
+ IF( NOUNIT ) THEN
+ D11 = A( K-1, K-1 )
+ D22 = A( K, K )
+ D12 = A( K-1, K )
+ D21 = CONJG( D12 )
+ DO 80 J = 1, NRHS
+ T1 = B( K-1, J )
+ T2 = B( K, J )
+ B( K-1, J ) = D11*T1 + D12*T2
+ B( K, J ) = D21*T1 + D22*T2
+ 80 CONTINUE
+ END IF
+ K = K - 2
+ END IF
+ GO TO 70
+ 90 CONTINUE
+*
+* Form B := L^H*B
+* where L = P(1)*inv(L(1))* ... *P(m)*inv(L(m))
+* and L^H = inv(L^H(m))*P(m)* ... *inv(L^H(1))*P(1)
+*
+ ELSE
+*
+* Loop forward applying the L-transformations.
+*
+ K = 1
+ 100 CONTINUE
+ IF( K.GT.N )
+ $ GO TO 120
+*
+* 1 x 1 pivot block
+*
+ IF( IPIV( K ).GT.0 ) THEN
+ IF( K.LT.N ) THEN
+*
+* Interchange if P(K) != I.
+*
+ KP = IPIV( K )
+ IF( KP.NE.K )
+ $ CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+*
+* Apply the transformation
+*
+ CALL CLACGV( NRHS, B( K, 1 ), LDB )
+ CALL CGEMV( 'Conjugate', N-K, NRHS, CONE, B( K+1, 1 ),
+ $ LDB, A( K+1, K ), 1, CONE, B( K, 1 ), LDB )
+ CALL CLACGV( NRHS, B( K, 1 ), LDB )
+ END IF
+ IF( NOUNIT )
+ $ CALL CSCAL( NRHS, A( K, K ), B( K, 1 ), LDB )
+ K = K + 1
+*
+* 2 x 2 pivot block.
+*
+ ELSE
+ IF( K.LT.N-1 ) THEN
+*
+* Swap the first of pair with Pth
+*
+ KP = ABS( IPIV( K ) )
+ IF( KP.NE.K )
+ $ CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+*
+* Now swap the second of pair with IMAX(r)th
+*
+ KP = ABS( IPIV( K+1 ) )
+ IF( KP.NE.K+1 )
+ $ CALL CSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ),
+ $ LDB )
+*
+* Apply the transformation
+*
+ CALL CLACGV( NRHS, B( K+1, 1 ), LDB )
+ CALL CGEMV( 'Conjugate', N-K-1, NRHS, CONE,
+ $ B( K+2, 1 ), LDB, A( K+2, K+1 ), 1, CONE,
+ $ B( K+1, 1 ), LDB )
+ CALL CLACGV( NRHS, B( K+1, 1 ), LDB )
+*
+ CALL CLACGV( NRHS, B( K, 1 ), LDB )
+ CALL CGEMV( 'Conjugate', N-K-1, NRHS, CONE,
+ $ B( K+2, 1 ), LDB, A( K+2, K ), 1, CONE,
+ $ B( K, 1 ), LDB )
+ CALL CLACGV( NRHS, B( K, 1 ), LDB )
+ END IF
+*
+* Multiply by the diagonal block if non-unit.
+*
+ IF( NOUNIT ) THEN
+ D11 = A( K, K )
+ D22 = A( K+1, K+1 )
+ D21 = A( K+1, K )
+ D12 = CONJG( D21 )
+ DO 110 J = 1, NRHS
+ T1 = B( K, J )
+ T2 = B( K+1, J )
+ B( K, J ) = D11*T1 + D12*T2
+ B( K+1, J ) = D21*T1 + D22*T2
+ 110 CONTINUE
+ END IF
+ K = K + 2
+ END IF
+ GO TO 100
+ 120 CONTINUE
+ END IF
+*
+ END IF
+ RETURN
+*
+* End of CLAVHE_ROOK
+*
+ END
--- /dev/null
+*> \brief \b ZCHKHE_ROOK
+*
+* =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+* http://www.netlib.org/lapack/explore-html/
+*
+* Definition:
+* ===========
+*
+* SUBROUTINE ZCHKHE_ROOK( DOTYPE, NN, NVAL, NNB, NBVAL, NNS, NSVAL,
+* THRESH, TSTERR, NMAX, A, AFAC, AINV, B, X,
+* XACT, WORK, RWORK, IWORK, NOUT )
+*
+* .. Scalar Arguments ..
+* LOGICAL TSTERR
+* INTEGER NMAX, NN, NNB, NNS, NOUT
+* DOUBLE PRECISION THRESH
+* ..
+* .. Array Arguments ..
+* LOGICAL DOTYPE( * )
+* INTEGER IWORK( * ), NBVAL( * ), NSVAL( * ), NVAL( * )
+* DOUBLE PRECISION RWORK( * )
+* COMPLEX*16 A( * ), AFAC( * ), AINV( * ), B( * ),
+* $ WORK( * ), X( * ), XACT( * )
+* ..
+*
+*
+*> \par Purpose:
+* =============
+*>
+*> \verbatim
+*>
+*> ZCHKHE_ROOK tests ZHETRF_ROOK, -TRI_ROOK, -TRS_ROOK,
+*> and -CON_ROOK.
+*> \endverbatim
+*
+* Arguments:
+* ==========
+*
+*> \param[in] DOTYPE
+*> \verbatim
+*> DOTYPE is LOGICAL array, dimension (NTYPES)
+*> The matrix types to be used for testing. Matrices of type j
+*> (for 1 <= j <= NTYPES) are used for testing if DOTYPE(j) =
+*> .TRUE.; if DOTYPE(j) = .FALSE., then type j is not used.
+*> \endverbatim
+*>
+*> \param[in] NN
+*> \verbatim
+*> NN is INTEGER
+*> The number of values of N contained in the vector NVAL.
+*> \endverbatim
+*>
+*> \param[in] NVAL
+*> \verbatim
+*> NVAL is INTEGER array, dimension (NN)
+*> The values of the matrix dimension N.
+*> \endverbatim
+*>
+*> \param[in] NNB
+*> \verbatim
+*> NNB is INTEGER
+*> The number of values of NB contained in the vector NBVAL.
+*> \endverbatim
+*>
+*> \param[in] NBVAL
+*> \verbatim
+*> NBVAL is INTEGER array, dimension (NBVAL)
+*> The values of the blocksize NB.
+*> \endverbatim
+*>
+*> \param[in] NNS
+*> \verbatim
+*> NNS is INTEGER
+*> The number of values of NRHS contained in the vector NSVAL.
+*> \endverbatim
+*>
+*> \param[in] NSVAL
+*> \verbatim
+*> NSVAL is INTEGER array, dimension (NNS)
+*> The values of the number of right hand sides NRHS.
+*> \endverbatim
+*>
+*> \param[in] THRESH
+*> \verbatim
+*> THRESH is DOUBLE PRECISION
+*> The threshold value for the test ratios. A result is
+*> included in the output file if RESULT >= THRESH. To have
+*> every test ratio printed, use THRESH = 0.
+*> \endverbatim
+*>
+*> \param[in] TSTERR
+*> \verbatim
+*> TSTERR is LOGICAL
+*> Flag that indicates whether error exits are to be tested.
+*> \endverbatim
+*>
+*> \param[in] NMAX
+*> \verbatim
+*> NMAX is INTEGER
+*> The maximum value permitted for N, used in dimensioning the
+*> work arrays.
+*> \endverbatim
+*>
+*> \param[out] A
+*> \verbatim
+*> A is COMPLEX*16 array, dimension (NMAX*NMAX)
+*> \endverbatim
+*>
+*> \param[out] AFAC
+*> \verbatim
+*> AFAC is COMPLEX*16 array, dimension (NMAX*NMAX)
+*> \endverbatim
+*>
+*> \param[out] AINV
+*> \verbatim
+*> AINV is COMPLEX*16 array, dimension (NMAX*NMAX)
+*> \endverbatim
+*>
+*> \param[out] B
+*> \verbatim
+*> B is COMPLEX*16 array, dimension (NMAX*NSMAX)
+*> where NSMAX is the largest entry in NSVAL.
+*> \endverbatim
+*>
+*> \param[out] X
+*> \verbatim
+*> X is COMPLEX*16 array, dimension (NMAX*NSMAX)
+*> \endverbatim
+*>
+*> \param[out] XACT
+*> \verbatim
+*> XACT is COMPLEX*16 array, dimension (NMAX*NSMAX)
+*> \endverbatim
+*>
+*> \param[out] WORK
+*> \verbatim
+*> WORK is COMPLEX*16 array, dimension
+*> (NMAX*max(3,NSMAX))
+*> \endverbatim
+*>
+*> \param[out] RWORK
+*> \verbatim
+*> RWORK is DOUBLE PRECISION array, dimension
+*> (max(NMAX,2*NSMAX))
+*> \endverbatim
+*>
+*> \param[out] IWORK
+*> \verbatim
+*> IWORK is INTEGER array, dimension (2*NMAX)
+*> \endverbatim
+*>
+*> \param[in] NOUT
+*> \verbatim
+*> NOUT is INTEGER
+*> The unit number for output.
+*> \endverbatim
+*
+* Authors:
+* ========
+*
+*> \author Univ. of Tennessee
+*> \author Univ. of California Berkeley
+*> \author Univ. of Colorado Denver
+*> \author NAG Ltd.
+*
+*> \date April 2013
+*
+*> \ingroup complex16_lin
+*
+* =====================================================================
+ SUBROUTINE ZCHKHE_ROOK( DOTYPE, NN, NVAL, NNB, NBVAL, NNS, NSVAL,
+ $ THRESH, TSTERR, NMAX, A, AFAC, AINV, B, X,
+ $ XACT, WORK, RWORK, IWORK, NOUT )
+*
+* -- LAPACK test routine (version 3.4.1) --
+* -- LAPACK is a software package provided by Univ. of Tennessee, --
+* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+* April 2012
+*
+* .. Scalar Arguments ..
+ LOGICAL TSTERR
+ INTEGER NMAX, NN, NNB, NNS, NOUT
+ DOUBLE PRECISION THRESH
+* ..
+* .. Array Arguments ..
+ LOGICAL DOTYPE( * )
+ INTEGER IWORK( * ), NBVAL( * ), NSVAL( * ), NVAL( * )
+ DOUBLE PRECISION RWORK( * )
+ COMPLEX*16 A( * ), AFAC( * ), AINV( * ), B( * ),
+ $ WORK( * ), X( * ), XACT( * )
+* ..
+*
+* =====================================================================
+*
+* .. Parameters ..
+ DOUBLE PRECISION ZERO, ONE
+ PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 )
+ DOUBLE PRECISION ONEHALF
+ PARAMETER ( ONEHALF = 0.5E+0 )
+ DOUBLE PRECISION EIGHT, SEVTEN
+ PARAMETER ( EIGHT = 8.0E+0, SEVTEN = 17.0E+0 )
+ COMPLEX CZERO
+ PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ) )
+ INTEGER NTYPES
+ PARAMETER ( NTYPES = 10 )
+ INTEGER NTESTS
+ PARAMETER ( NTESTS = 7 )
+* ..
+* .. Local Scalars ..
+ LOGICAL TRFCON, ZEROT
+ CHARACTER DIST, TYPE, UPLO, XTYPE
+ CHARACTER*3 PATH, MATPATH
+ INTEGER I, I1, I2, IMAT, IN, INB, INFO, IOFF, IRHS,
+ $ ITEMP, ITEMP2, IUPLO, IZERO, J, K, KL, KU, LDA,
+ $ LWORK, MODE, N, NB, NERRS, NFAIL, NIMAT, NRHS,
+ $ NRUN, NT
+ DOUBLE PRECISION ALPHA, ANORM, CNDNUM, CONST, LAM_MAX, LAM_MIN,
+ $ RCOND, RCONDC, STEMP
+* ..
+* .. Local Arrays ..
+ CHARACTER UPLOS( 2 )
+ INTEGER ISEED( 4 ), ISEEDY( 4 )
+ DOUBLE PRECISION RESULT( NTESTS )
+ COMPLEX*16 BLOCK( 2, 2 ), CDUMMY( 1 )
+* ..
+* .. External Functions ..
+ DOUBLE PRECISION ZLANGE, ZLANHE, DGET06
+ EXTERNAL ZLANGE, ZLANHE, DGET06
+* ..
+* .. External Subroutines ..
+ EXTERNAL ALAERH, ALAHD, ALASUM, ZERRHE, ZHEEVX, ZGET04,
+ $ ZLACPY, ZLARHS, ZLATB4, ZLATMS, ZPOT02,
+ $ ZPOT03, ZHECON_ROOK, ZHET01_ROOK, ZHETRF_ROOK,
+ $ ZHETRI_ROOK, ZHETRS_ROOK, XLAENV
+* ..
+* .. Intrinsic Functions ..
+ INTRINSIC ABS, MAX, MIN, SQRT
+* ..
+* .. Scalars in Common ..
+ LOGICAL LERR, OK
+ CHARACTER*32 SRNAMT
+ INTEGER INFOT, NUNIT
+* ..
+* .. Common blocks ..
+ COMMON / INFOC / INFOT, NUNIT, OK, LERR
+ COMMON / SRNAMC / SRNAMT
+* ..
+* .. Data statements ..
+ DATA ISEEDY / 1988, 1989, 1990, 1991 /
+ DATA UPLOS / 'U', 'L' /
+* ..
+* .. Executable Statements ..
+*
+* Initialize constants and the random number seed.
+*
+ ALPHA = ( ONE+SQRT( SEVTEN ) ) / EIGHT
+*
+* Test path
+*
+ PATH( 1: 1 ) = 'Zomplex precision'
+ PATH( 2: 3 ) = 'HR'
+*
+* Path to generate matrices
+*
+ MATPATH( 1: 1 ) = 'Zomplex precision'
+ MATPATH( 2: 3 ) = 'HE'
+*
+ NRUN = 0
+ NFAIL = 0
+ NERRS = 0
+ DO 10 I = 1, 4
+ ISEED( I ) = ISEEDY( I )
+ 10 CONTINUE
+*
+* Test the error exits
+*
+ IF( TSTERR )
+ $ CALL ZERRHE( PATH, NOUT )
+ INFOT = 0
+*
+* Set the minimum block size for which the block routine should
+* be used, which will be later returned by ILAENV
+*
+ CALL XLAENV( 2, 2 )
+*
+* Do for each value of N in NVAL
+*
+ DO 270 IN = 1, NN
+ N = NVAL( IN )
+ LDA = MAX( N, 1 )
+ XTYPE = 'N'
+ NIMAT = NTYPES
+ IF( N.LE.0 )
+ $ NIMAT = 1
+*
+ IZERO = 0
+*
+* Do for each value of matrix type IMAT
+*
+ DO 260 IMAT = 1, NIMAT
+*
+* Do the tests only if DOTYPE( IMAT ) is true.
+*
+ IF( .NOT.DOTYPE( IMAT ) )
+ $ GO TO 260
+*
+* Skip types 3, 4, 5, or 6 if the matrix size is too small.
+*
+ ZEROT = IMAT.GE.3 .AND. IMAT.LE.6
+ IF( ZEROT .AND. N.LT.IMAT-2 )
+ $ GO TO 260
+*
+* Do first for UPLO = 'U', then for UPLO = 'L'
+*
+ DO 250 IUPLO = 1, 2
+ UPLO = UPLOS( IUPLO )
+*
+* Begin generate the test matrix A.
+*
+* Set up parameters with ZLATB4 for the matrix generator
+* based on the type of matrix to be generated.
+*
+ CALL ZLATB4( MATPATH, IMAT, N, N, TYPE, KL, KU, ANORM,
+ $ MODE, CNDNUM, DIST )
+*
+* Generate a matrix with ZLATMS.
+*
+ SRNAMT = 'ZLATMS'
+ CALL ZLATMS( N, N, DIST, ISEED, TYPE, RWORK, MODE,
+ $ CNDNUM, ANORM, KL, KU, UPLO, A, LDA,
+ $ WORK, INFO )
+*
+* Check error code from ZLATMS and handle error.
+*
+ IF( INFO.NE.0 ) THEN
+ CALL ALAERH( PATH, 'ZLATMS', INFO, 0, UPLO, N, N,
+ $ -1, -1, -1, IMAT, NFAIL, NERRS, NOUT )
+*
+* Skip all tests for this generated matrix
+*
+ GO TO 250
+ END IF
+*
+* For matrix types 3-6, zero one or more rows and
+* columns of the matrix to test that INFO is returned
+* correctly.
+*
+ IF( ZEROT ) THEN
+ IF( IMAT.EQ.3 ) THEN
+ IZERO = 1
+ ELSE IF( IMAT.EQ.4 ) THEN
+ IZERO = N
+ ELSE
+ IZERO = N / 2 + 1
+ END IF
+*
+ IF( IMAT.LT.6 ) THEN
+*
+* Set row and column IZERO to zero.
+*
+ IF( IUPLO.EQ.1 ) THEN
+ IOFF = ( IZERO-1 )*LDA
+ DO 20 I = 1, IZERO - 1
+ A( IOFF+I ) = CZERO
+ 20 CONTINUE
+ IOFF = IOFF + IZERO
+ DO 30 I = IZERO, N
+ A( IOFF ) = CZERO
+ IOFF = IOFF + LDA
+ 30 CONTINUE
+ ELSE
+ IOFF = IZERO
+ DO 40 I = 1, IZERO - 1
+ A( IOFF ) = CZERO
+ IOFF = IOFF + LDA
+ 40 CONTINUE
+ IOFF = IOFF - IZERO
+ DO 50 I = IZERO, N
+ A( IOFF+I ) = CZERO
+ 50 CONTINUE
+ END IF
+ ELSE
+ IOFF = 0
+ IF( IUPLO.EQ.1 ) THEN
+*
+* Set the first IZERO rows and columns to zero.
+*
+ DO 70 J = 1, N
+ I2 = MIN( J, IZERO )
+ DO 60 I = 1, I2
+ A( IOFF+I ) = CZERO
+ 60 CONTINUE
+ IOFF = IOFF + LDA
+ 70 CONTINUE
+ ELSE
+*
+* Set the last IZERO rows and columns to zero.
+*
+ DO 90 J = 1, N
+ I1 = MAX( J, IZERO )
+ DO 80 I = I1, N
+ A( IOFF+I ) = CZERO
+ 80 CONTINUE
+ IOFF = IOFF + LDA
+ 90 CONTINUE
+ END IF
+ END IF
+ ELSE
+ IZERO = 0
+ END IF
+*
+* End generate the test matrix A.
+*
+*
+* Do for each value of NB in NBVAL
+*
+ DO 240 INB = 1, NNB
+*
+* Set the optimal blocksize, which will be later
+* returned by ILAENV.
+*
+ NB = NBVAL( INB )
+ CALL XLAENV( 1, NB )
+*
+* Copy the test matrix A into matrix AFAC which
+* will be factorized in place. This is needed to
+* preserve the test matrix A for subsequent tests.
+*
+ CALL ZLACPY( UPLO, N, N, A, LDA, AFAC, LDA )
+*
+* Compute the L*D*L**T or U*D*U**T factorization of the
+* matrix. IWORK stores details of the interchanges and
+* the block structure of D. AINV is a work array for
+* block factorization, LWORK is the length of AINV.
+*
+ LWORK = MAX( 2, NB )*LDA
+ SRNAMT = 'ZHETRF_ROOK'
+ CALL ZHETRF_ROOK( UPLO, N, AFAC, LDA, IWORK, AINV,
+ $ LWORK, INFO )
+*
+* Adjust the expected value of INFO to account for
+* pivoting.
+*
+ K = IZERO
+ IF( K.GT.0 ) THEN
+ 100 CONTINUE
+ IF( IWORK( K ).LT.0 ) THEN
+ IF( IWORK( K ).NE.-K ) THEN
+ K = -IWORK( K )
+ GO TO 100
+ END IF
+ ELSE IF( IWORK( K ).NE.K ) THEN
+ K = IWORK( K )
+ GO TO 100
+ END IF
+ END IF
+*
+* Check error code from ZHETRF_ROOK and handle error.
+*
+ IF( INFO.NE.K)
+ $ CALL ALAERH( PATH, 'ZHETRF_ROOK', INFO, K,
+ $ UPLO, N, N, -1, -1, NB, IMAT,
+ $ NFAIL, NERRS, NOUT )
+*
+* Set the condition estimate flag if the INFO is not 0.
+*
+ IF( INFO.NE.0 ) THEN
+ TRFCON = .TRUE.
+ ELSE
+ TRFCON = .FALSE.
+ END IF
+*
+*+ TEST 1
+* Reconstruct matrix from factors and compute residual.
+*
+ CALL ZHET01_ROOK( UPLO, N, A, LDA, AFAC, LDA, IWORK,
+ $ AINV, LDA, RWORK, RESULT( 1 ) )
+ NT = 1
+*
+*+ TEST 2
+* Form the inverse and compute the residual,
+* if the factorization was competed without INFO > 0
+* (i.e. there is no zero rows and columns).
+* Do it only for the first block size.
+*
+ IF( INB.EQ.1 .AND. .NOT.TRFCON ) THEN
+ CALL ZLACPY( UPLO, N, N, AFAC, LDA, AINV, LDA )
+ SRNAMT = 'ZHETRI_ROOK'
+ CALL ZHETRI_ROOK( UPLO, N, AINV, LDA, IWORK, WORK,
+ $ INFO )
+*
+* Check error code from ZHETRI_ROOK and handle error.
+*
+ IF( INFO.NE.0 )
+ $ CALL ALAERH( PATH, 'ZHETRI_ROOK', INFO, -1,
+ $ UPLO, N, N, -1, -1, -1, IMAT,
+ $ NFAIL, NERRS, NOUT )
+*
+* Compute the residual for a Hermitian matrix times
+* its inverse.
+*
+ CALL ZPOT03( UPLO, N, A, LDA, AINV, LDA, WORK, LDA,
+ $ RWORK, RCONDC, RESULT( 2 ) )
+ NT = 2
+ END IF
+*
+* Print information about the tests that did not pass
+* the threshold.
+*
+ DO 110 K = 1, NT
+ IF( RESULT( K ).GE.THRESH ) THEN
+ IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
+ $ CALL ALAHD( NOUT, PATH )
+ WRITE( NOUT, FMT = 9999 )UPLO, N, NB, IMAT, K,
+ $ RESULT( K )
+ NFAIL = NFAIL + 1
+ END IF
+ 110 CONTINUE
+ NRUN = NRUN + NT
+*
+*+ TEST 3
+* Compute largest element in U or L
+*
+ RESULT( 3 ) = ZERO
+ STEMP = ZERO
+*
+ CONST = ( ( ALPHA**2-ONE ) / ( ALPHA**2-ONEHALF ) ) /
+ $ ( ONE-ALPHA )
+*
+ IF( IUPLO.EQ.1 ) THEN
+*
+* Compute largest element in U
+*
+ K = N
+ 120 CONTINUE
+ IF( K.LE.1 )
+ $ GO TO 130
+*
+ IF( IWORK( K ).GT.ZERO ) THEN
+*
+* Get max absolute value from elements
+* in column k in in U
+*
+ STEMP = ZLANGE( 'M', K-1, 1,
+ $ AFAC( ( K-1 )*LDA+1 ), LDA, RWORK )
+ ELSE
+*
+* Get max absolute value from elements
+* in columns k and k-1 in U
+*
+ STEMP = ZLANGE( 'M', K-2, 2,
+ $ AFAC( ( K-2 )*LDA+1 ), LDA, RWORK )
+ K = K - 1
+*
+ END IF
+*
+* STEMP should be bounded by CONST
+*
+ STEMP = STEMP - CONST + THRESH
+ IF( STEMP.GT.RESULT( 3 ) )
+ $ RESULT( 3 ) = STEMP
+*
+ K = K - 1
+*
+ GO TO 120
+ 130 CONTINUE
+*
+ ELSE
+*
+* Compute largest element in L
+*
+ K = 1
+ 140 CONTINUE
+ IF( K.GE.N )
+ $ GO TO 150
+*
+ IF( IWORK( K ).GT.ZERO ) THEN
+*
+* Get max absolute value from elements
+* in column k in L
+*
+ STEMP = ZLANGE( 'M', N-K, 1,
+ $ AFAC( ( K-1 )*LDA+K+1 ), LDA, RWORK )
+ ELSE
+*
+* Get max absolute value from elements
+* in columns k and k+1 in L
+*
+ STEMP = ZLANGE( 'M', N-K-1, 2,
+ $ AFAC( ( K-1 )*LDA+K+2 ), LDA, RWORK )
+ K = K + 1
+*
+ END IF
+*
+* STEMP should be bounded by CONST
+*
+ STEMP = STEMP - CONST + THRESH
+ IF( STEMP.GT.RESULT( 3 ) )
+ $ RESULT( 3 ) = STEMP
+*
+ K = K + 1
+*
+ GO TO 140
+ 150 CONTINUE
+ END IF
+*
+*
+*+ TEST 4
+* Compute largest 2-Norm of 2-by-2 diag blocks
+*
+ RESULT( 4 ) = ZERO
+ STEMP = ZERO
+*
+ CONST = ( ( ALPHA**2-ONE ) / ( ALPHA**2-ONEHALF ) )*
+ $ ( ( ONE + ALPHA ) / ( ONE - ALPHA ) )
+*
+ IF( IUPLO.EQ.1 ) THEN
+*
+* Loop backward for UPLO = 'U'
+*
+ K = N
+ 160 CONTINUE
+ IF( K.LE.1 )
+ $ GO TO 170
+*
+ IF( IWORK( K ).LT.ZERO ) THEN
+*
+* Get the two eigenvalues of a 2-by-2 block,
+* store them in WORK array
+*
+ BLOCK( 1, 1 ) = AFAC( ( K-2 )*LDA+K-1 )
+ BLOCK( 2, 1 ) = AFAC( ( K-2 )*LDA+K )
+ BLOCK( 1, 2 ) = BLOCK( 2, 1 )
+ BLOCK( 2, 2 ) = AFAC( (K-1)*LDA+K )
+*
+ CALL ZHEEVX( 'N', 'N', 'N', 'N', 2, BLOCK,
+ $ 2, WORK, CDUMMY, 1, CDUMMY, 1,
+ $ ITEMP, ITEMP2, RWORK, STEMP,
+ $ RWORK( 3 ), RWORK( 5 ), WORK( 3 ),
+ $ 4, RWORK( 7 ), INFO )
+*
+ LAM_MAX = MAX( ABS( WORK( 1 ) ),
+ $ ABS( WORK( 2 ) ) )
+ LAM_MIN = MIN( ABS( WORK( 1 ) ),
+ $ ABS( WORK( 2 ) ) )
+*
+ STEMP = LAM_MAX / LAM_MIN
+*
+* STEMP should be bounded by CONST
+*
+ STEMP = ABS( STEMP ) - CONST + THRESH
+ IF( STEMP.GT.RESULT( 4 ) )
+ $ RESULT( 4 ) = STEMP
+ K = K - 1
+*
+ END IF
+*
+ K = K - 1
+*
+ GO TO 160
+ 170 CONTINUE
+*
+ ELSE
+*
+* Loop forward for UPLO = 'L'
+*
+ K = 1
+ 180 CONTINUE
+ IF( K.GE.N )
+ $ GO TO 190
+*
+ IF( IWORK( K ).LT.ZERO ) THEN
+*
+* Get the two eigenvalues of a 2-by-2 block,
+* store them in WORK array
+*
+ BLOCK( 1, 1 ) = AFAC( ( K-1 )*LDA+K )
+ BLOCK( 2, 1 ) = AFAC( ( K-1 )*LDA+K+1 )
+ BLOCK( 1, 2 ) = BLOCK( 2, 1 )
+ BLOCK( 2, 2 ) = AFAC( K*LDA+K+1 )
+*
+ CALL ZHEEVX( 'N', 'N', 'N', 'N', 2, BLOCK,
+ $ 2, WORK, CDUMMY, 1, CDUMMY, 1,
+ $ ITEMP, ITEMP2, RWORK, STEMP,
+ $ RWORK( 3 ), RWORK( 5 ), WORK( 3 ),
+ $ 4, RWORK( 7 ), INFO )
+*
+ LAM_MAX = MAX( ABS( WORK( 1 ) ),
+ $ ABS( WORK( 2 ) ) )
+ LAM_MIN = MIN( ABS( WORK( 1 ) ),
+ $ ABS( WORK( 2 ) ) )
+*
+ STEMP = LAM_MAX / LAM_MIN
+*
+* STEMP should be bounded by CONST
+*
+ STEMP = ABS( STEMP ) - CONST + THRESH
+ IF( STEMP.GT.RESULT( 4 ) )
+ $ RESULT( 4 ) = STEMP
+ K = K + 1
+*
+ END IF
+*
+ K = K + 1
+*
+ GO TO 180
+ 190 CONTINUE
+ END IF
+*
+* Print information about the tests that did not pass
+* the threshold.
+*
+ DO 200 K = 3, 4
+ IF( RESULT( K ).GE.THRESH ) THEN
+ IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
+ $ CALL ALAHD( NOUT, PATH )
+ WRITE( NOUT, FMT = 9999 )UPLO, N, NB, IMAT, K,
+ $ RESULT( K )
+ NFAIL = NFAIL + 1
+ END IF
+ 200 CONTINUE
+ NRUN = NRUN + NT
+*
+* Skip the other tests if this is not the first block
+* size.
+*
+ IF( INB.GT.1 )
+ $ GO TO 240
+*
+* Do only the condition estimate if INFO is not 0.
+*
+ IF( TRFCON ) THEN
+ RCONDC = ZERO
+ GO TO 230
+ END IF
+*
+ DO 220 IRHS = 1, NNS
+ NRHS = NSVAL( IRHS )
+*
+* Begin loop over NRHS values
+*
+*
+*+ TEST 5 ( Using TRS_ROOK)
+* Solve and compute residual for A * X = B.
+*
+* Choose a set of NRHS random solution vectors
+* stored in XACT and set up the right hand side B
+*
+ SRNAMT = 'ZLARHS'
+ CALL ZLARHS( MATPATH, XTYPE, UPLO, ' ', N, N,
+ $ KL, KU, NRHS, A, LDA, XACT, LDA,
+ $ B, LDA, ISEED, INFO )
+ CALL ZLACPY( 'Full', N, NRHS, B, LDA, X, LDA )
+*
+ SRNAMT = 'ZHETRS_ROOK'
+ CALL ZHETRS_ROOK( UPLO, N, NRHS, AFAC, LDA, IWORK,
+ $ X, LDA, INFO )
+*
+* Check error code from ZHETRS_ROOK and handle error.
+*
+ IF( INFO.NE.0 )
+ $ CALL ALAERH( PATH, 'ZHETRS_ROOK', INFO, 0,
+ $ UPLO, N, N, -1, -1, NRHS, IMAT,
+ $ NFAIL, NERRS, NOUT )
+*
+ CALL ZLACPY( 'Full', N, NRHS, B, LDA, WORK, LDA )
+*
+* Compute the residual for the solution
+*
+ CALL ZPOT02( UPLO, N, NRHS, A, LDA, X, LDA, WORK,
+ $ LDA, RWORK, RESULT( 5 ) )
+*
+*+ TEST 6
+* Check solution from generated exact solution.
+*
+ CALL ZGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC,
+ $ RESULT( 6 ) )
+*
+* Print information about the tests that did not pass
+* the threshold.
+*
+ DO 210 K = 5, 6
+ IF( RESULT( K ).GE.THRESH ) THEN
+ IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
+ $ CALL ALAHD( NOUT, PATH )
+ WRITE( NOUT, FMT = 9998 )UPLO, N, NRHS,
+ $ IMAT, K, RESULT( K )
+ NFAIL = NFAIL + 1
+ END IF
+ 210 CONTINUE
+ NRUN = NRUN + 2
+*
+* End loop over NRHS values
+*
+ 220 CONTINUE
+*
+*+ TEST 7
+* Get an estimate of RCOND = 1/CNDNUM.
+*
+ 230 CONTINUE
+ ANORM = ZLANHE( '1', UPLO, N, A, LDA, RWORK )
+ SRNAMT = 'ZHECON_ROOK'
+ CALL ZHECON_ROOK( UPLO, N, AFAC, LDA, IWORK, ANORM,
+ $ RCOND, WORK, INFO )
+*
+* Check error code from ZHECON_ROOK and handle error.
+*
+ IF( INFO.NE.0 )
+ $ CALL ALAERH( PATH, 'ZHECON_ROOK', INFO, 0,
+ $ UPLO, N, N, -1, -1, -1, IMAT,
+ $ NFAIL, NERRS, NOUT )
+*
+* Compute the test ratio to compare to values of RCOND
+*
+ RESULT( 7 ) = DGET06( RCOND, RCONDC )
+*
+* Print information about the tests that did not pass
+* the threshold.
+*
+ IF( RESULT( 7 ).GE.THRESH ) THEN
+ IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
+ $ CALL ALAHD( NOUT, PATH )
+ WRITE( NOUT, FMT = 9997 )UPLO, N, IMAT, 7,
+ $ RESULT( 7 )
+ NFAIL = NFAIL + 1
+ END IF
+ NRUN = NRUN + 1
+ 240 CONTINUE
+*
+ 250 CONTINUE
+ 260 CONTINUE
+ 270 CONTINUE
+*
+* Print a summary of the results.
+*
+ CALL ALASUM( PATH, NOUT, NFAIL, NRUN, NERRS )
+*
+ 9999 FORMAT( ' UPLO = ''', A1, ''', N =', I5, ', NB =', I4, ', type ',
+ $ I2, ', test ', I2, ', ratio =', G12.5 )
+ 9998 FORMAT( ' UPLO = ''', A1, ''', N =', I5, ', NRHS=', I3, ', type ',
+ $ I2, ', test(', I2, ') =', G12.5 )
+ 9997 FORMAT( ' UPLO = ''', A1, ''', N =', I5, ',', 10X, ' type ', I2,
+ $ ', test(', I2, ') =', G12.5 )
+ RETURN
+*
+* End of ZCHKHE_ROOK
+*
+ END
--- /dev/null
+*> \brief \b ZHET01_ROOK
+*
+* =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+* http://www.netlib.org/lapack/explore-html/
+*
+* Definition:
+* ===========
+*
+* SUBROUTINE ZHET01_ROOK( UPLO, N, A, LDA, AFAC, LDAFAC, IPIV, C, LDC,
+* RWORK, RESID )
+*
+* .. Scalar Arguments ..
+* CHARACTER UPLO
+* INTEGER LDA, LDAFAC, LDC, N
+* DOUBLE PRECISION RESID
+* ..
+* .. Array Arguments ..
+* INTEGER IPIV( * )
+* DOUBLE PRECISION RWORK( * )
+* COMPLEX*16 A( LDA, * ), AFAC( LDAFAC, * ), C( LDC, * )
+* ..
+*
+*
+*> \par Purpose:
+* =============
+*>
+*> \verbatim
+*>
+*> ZHET01_ROOK reconstructs a complex Hermitian indefinite matrix A from its
+*> block L*D*L' or U*D*U' factorization and computes the residual
+*> norm( C - A ) / ( N * norm(A) * EPS ),
+*> where C is the reconstructed matrix, EPS is the machine epsilon,
+*> L' is the transpose of L, and U' is the transpose of U.
+*> \endverbatim
+*
+* Arguments:
+* ==========
+*
+*> \param[in] UPLO
+*> \verbatim
+*> UPLO is CHARACTER*1
+*> Specifies whether the upper or lower triangular part of the
+*> complex Hermitian matrix A is stored:
+*> = 'U': Upper triangular
+*> = 'L': Lower triangular
+*> \endverbatim
+*>
+*> \param[in] N
+*> \verbatim
+*> N is INTEGER
+*> The number of rows and columns of the matrix A. N >= 0.
+*> \endverbatim
+*>
+*> \param[in] A
+*> \verbatim
+*> A is COMPLEX*16 array, dimension (LDA,N)
+*> The original complex Hermitian matrix A.
+*> \endverbatim
+*>
+*> \param[in] LDA
+*> \verbatim
+*> LDA is INTEGER
+*> The leading dimension of the array A. LDA >= max(1,N)
+*> \endverbatim
+*>
+*> \param[in] AFAC
+*> \verbatim
+*> AFAC is COMPLEX*16 array, dimension (LDAFAC,N)
+*> The factored form of the matrix A. AFAC contains the block
+*> diagonal matrix D and the multipliers used to obtain the
+*> factor L or U from the block L*D*L' or U*D*U' factorization
+*> as computed by CSYTRF_ROOK.
+*> \endverbatim
+*>
+*> \param[in] LDAFAC
+*> \verbatim
+*> LDAFAC is INTEGER
+*> The leading dimension of the array AFAC. LDAFAC >= max(1,N).
+*> \endverbatim
+*>
+*> \param[in] IPIV
+*> \verbatim
+*> IPIV is INTEGER array, dimension (N)
+*> The pivot indices from CSYTRF_ROOK.
+*> \endverbatim
+*>
+*> \param[out] C
+*> \verbatim
+*> C is COMPLEX*16 array, dimension (LDC,N)
+*> \endverbatim
+*>
+*> \param[in] LDC
+*> \verbatim
+*> LDC is INTEGER
+*> The leading dimension of the array C. LDC >= max(1,N).
+*> \endverbatim
+*>
+*> \param[out] RWORK
+*> \verbatim
+*> RWORK is DOUBLE PRECISION array, dimension (N)
+*> \endverbatim
+*>
+*> \param[out] RESID
+*> \verbatim
+*> RESID is DOUBLE PRECISION
+*> If UPLO = 'L', norm(L*D*L' - A) / ( N * norm(A) * EPS )
+*> If UPLO = 'U', norm(U*D*U' - A) / ( N * norm(A) * EPS )
+*> \endverbatim
+*
+* Authors:
+* ========
+*
+*> \author Univ. of Tennessee
+*> \author Univ. of California Berkeley
+*> \author Univ. of Colorado Denver
+*> \author NAG Ltd.
+*
+*> \date April 2013
+*
+*> \ingroup complex16_lin
+*
+* =====================================================================
+ SUBROUTINE ZHET01_ROOK( UPLO, N, A, LDA, AFAC, LDAFAC, IPIV, C,
+ $ LDC, RWORK, RESID )
+*
+* -- LAPACK test 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 LDA, LDAFAC, LDC, N
+ DOUBLE PRECISION RESID
+* ..
+* .. Array Arguments ..
+ INTEGER IPIV( * )
+ DOUBLE PRECISION RWORK( * )
+ COMPLEX*16 A( LDA, * ), AFAC( LDAFAC, * ), C( LDC, * )
+* ..
+*
+* =====================================================================
+*
+* .. Parameters ..
+ DOUBLE PRECISION ZERO, ONE
+ PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 )
+ COMPLEX*16 CZERO, CONE
+ PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ),
+ $ CONE = ( 1.0E+0, 0.0E+0 ) )
+* ..
+* .. Local Scalars ..
+ INTEGER I, INFO, J
+ DOUBLE PRECISION ANORM, EPS
+* ..
+* .. External Functions ..
+ LOGICAL LSAME
+ DOUBLE PRECISION ZLANHE, DLAMCH
+ EXTERNAL LSAME, ZLANHE, DLAMCH
+* ..
+* .. External Subroutines ..
+ EXTERNAL ZLASET, ZLAVHE_ROOK
+* ..
+* .. Intrinsic Functions ..
+ INTRINSIC DIMAG, DBLE
+* ..
+* .. Executable Statements ..
+*
+* Quick exit if N = 0.
+*
+ IF( N.LE.0 ) THEN
+ RESID = ZERO
+ RETURN
+ END IF
+*
+* Determine EPS and the norm of A.
+*
+ EPS = DLAMCH( 'Epsilon' )
+ ANORM = ZLANHE( '1', UPLO, N, A, LDA, RWORK )
+*
+* Check the imaginary parts of the diagonal elements and return with
+* an error code if any are nonzero.
+*
+ DO 10 J = 1, N
+ IF( DIMAG( AFAC( J, J ) ).NE.ZERO ) THEN
+ RESID = ONE / EPS
+ RETURN
+ END IF
+ 10 CONTINUE
+*
+* Initialize C to the identity matrix.
+*
+ CALL ZLASET( 'Full', N, N, CZERO, CONE, C, LDC )
+*
+* Call ZLAVHE_ROOK to form the product D * U' (or D * L' ).
+*
+ CALL ZLAVHE_ROOK( UPLO, 'Conjugate', 'Non-unit', N, N, AFAC,
+ $ LDAFAC, IPIV, C, LDC, INFO )
+*
+* Call ZLAVHE_ROOK again to multiply by U (or L ).
+*
+ CALL ZLAVHE_ROOK( UPLO, 'No transpose', 'Unit', N, N, AFAC,
+ $ LDAFAC, IPIV, C, LDC, INFO )
+*
+* Compute the difference C - A .
+*
+ IF( LSAME( UPLO, 'U' ) ) THEN
+ DO 30 J = 1, N
+ DO 20 I = 1, J - 1
+ C( I, J ) = C( I, J ) - A( I, J )
+ 20 CONTINUE
+ C( J, J ) = C( J, J ) - DBLE( A( J, J ) )
+ 30 CONTINUE
+ ELSE
+ DO 50 J = 1, N
+ C( J, J ) = C( J, J ) - DBLE( A( J, J ) )
+ DO 40 I = J + 1, N
+ C( I, J ) = C( I, J ) - A( I, J )
+ 40 CONTINUE
+ 50 CONTINUE
+ END IF
+*
+* Compute norm( C - A ) / ( N * norm(A) * EPS )
+*
+ RESID = ZLANHE( '1', UPLO, N, C, LDC, RWORK )
+*
+ IF( ANORM.LE.ZERO ) THEN
+ IF( RESID.NE.ZERO )
+ $ RESID = ONE / EPS
+ ELSE
+ RESID = ( ( RESID/DBLE( N ) )/ANORM ) / EPS
+ END IF
+*
+ RETURN
+*
+* End of ZHET01_ROOK
+*
+ END
--- /dev/null
+*> \brief \b ZLAVHE_ROOK
+*
+* =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+* http://www.netlib.org/lapack/explore-html/
+*
+* Definition:
+* ===========
+*
+* SUBROUTINE ZLAVHE_ROOK( UPLO, TRANS, DIAG, N, NRHS, A, LDA, IPIV, B,
+* LDB, INFO )
+*
+* .. Scalar Arguments ..
+* CHARACTER DIAG, TRANS, UPLO
+* INTEGER INFO, LDA, LDB, N, NRHS
+* ..
+* .. Array Arguments ..
+* INTEGER IPIV( * )
+* COMPLEX*16 A( LDA, * ), B( LDB, * )
+* ..
+*
+*
+*> \par Purpose:
+* =============
+*>
+*> \verbatim
+*>
+*> ZLAVHE_ROOK performs one of the matrix-vector operations
+*> x := A*x or x := A^H*x,
+*> where x is an N element vector and A is one of the factors
+*> from the Hermitian factorization computed by CHETRF_ROOK.
+*>
+*> ZHETRF_ROOK produces a factorization of the form
+*> U * D * U^H or L * D * L^H,
+*> where U (or L) is a product of permutation and unit upper (lower)
+*> triangular matrices, U^H (or L^H) is the conjugate transpose of
+*> U (or L), and D is Hermitian and block diagonal with 1 x 1 and
+*> 2 x 2 diagonal blocks. The multipliers for the transformations
+*> and the upper or lower triangular parts of the diagonal blocks
+*> are stored in the leading upper or lower triangle of the 2-D
+*> array A.
+*>
+*> If TRANS = 'N', multiplies by U or U * D (or L or L * D)
+*> If TRANS = 'T', multiplies by U' or D * U' (or L' or D * L')
+*> If TRANS = 'C', multiplies by U' or D * U' (or L' or D * L')
+*> \endverbatim
+*
+* Arguments:
+* ==========
+*
+*> \param[in] UPLO
+*> \verbatim
+*> UPLO is CHARACTER*1
+*> Specifies whether the factor stored in A is upper or lower
+*> triangular.
+*> = 'U': Upper triangular
+*> = 'L': Lower triangular
+*> \endverbatim
+*>
+*> \param[in] TRANS
+*> \verbatim
+*> TRANS is CHARACTER*1
+*> Specifies the operation to be performed:
+*> = 'N': x := A*x
+*> = 'T': x := A^H*x
+*> \endverbatim
+*>
+*> \param[in] DIAG
+*> \verbatim
+*> DIAG is CHARACTER*1
+*> Specifies whether or not the diagonal blocks are unit
+*> matrices. If the diagonal blocks are assumed to be unit,
+*> then A = U or A = L, otherwise A = U*D or A = L*D.
+*> = 'U': Diagonal blocks are assumed to be unit matrices.
+*> = 'N': Diagonal blocks are assumed to be non-unit matrices.
+*> \endverbatim
+*>
+*> \param[in] N
+*> \verbatim
+*> N is INTEGER
+*> The number of rows and columns 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 vectors
+*> x to be multiplied by A. NRHS >= 0.
+*> \endverbatim
+*>
+*> \param[in] A
+*> \verbatim
+*> A is COMPLEX*16 array, dimension (LDA,N)
+*> The block diagonal matrix D and the multipliers used to
+*> obtain the factor U or L as computed by ZHETRF_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 ZHETRF_ROOK.
+*> 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 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: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, B contains NRHS vectors of length N.
+*> On exit, B is overwritten with the product A * B.
+*> \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 = -k, the k-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 April 2013
+*
+*> \ingroup complex16_lin
+*
+* =====================================================================
+ SUBROUTINE ZLAVHE_ROOK( UPLO, TRANS, DIAG, N, NRHS, A, LDA, IPIV,
+ $ B, LDB, INFO )
+*
+* -- LAPACK test 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 DIAG, TRANS, UPLO
+ INTEGER INFO, LDA, LDB, N, NRHS
+* ..
+* .. Array Arguments ..
+ INTEGER IPIV( * )
+ COMPLEX*16 A( LDA, * ), B( LDB, * )
+* ..
+*
+* =====================================================================
+*
+* .. Parameters ..
+ COMPLEX*16 CONE
+ PARAMETER ( CONE = ( 1.0E+0, 0.0E+0 ) )
+* ..
+* .. Local Scalars ..
+ LOGICAL NOUNIT
+ INTEGER J, K, KP
+ COMPLEX*16 D11, D12, D21, D22, T1, T2
+* ..
+* .. External Functions ..
+ LOGICAL LSAME
+ EXTERNAL LSAME
+* ..
+* .. External Subroutines ..
+ EXTERNAL ZGEMV, ZGERU, ZLACGV, ZSCAL, ZSWAP, XERBLA
+* ..
+* .. Intrinsic Functions ..
+ INTRINSIC ABS, DCONJG, MAX
+* ..
+* .. Executable Statements ..
+*
+* Test the input parameters.
+*
+ INFO = 0
+ IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
+ INFO = -1
+ ELSE IF( .NOT.LSAME( TRANS, 'N' ) .AND. .NOT.LSAME( TRANS, 'C' ) )
+ $ THEN
+ INFO = -2
+ ELSE IF( .NOT.LSAME( DIAG, 'U' ) .AND. .NOT.LSAME( DIAG, 'N' ) )
+ $ THEN
+ INFO = -3
+ ELSE IF( N.LT.0 ) THEN
+ INFO = -4
+ ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
+ INFO = -6
+ ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
+ INFO = -9
+ END IF
+ IF( INFO.NE.0 ) THEN
+ CALL XERBLA( 'ZLAVHE_ROOK ', -INFO )
+ RETURN
+ END IF
+*
+* Quick return if possible.
+*
+ IF( N.EQ.0 )
+ $ RETURN
+*
+ NOUNIT = LSAME( DIAG, 'N' )
+*------------------------------------------
+*
+* Compute B := A * B (No transpose)
+*
+*------------------------------------------
+ IF( LSAME( TRANS, 'N' ) ) THEN
+*
+* Compute B := U*B
+* where U = P(m)*inv(U(m))* ... *P(1)*inv(U(1))
+*
+ IF( LSAME( UPLO, 'U' ) ) THEN
+*
+* Loop forward applying the transformations.
+*
+ K = 1
+ 10 CONTINUE
+ IF( K.GT.N )
+ $ GO TO 30
+ IF( IPIV( K ).GT.0 ) THEN
+*
+* 1 x 1 pivot block
+*
+* Multiply by the diagonal element if forming U * D.
+*
+ IF( NOUNIT )
+ $ CALL ZSCAL( NRHS, A( K, K ), B( K, 1 ), LDB )
+*
+* Multiply by P(K) * inv(U(K)) if K > 1.
+*
+ IF( K.GT.1 ) THEN
+*
+* Apply the transformation.
+*
+ CALL ZGERU( K-1, NRHS, CONE, A( 1, K ), 1, B( K, 1 ),
+ $ LDB, B( 1, 1 ), LDB )
+*
+* Interchange if P(K) != I.
+*
+ KP = IPIV( K )
+ IF( KP.NE.K )
+ $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+ END IF
+ K = K + 1
+ ELSE
+*
+* 2 x 2 pivot block
+*
+* Multiply by the diagonal block if forming U * D.
+*
+ IF( NOUNIT ) THEN
+ D11 = A( K, K )
+ D22 = A( K+1, K+1 )
+ D12 = A( K, K+1 )
+ D21 = DCONJG( D12 )
+ DO 20 J = 1, NRHS
+ T1 = B( K, J )
+ T2 = B( K+1, J )
+ B( K, J ) = D11*T1 + D12*T2
+ B( K+1, J ) = D21*T1 + D22*T2
+ 20 CONTINUE
+ END IF
+*
+* Multiply by P(K) * inv(U(K)) if K > 1.
+*
+ IF( K.GT.1 ) THEN
+*
+* Apply the transformations.
+*
+ CALL ZGERU( K-1, NRHS, CONE, A( 1, K ), 1, B( K, 1 ),
+ $ LDB, B( 1, 1 ), LDB )
+ CALL ZGERU( K-1, NRHS, CONE, A( 1, K+1 ), 1,
+ $ B( K+1, 1 ), LDB, B( 1, 1 ), LDB )
+*
+* Interchange if a permutation was applied at the
+* K-th step of the factorization.
+*
+* Swap the first of pair with IMAXth
+*
+ KP = ABS( IPIV( K ) )
+ IF( KP.NE.K )
+ $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+*
+* NOW swap the first of pair with Pth
+*
+ KP = ABS( IPIV( K+1 ) )
+ IF( KP.NE.K+1 )
+ $ CALL ZSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ),
+ $ LDB )
+ END IF
+ K = K + 2
+ END IF
+ GO TO 10
+ 30 CONTINUE
+*
+* Compute B := L*B
+* where L = P(1)*inv(L(1))* ... *P(m)*inv(L(m)) .
+*
+ ELSE
+*
+* Loop backward applying the transformations to B.
+*
+ K = N
+ 40 CONTINUE
+ IF( K.LT.1 )
+ $ GO TO 60
+*
+* Test the pivot index. If greater than zero, a 1 x 1
+* pivot was used, otherwise a 2 x 2 pivot was used.
+*
+ IF( IPIV( K ).GT.0 ) THEN
+*
+* 1 x 1 pivot block:
+*
+* Multiply by the diagonal element if forming L * D.
+*
+ IF( NOUNIT )
+ $ CALL ZSCAL( NRHS, A( K, K ), B( K, 1 ), LDB )
+*
+* Multiply by P(K) * inv(L(K)) if K < N.
+*
+ IF( K.NE.N ) THEN
+ KP = IPIV( K )
+*
+* Apply the transformation.
+*
+ CALL ZGERU( N-K, NRHS, CONE, A( K+1, K ), 1,
+ $ B( K, 1 ), LDB, B( K+1, 1 ), LDB )
+*
+* Interchange if a permutation was applied at the
+* K-th step of the factorization.
+*
+ IF( KP.NE.K )
+ $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+ END IF
+ K = K - 1
+*
+ ELSE
+*
+* 2 x 2 pivot block:
+*
+* Multiply by the diagonal block if forming L * D.
+*
+ IF( NOUNIT ) THEN
+ D11 = A( K-1, K-1 )
+ D22 = A( K, K )
+ D21 = A( K, K-1 )
+ D12 = DCONJG( D21 )
+ DO 50 J = 1, NRHS
+ T1 = B( K-1, J )
+ T2 = B( K, J )
+ B( K-1, J ) = D11*T1 + D12*T2
+ B( K, J ) = D21*T1 + D22*T2
+ 50 CONTINUE
+ END IF
+*
+* Multiply by P(K) * inv(L(K)) if K < N.
+*
+ IF( K.NE.N ) THEN
+*
+* Apply the transformation.
+*
+ CALL ZGERU( N-K, NRHS, CONE, A( K+1, K ), 1,
+ $ B( K, 1 ), LDB, B( K+1, 1 ), LDB )
+ CALL ZGERU( N-K, NRHS, CONE, A( K+1, K-1 ), 1,
+ $ B( K-1, 1 ), LDB, B( K+1, 1 ), LDB )
+*
+* Interchange if a permutation was applied at the
+* K-th step of the factorization.
+*
+*
+* Swap the second of pair with IMAXth
+*
+ KP = ABS( IPIV( K ) )
+ IF( KP.NE.K )
+ $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+*
+* NOW swap the first of pair with Pth
+*
+ KP = ABS( IPIV( K-1 ) )
+ IF( KP.NE.K-1 )
+ $ CALL ZSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ),
+ $ LDB )
+*
+ END IF
+ K = K - 2
+ END IF
+ GO TO 40
+ 60 CONTINUE
+ END IF
+*--------------------------------------------------
+*
+* Compute B := A^H * B (conjugate transpose)
+*
+*--------------------------------------------------
+ ELSE
+*
+* Form B := U^H*B
+* where U = P(m)*inv(U(m))* ... *P(1)*inv(U(1))
+* and U^H = inv(U^H(1))*P(1)* ... *inv(U^H(m))*P(m)
+*
+ IF( LSAME( UPLO, 'U' ) ) THEN
+*
+* Loop backward applying the transformations.
+*
+ K = N
+ 70 IF( K.LT.1 )
+ $ GO TO 90
+*
+* 1 x 1 pivot block.
+*
+ IF( IPIV( K ).GT.0 ) THEN
+ IF( K.GT.1 ) THEN
+*
+* Interchange if P(K) != I.
+*
+ KP = IPIV( K )
+ IF( KP.NE.K )
+ $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+*
+* Apply the transformation
+* y = y - B' DCONJG(x),
+* where x is a column of A and y is a row of B.
+*
+ CALL ZLACGV( NRHS, B( K, 1 ), LDB )
+ CALL ZGEMV( 'Conjugate', K-1, NRHS, CONE, B, LDB,
+ $ A( 1, K ), 1, CONE, B( K, 1 ), LDB )
+ CALL ZLACGV( NRHS, B( K, 1 ), LDB )
+ END IF
+ IF( NOUNIT )
+ $ CALL ZSCAL( NRHS, A( K, K ), B( K, 1 ), LDB )
+ K = K - 1
+*
+* 2 x 2 pivot block.
+*
+ ELSE
+ IF( K.GT.2 ) THEN
+*
+* Swap the second of pair with Pth
+*
+ KP = ABS( IPIV( K ) )
+ IF( KP.NE.K )
+ $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+*
+* Now swap the first of pair with IMAX(r)th
+*
+ KP = ABS( IPIV( K-1 ) )
+ IF( KP.NE.K-1 )
+ $ CALL ZSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ),
+ $ LDB )
+*
+* Apply the transformations
+* y = y - B' DCONJG(x),
+* where x is a block column of A and y is a block
+* row of B.
+*
+ CALL ZLACGV( NRHS, B( K, 1 ), LDB )
+ CALL ZGEMV( 'Conjugate', K-2, NRHS, CONE, B, LDB,
+ $ A( 1, K ), 1, CONE, B( K, 1 ), LDB )
+ CALL ZLACGV( NRHS, B( K, 1 ), LDB )
+*
+ CALL ZLACGV( NRHS, B( K-1, 1 ), LDB )
+ CALL ZGEMV( 'Conjugate', K-2, NRHS, CONE, B, LDB,
+ $ A( 1, K-1 ), 1, CONE, B( K-1, 1 ), LDB )
+ CALL ZLACGV( NRHS, B( K-1, 1 ), LDB )
+ END IF
+*
+* Multiply by the diagonal block if non-unit.
+*
+ IF( NOUNIT ) THEN
+ D11 = A( K-1, K-1 )
+ D22 = A( K, K )
+ D12 = A( K-1, K )
+ D21 = DCONJG( D12 )
+ DO 80 J = 1, NRHS
+ T1 = B( K-1, J )
+ T2 = B( K, J )
+ B( K-1, J ) = D11*T1 + D12*T2
+ B( K, J ) = D21*T1 + D22*T2
+ 80 CONTINUE
+ END IF
+ K = K - 2
+ END IF
+ GO TO 70
+ 90 CONTINUE
+*
+* Form B := L^H*B
+* where L = P(1)*inv(L(1))* ... *P(m)*inv(L(m))
+* and L^H = inv(L^H(m))*P(m)* ... *inv(L^H(1))*P(1)
+*
+ ELSE
+*
+* Loop forward applying the L-transformations.
+*
+ K = 1
+ 100 CONTINUE
+ IF( K.GT.N )
+ $ GO TO 120
+*
+* 1 x 1 pivot block
+*
+ IF( IPIV( K ).GT.0 ) THEN
+ IF( K.LT.N ) THEN
+*
+* Interchange if P(K) != I.
+*
+ KP = IPIV( K )
+ IF( KP.NE.K )
+ $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+*
+* Apply the transformation
+*
+ CALL ZLACGV( NRHS, B( K, 1 ), LDB )
+ CALL ZGEMV( 'Conjugate', N-K, NRHS, CONE, B( K+1, 1 ),
+ $ LDB, A( K+1, K ), 1, CONE, B( K, 1 ), LDB )
+ CALL ZLACGV( NRHS, B( K, 1 ), LDB )
+ END IF
+ IF( NOUNIT )
+ $ CALL ZSCAL( NRHS, A( K, K ), B( K, 1 ), LDB )
+ K = K + 1
+*
+* 2 x 2 pivot block.
+*
+ ELSE
+ IF( K.LT.N-1 ) THEN
+*
+* Swap the first of pair with Pth
+*
+ KP = ABS( IPIV( K ) )
+ IF( KP.NE.K )
+ $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+*
+* Now swap the second of pair with IMAX(r)th
+*
+ KP = ABS( IPIV( K+1 ) )
+ IF( KP.NE.K+1 )
+ $ CALL ZSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ),
+ $ LDB )
+*
+* Apply the transformation
+*
+ CALL ZLACGV( NRHS, B( K+1, 1 ), LDB )
+ CALL ZGEMV( 'Conjugate', N-K-1, NRHS, CONE,
+ $ B( K+2, 1 ), LDB, A( K+2, K+1 ), 1, CONE,
+ $ B( K+1, 1 ), LDB )
+ CALL ZLACGV( NRHS, B( K+1, 1 ), LDB )
+*
+ CALL ZLACGV( NRHS, B( K, 1 ), LDB )
+ CALL ZGEMV( 'Conjugate', N-K-1, NRHS, CONE,
+ $ B( K+2, 1 ), LDB, A( K+2, K ), 1, CONE,
+ $ B( K, 1 ), LDB )
+ CALL ZLACGV( NRHS, B( K, 1 ), LDB )
+ END IF
+*
+* Multiply by the diagonal block if non-unit.
+*
+ IF( NOUNIT ) THEN
+ D11 = A( K, K )
+ D22 = A( K+1, K+1 )
+ D21 = A( K+1, K )
+ D12 = DCONJG( D21 )
+ DO 110 J = 1, NRHS
+ T1 = B( K, J )
+ T2 = B( K+1, J )
+ B( K, J ) = D11*T1 + D12*T2
+ B( K+1, J ) = D21*T1 + D22*T2
+ 110 CONTINUE
+ END IF
+ K = K + 2
+ END IF
+ GO TO 100
+ 120 CONTINUE
+ END IF
+*
+ END IF
+ RETURN
+*
+* End of ZLAVHE_ROOK
+*
+ END
projects / platform / upstream / lapack.git / commitdiff