LAPACK 3.3.1
Linear Algebra PACKage

cppcon.f

Go to the documentation of this file.
00001       SUBROUTINE CPPCON( UPLO, N, AP, ANORM, RCOND, WORK, RWORK, INFO )
00002 *
00003 *  -- LAPACK routine (version 3.3.1) --
00004 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00005 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00006 *  -- April 2011                                                      --
00007 *
00008 *     Modified to call CLACN2 in place of CLACON, 10 Feb 03, SJH.
00009 *
00010 *     .. Scalar Arguments ..
00011       CHARACTER          UPLO
00012       INTEGER            INFO, N
00013       REAL               ANORM, RCOND
00014 *     ..
00015 *     .. Array Arguments ..
00016       REAL               RWORK( * )
00017       COMPLEX            AP( * ), WORK( * )
00018 *     ..
00019 *
00020 *  Purpose
00021 *  =======
00022 *
00023 *  CPPCON estimates the reciprocal of the condition number (in the 
00024 *  1-norm) of a complex Hermitian positive definite packed matrix using
00025 *  the Cholesky factorization A = U**H*U or A = L*L**H computed by
00026 *  CPPTRF.
00027 *
00028 *  An estimate is obtained for norm(inv(A)), and the reciprocal of the
00029 *  condition number is computed as RCOND = 1 / (ANORM * norm(inv(A))).
00030 *
00031 *  Arguments
00032 *  =========
00033 *
00034 *  UPLO    (input) CHARACTER*1
00035 *          = 'U':  Upper triangle of A is stored;
00036 *          = 'L':  Lower triangle of A is stored.
00037 *
00038 *  N       (input) INTEGER
00039 *          The order of the matrix A.  N >= 0.
00040 *
00041 *  AP      (input) COMPLEX array, dimension (N*(N+1)/2)
00042 *          The triangular factor U or L from the Cholesky factorization
00043 *          A = U**H*U or A = L*L**H, packed columnwise in a linear
00044 *          array.  The j-th column of U or L is stored in the array AP
00045 *          as follows:
00046 *          if UPLO = 'U', AP(i + (j-1)*j/2) = U(i,j) for 1<=i<=j;
00047 *          if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = L(i,j) for j<=i<=n.
00048 *
00049 *  ANORM   (input) REAL
00050 *          The 1-norm (or infinity-norm) of the Hermitian matrix A.
00051 *
00052 *  RCOND   (output) REAL
00053 *          The reciprocal of the condition number of the matrix A,
00054 *          computed as RCOND = 1/(ANORM * AINVNM), where AINVNM is an
00055 *          estimate of the 1-norm of inv(A) computed in this routine.
00056 *
00057 *  WORK    (workspace) COMPLEX array, dimension (2*N)
00058 *
00059 *  RWORK   (workspace) REAL array, dimension (N)
00060 *
00061 *  INFO    (output) INTEGER
00062 *          = 0:  successful exit
00063 *          < 0:  if INFO = -i, the i-th argument had an illegal value
00064 *
00065 *  =====================================================================
00066 *
00067 *     .. Parameters ..
00068       REAL               ONE, ZERO
00069       PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
00070 *     ..
00071 *     .. Local Scalars ..
00072       LOGICAL            UPPER
00073       CHARACTER          NORMIN
00074       INTEGER            IX, KASE
00075       REAL               AINVNM, SCALE, SCALEL, SCALEU, SMLNUM
00076       COMPLEX            ZDUM
00077 *     ..
00078 *     .. Local Arrays ..
00079       INTEGER            ISAVE( 3 )
00080 *     ..
00081 *     .. External Functions ..
00082       LOGICAL            LSAME
00083       INTEGER            ICAMAX
00084       REAL               SLAMCH
00085       EXTERNAL           LSAME, ICAMAX, SLAMCH
00086 *     ..
00087 *     .. External Subroutines ..
00088       EXTERNAL           CLACN2, CLATPS, CSRSCL, XERBLA
00089 *     ..
00090 *     .. Intrinsic Functions ..
00091       INTRINSIC          ABS, AIMAG, REAL
00092 *     ..
00093 *     .. Statement Functions ..
00094       REAL               CABS1
00095 *     ..
00096 *     .. Statement Function definitions ..
00097       CABS1( ZDUM ) = ABS( REAL( ZDUM ) ) + ABS( AIMAG( ZDUM ) )
00098 *     ..
00099 *     .. Executable Statements ..
00100 *
00101 *     Test the input parameters.
00102 *
00103       INFO = 0
00104       UPPER = LSAME( UPLO, 'U' )
00105       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
00106          INFO = -1
00107       ELSE IF( N.LT.0 ) THEN
00108          INFO = -2
00109       ELSE IF( ANORM.LT.ZERO ) THEN
00110          INFO = -4
00111       END IF
00112       IF( INFO.NE.0 ) THEN
00113          CALL XERBLA( 'CPPCON', -INFO )
00114          RETURN
00115       END IF
00116 *
00117 *     Quick return if possible
00118 *
00119       RCOND = ZERO
00120       IF( N.EQ.0 ) THEN
00121          RCOND = ONE
00122          RETURN
00123       ELSE IF( ANORM.EQ.ZERO ) THEN
00124          RETURN
00125       END IF
00126 *
00127       SMLNUM = SLAMCH( 'Safe minimum' )
00128 *
00129 *     Estimate the 1-norm of the inverse.
00130 *
00131       KASE = 0
00132       NORMIN = 'N'
00133    10 CONTINUE
00134       CALL CLACN2( N, WORK( N+1 ), WORK, AINVNM, KASE, ISAVE )
00135       IF( KASE.NE.0 ) THEN
00136          IF( UPPER ) THEN
00137 *
00138 *           Multiply by inv(U**H).
00139 *
00140             CALL CLATPS( 'Upper', 'Conjugate transpose', 'Non-unit',
00141      $                   NORMIN, N, AP, WORK, SCALEL, RWORK, INFO )
00142             NORMIN = 'Y'
00143 *
00144 *           Multiply by inv(U).
00145 *
00146             CALL CLATPS( 'Upper', 'No transpose', 'Non-unit', NORMIN, N,
00147      $                   AP, WORK, SCALEU, RWORK, INFO )
00148          ELSE
00149 *
00150 *           Multiply by inv(L).
00151 *
00152             CALL CLATPS( 'Lower', 'No transpose', 'Non-unit', NORMIN, N,
00153      $                   AP, WORK, SCALEL, RWORK, INFO )
00154             NORMIN = 'Y'
00155 *
00156 *           Multiply by inv(L**H).
00157 *
00158             CALL CLATPS( 'Lower', 'Conjugate transpose', 'Non-unit',
00159      $                   NORMIN, N, AP, WORK, SCALEU, RWORK, INFO )
00160          END IF
00161 *
00162 *        Multiply by 1/SCALE if doing so will not cause overflow.
00163 *
00164          SCALE = SCALEL*SCALEU
00165          IF( SCALE.NE.ONE ) THEN
00166             IX = ICAMAX( N, WORK, 1 )
00167             IF( SCALE.LT.CABS1( WORK( IX ) )*SMLNUM .OR. SCALE.EQ.ZERO )
00168      $         GO TO 20
00169             CALL CSRSCL( N, SCALE, WORK, 1 )
00170          END IF
00171          GO TO 10
00172       END IF
00173 *
00174 *     Compute the estimate of the reciprocal condition number.
00175 *
00176       IF( AINVNM.NE.ZERO )
00177      $   RCOND = ( ONE / AINVNM ) / ANORM
00178 *
00179    20 CONTINUE
00180       RETURN
00181 *
00182 *     End of CPPCON
00183 *
00184       END
 All Files Functions