LAPACK: SRC/sppcon.f Source File

Go to the documentation of this file.
00001       SUBROUTINE SPPCON( UPLO, N, AP, ANORM, RCOND, WORK, IWORK, INFO )
00002 *
00003 *  -- LAPACK routine (version 3.2) --
00004 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00005 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00006 *     November 2006
00007 *
00008 *     Modified to call SLACN2 in place of SLACON, 7 Feb 03, SJH.
00009 *
00010 *     .. Scalar Arguments ..
00011       CHARACTER          UPLO
00012       INTEGER            INFO, N
00013       REAL               ANORM, RCOND
00014 *     ..
00015 *     .. Array Arguments ..
00016       INTEGER            IWORK( * )
00017       REAL               AP( * ), WORK( * )
00018 *     ..
00019 *
00020 *  Purpose
00021 *  =======
00022 *
00023 *  SPPCON estimates the reciprocal of the condition number (in the
00024 *  1-norm) of a real symmetric positive definite packed matrix using
00025 *  the Cholesky factorization A = U**T*U or A = L*L**T computed by
00026 *  SPPTRF.
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) REAL array, dimension (N*(N+1)/2)
00042 *          The triangular factor U or L from the Cholesky factorization
00043 *          A = U**T*U or A = L*L**T, 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 symmetric 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) REAL array, dimension (3*N)
00058 *
00059 *  IWORK   (workspace) INTEGER 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 *     ..
00077 *     .. Local Arrays ..
00078       INTEGER            ISAVE( 3 )
00079 *     ..
00080 *     .. External Functions ..
00081       LOGICAL            LSAME
00082       INTEGER            ISAMAX
00083       REAL               SLAMCH
00084       EXTERNAL           LSAME, ISAMAX, SLAMCH
00085 *     ..
00086 *     .. External Subroutines ..
00087       EXTERNAL           SLACN2, SLATPS, SRSCL, XERBLA
00088 *     ..
00089 *     .. Intrinsic Functions ..
00090       INTRINSIC          ABS
00091 *     ..
00092 *     .. Executable Statements ..
00093 *
00094 *     Test the input parameters.
00095 *
00096       INFO = 0
00097       UPPER = LSAME( UPLO, 'U' )
00098       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
00099          INFO = -1
00100       ELSE IF( N.LT.0 ) THEN
00101          INFO = -2
00102       ELSE IF( ANORM.LT.ZERO ) THEN
00103          INFO = -4
00104       END IF
00105       IF( INFO.NE.0 ) THEN
00106          CALL XERBLA( 'SPPCON', -INFO )
00107          RETURN
00108       END IF
00109 *
00110 *     Quick return if possible
00111 *
00112       RCOND = ZERO
00113       IF( N.EQ.0 ) THEN
00114          RCOND = ONE
00115          RETURN
00116       ELSE IF( ANORM.EQ.ZERO ) THEN
00117          RETURN
00118       END IF
00119 *
00120       SMLNUM = SLAMCH( 'Safe minimum' )
00121 *
00122 *     Estimate the 1-norm of the inverse.
00123 *
00124       KASE = 0
00125       NORMIN = 'N'
00126    10 CONTINUE
00127       CALL SLACN2( N, WORK( N+1 ), WORK, IWORK, AINVNM, KASE, ISAVE )
00128       IF( KASE.NE.0 ) THEN
00129          IF( UPPER ) THEN
00130 *
00131 *           Multiply by inv(U').
00132 *
00133             CALL SLATPS( 'Upper', 'Transpose', 'Non-unit', NORMIN, N,
00134      $                   AP, WORK, SCALEL, WORK( 2*N+1 ), INFO )
00135             NORMIN = 'Y'
00136 *
00137 *           Multiply by inv(U).
00138 *
00139             CALL SLATPS( 'Upper', 'No transpose', 'Non-unit', NORMIN, N,
00140      $                   AP, WORK, SCALEU, WORK( 2*N+1 ), INFO )
00141          ELSE
00142 *
00143 *           Multiply by inv(L).
00144 *
00145             CALL SLATPS( 'Lower', 'No transpose', 'Non-unit', NORMIN, N,
00146      $                   AP, WORK, SCALEL, WORK( 2*N+1 ), INFO )
00147             NORMIN = 'Y'
00148 *
00149 *           Multiply by inv(L').
00150 *
00151             CALL SLATPS( 'Lower', 'Transpose', 'Non-unit', NORMIN, N,
00152      $                   AP, WORK, SCALEU, WORK( 2*N+1 ), INFO )
00153          END IF
00154 *
00155 *        Multiply by 1/SCALE if doing so will not cause overflow.
00156 *
00157          SCALE = SCALEL*SCALEU
00158          IF( SCALE.NE.ONE ) THEN
00159             IX = ISAMAX( N, WORK, 1 )
00160             IF( SCALE.LT.ABS( WORK( IX ) )*SMLNUM .OR. SCALE.EQ.ZERO )
00161      $         GO TO 20
00162             CALL SRSCL( N, SCALE, WORK, 1 )
00163          END IF
00164          GO TO 10
00165       END IF
00166 *
00167 *     Compute the estimate of the reciprocal condition number.
00168 *
00169       IF( AINVNM.NE.ZERO )
00170      $   RCOND = ( ONE / AINVNM ) / ANORM
00171 *
00172    20 CONTINUE
00173       RETURN
00174 *
00175 *     End of SPPCON
00176 *
00177       END