Julie
bc6a755be4
This is really old school, but a lot of times we have users sending us copy pasting of codes, and that is the only way to know the version of the code.
256 lines
7.0 KiB
FortranFixed
256 lines
7.0 KiB
FortranFixed
*> \brief \b CPPCON
|
|
*
|
|
* =========== DOCUMENTATION ===========
|
|
*
|
|
* Online html documentation available at
|
|
* http://www.netlib.org/lapack/explore-html/
|
|
*
|
|
*> \htmlonly
|
|
*> Download CPPCON + dependencies
|
|
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cppcon.f">
|
|
*> [TGZ]</a>
|
|
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cppcon.f">
|
|
*> [ZIP]</a>
|
|
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cppcon.f">
|
|
*> [TXT]</a>
|
|
*> \endhtmlonly
|
|
*
|
|
* Definition:
|
|
* ===========
|
|
*
|
|
* SUBROUTINE CPPCON( UPLO, N, AP, ANORM, RCOND, WORK, RWORK, INFO )
|
|
*
|
|
* .. Scalar Arguments ..
|
|
* CHARACTER UPLO
|
|
* INTEGER INFO, N
|
|
* REAL ANORM, RCOND
|
|
* ..
|
|
* .. Array Arguments ..
|
|
* REAL RWORK( * )
|
|
* COMPLEX AP( * ), WORK( * )
|
|
* ..
|
|
*
|
|
*
|
|
*> \par Purpose:
|
|
* =============
|
|
*>
|
|
*> \verbatim
|
|
*>
|
|
*> CPPCON estimates the reciprocal of the condition number (in the
|
|
*> 1-norm) of a complex Hermitian positive definite packed matrix using
|
|
*> the Cholesky factorization A = U**H*U or A = L*L**H computed by
|
|
*> CPPTRF.
|
|
*>
|
|
*> 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
|
|
*> = '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] AP
|
|
*> \verbatim
|
|
*> AP is COMPLEX array, dimension (N*(N+1)/2)
|
|
*> The triangular factor U or L from the Cholesky factorization
|
|
*> A = U**H*U or A = L*L**H, packed columnwise in a linear
|
|
*> array. The j-th column of U or L is stored in the array AP
|
|
*> as follows:
|
|
*> if UPLO = 'U', AP(i + (j-1)*j/2) = U(i,j) for 1<=i<=j;
|
|
*> if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = L(i,j) for j<=i<=n.
|
|
*> \endverbatim
|
|
*>
|
|
*> \param[in] ANORM
|
|
*> \verbatim
|
|
*> ANORM is REAL
|
|
*> The 1-norm (or infinity-norm) of the Hermitian 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] RWORK
|
|
*> \verbatim
|
|
*> RWORK is REAL 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 December 2016
|
|
*
|
|
*> \ingroup complexOTHERcomputational
|
|
*
|
|
* =====================================================================
|
|
SUBROUTINE CPPCON( UPLO, N, AP, ANORM, RCOND, WORK, RWORK, INFO )
|
|
*
|
|
* -- LAPACK computational routine (version 3.7.0) --
|
|
* -- LAPACK is a software package provided by Univ. of Tennessee, --
|
|
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
|
|
* December 2016
|
|
*
|
|
* .. Scalar Arguments ..
|
|
CHARACTER UPLO
|
|
INTEGER INFO, N
|
|
REAL ANORM, RCOND
|
|
* ..
|
|
* .. Array Arguments ..
|
|
REAL RWORK( * )
|
|
COMPLEX AP( * ), WORK( * )
|
|
* ..
|
|
*
|
|
* =====================================================================
|
|
*
|
|
* .. Parameters ..
|
|
REAL ONE, ZERO
|
|
PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 )
|
|
* ..
|
|
* .. Local Scalars ..
|
|
LOGICAL UPPER
|
|
CHARACTER NORMIN
|
|
INTEGER IX, KASE
|
|
REAL AINVNM, SCALE, SCALEL, SCALEU, SMLNUM
|
|
COMPLEX ZDUM
|
|
* ..
|
|
* .. Local Arrays ..
|
|
INTEGER ISAVE( 3 )
|
|
* ..
|
|
* .. External Functions ..
|
|
LOGICAL LSAME
|
|
INTEGER ICAMAX
|
|
REAL SLAMCH
|
|
EXTERNAL LSAME, ICAMAX, SLAMCH
|
|
* ..
|
|
* .. External Subroutines ..
|
|
EXTERNAL CLACN2, CLATPS, CSRSCL, XERBLA
|
|
* ..
|
|
* .. Intrinsic Functions ..
|
|
INTRINSIC ABS, AIMAG, REAL
|
|
* ..
|
|
* .. Statement Functions ..
|
|
REAL CABS1
|
|
* ..
|
|
* .. Statement Function definitions ..
|
|
CABS1( ZDUM ) = ABS( REAL( ZDUM ) ) + ABS( AIMAG( ZDUM ) )
|
|
* ..
|
|
* .. 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( ANORM.LT.ZERO ) THEN
|
|
INFO = -4
|
|
END IF
|
|
IF( INFO.NE.0 ) THEN
|
|
CALL XERBLA( 'CPPCON', -INFO )
|
|
RETURN
|
|
END IF
|
|
*
|
|
* Quick return if possible
|
|
*
|
|
RCOND = ZERO
|
|
IF( N.EQ.0 ) THEN
|
|
RCOND = ONE
|
|
RETURN
|
|
ELSE IF( ANORM.EQ.ZERO ) THEN
|
|
RETURN
|
|
END IF
|
|
*
|
|
SMLNUM = SLAMCH( 'Safe minimum' )
|
|
*
|
|
* Estimate the 1-norm of the inverse.
|
|
*
|
|
KASE = 0
|
|
NORMIN = 'N'
|
|
10 CONTINUE
|
|
CALL CLACN2( N, WORK( N+1 ), WORK, AINVNM, KASE, ISAVE )
|
|
IF( KASE.NE.0 ) THEN
|
|
IF( UPPER ) THEN
|
|
*
|
|
* Multiply by inv(U**H).
|
|
*
|
|
CALL CLATPS( 'Upper', 'Conjugate transpose', 'Non-unit',
|
|
$ NORMIN, N, AP, WORK, SCALEL, RWORK, INFO )
|
|
NORMIN = 'Y'
|
|
*
|
|
* Multiply by inv(U).
|
|
*
|
|
CALL CLATPS( 'Upper', 'No transpose', 'Non-unit', NORMIN, N,
|
|
$ AP, WORK, SCALEU, RWORK, INFO )
|
|
ELSE
|
|
*
|
|
* Multiply by inv(L).
|
|
*
|
|
CALL CLATPS( 'Lower', 'No transpose', 'Non-unit', NORMIN, N,
|
|
$ AP, WORK, SCALEL, RWORK, INFO )
|
|
NORMIN = 'Y'
|
|
*
|
|
* Multiply by inv(L**H).
|
|
*
|
|
CALL CLATPS( 'Lower', 'Conjugate transpose', 'Non-unit',
|
|
$ NORMIN, N, AP, WORK, SCALEU, RWORK, INFO )
|
|
END IF
|
|
*
|
|
* Multiply by 1/SCALE if doing so will not cause overflow.
|
|
*
|
|
SCALE = SCALEL*SCALEU
|
|
IF( SCALE.NE.ONE ) THEN
|
|
IX = ICAMAX( N, WORK, 1 )
|
|
IF( SCALE.LT.CABS1( WORK( IX ) )*SMLNUM .OR. SCALE.EQ.ZERO )
|
|
$ GO TO 20
|
|
CALL CSRSCL( N, SCALE, WORK, 1 )
|
|
END IF
|
|
GO TO 10
|
|
END IF
|
|
*
|
|
* Compute the estimate of the reciprocal condition number.
|
|
*
|
|
IF( AINVNM.NE.ZERO )
|
|
$ RCOND = ( ONE / AINVNM ) / ANORM
|
|
*
|
|
20 CONTINUE
|
|
RETURN
|
|
*
|
|
* End of CPPCON
|
|
*
|
|
END
|