SRC/dla_syrcond.f

   1 *> \brief \b DLA_SYRCOND estimates the Skeel condition number for a symmetric indefinite matrix.
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download DLA_SYRCOND + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dla_syrcond.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dla_syrcond.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dla_syrcond.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       DOUBLE PRECISION FUNCTION DLA_SYRCOND( UPLO, N, A, LDA, AF, LDAF,
  22 *                                              IPIV, CMODE, C, INFO, WORK,
  23 *                                              IWORK )
  24 *
  25 *       .. Scalar Arguments ..
  26 *       CHARACTER          UPLO
  27 *       INTEGER            N, LDA, LDAF, INFO, CMODE
  28 *       ..
  29 *       .. Array Arguments
  30 *       INTEGER            IWORK( * ), IPIV( * )
  31 *       DOUBLE PRECISION   A( LDA, * ), AF( LDAF, * ), WORK( * ), C( * )
  32 *       ..
  33 *
  34 *
  35 *> \par Purpose:
  36 *  =============
  37 *>
  38 *> \verbatim
  39 *>
  40 *>    DLA_SYRCOND estimates the Skeel condition number of  op(A) * op2(C)
  41 *>    where op2 is determined by CMODE as follows
  42 *>    CMODE =  1    op2(C) = C
  43 *>    CMODE =  0    op2(C) = I
  44 *>    CMODE = -1    op2(C) = inv(C)
  45 *>    The Skeel condition number cond(A) = norminf( |inv(A)||A| )
  46 *>    is computed by computing scaling factors R such that
  47 *>    diag(R)*A*op2(C) is row equilibrated and computing the standard
  48 *>    infinity-norm condition number.
  49 *> \endverbatim
  50 *
  51 *  Arguments:
  52 *  ==========
  53 *
  54 *> \param[in] UPLO
  55 *> \verbatim
  56 *>          UPLO is CHARACTER*1
  57 *>       = 'U':  Upper triangle of A is stored;
  58 *>       = 'L':  Lower triangle of A is stored.
  59 *> \endverbatim
  60 *>
  61 *> \param[in] N
  62 *> \verbatim
  63 *>          N is INTEGER
  64 *>     The number of linear equations, i.e., the order of the
  65 *>     matrix A.  N >= 0.
  66 *> \endverbatim
  67 *>
  68 *> \param[in] A
  69 *> \verbatim
  70 *>          A is DOUBLE PRECISION array, dimension (LDA,N)
  71 *>     On entry, the N-by-N matrix A.
  72 *> \endverbatim
  73 *>
  74 *> \param[in] LDA
  75 *> \verbatim
  76 *>          LDA is INTEGER
  77 *>     The leading dimension of the array A.  LDA >= max(1,N).
  78 *> \endverbatim
  79 *>
  80 *> \param[in] AF
  81 *> \verbatim
  82 *>          AF is DOUBLE PRECISION array, dimension (LDAF,N)
  83 *>     The block diagonal matrix D and the multipliers used to
  84 *>     obtain the factor U or L as computed by DSYTRF.
  85 *> \endverbatim
  86 *>
  87 *> \param[in] LDAF
  88 *> \verbatim
  89 *>          LDAF is INTEGER
  90 *>     The leading dimension of the array AF.  LDAF >= max(1,N).
  91 *> \endverbatim
  92 *>
  93 *> \param[in] IPIV
  94 *> \verbatim
  95 *>          IPIV is INTEGER array, dimension (N)
  96 *>     Details of the interchanges and the block structure of D
  97 *>     as determined by DSYTRF.
  98 *> \endverbatim
  99 *>
 100 *> \param[in] CMODE
 101 *> \verbatim
 102 *>          CMODE is INTEGER
 103 *>     Determines op2(C) in the formula op(A) * op2(C) as follows:
 104 *>     CMODE =  1    op2(C) = C
 105 *>     CMODE =  0    op2(C) = I
 106 *>     CMODE = -1    op2(C) = inv(C)
 107 *> \endverbatim
 108 *>
 109 *> \param[in] C
 110 *> \verbatim
 111 *>          C is DOUBLE PRECISION array, dimension (N)
 112 *>     The vector C in the formula op(A) * op2(C).
 113 *> \endverbatim
 114 *>
 115 *> \param[out] INFO
 116 *> \verbatim
 117 *>          INFO is INTEGER
 118 *>       = 0:  Successful exit.
 119 *>     i > 0:  The ith argument is invalid.
 120 *> \endverbatim
 121 *>
 122 *> \param[in] WORK
 123 *> \verbatim
 124 *>          WORK is DOUBLE PRECISION array, dimension (3*N).
 125 *>     Workspace.
 126 *> \endverbatim
 127 *>
 128 *> \param[in] IWORK
 129 *> \verbatim
 130 *>          IWORK is INTEGER array, dimension (N).
 131 *>     Workspace.
 132 *> \endverbatim
 133 *
 134 *  Authors:
 135 *  ========
 136 *
 137 *> \author Univ. of Tennessee
 138 *> \author Univ. of California Berkeley
 139 *> \author Univ. of Colorado Denver
 140 *> \author NAG Ltd.
 141 *
 142 *> \date September 2012
 143 *
 144 *> \ingroup doubleSYcomputational
 145 *
 146 *  =====================================================================
 147       DOUBLE PRECISION FUNCTION DLA_SYRCOND( UPLO, N, A, LDA, AF, LDAF,
 148      $                                       IPIV, CMODE, C, INFO, WORK,
 149      $                                       IWORK )
 150 *
 151 *  -- LAPACK computational routine (version 3.4.2) --
 152 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 153 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 154 *     September 2012
 155 *
 156 *     .. Scalar Arguments ..
 157       CHARACTER          UPLO
 158       INTEGER            N, LDA, LDAF, INFO, CMODE
 159 *     ..
 160 *     .. Array Arguments
 161       INTEGER            IWORK( * ), IPIV( * )
 162       DOUBLE PRECISION   A( LDA, * ), AF( LDAF, * ), WORK( * ), C( * )
 163 *     ..
 164 *
 165 *  =====================================================================
 166 *
 167 *     .. Local Scalars ..
 168       CHARACTER          NORMIN
 169       INTEGER            KASE, I, J
 170       DOUBLE PRECISION   AINVNM, SMLNUM, TMP
 171       LOGICAL            UP
 172 *     ..
 173 *     .. Local Arrays ..
 174       INTEGER            ISAVE( 3 )
 175 *     ..
 176 *     .. External Functions ..
 177       LOGICAL            LSAME
 178       INTEGER            IDAMAX
 179       DOUBLE PRECISION   DLAMCH
 180       EXTERNAL           LSAME, IDAMAX, DLAMCH
 181 *     ..
 182 *     .. External Subroutines ..
 183       EXTERNAL           DLACN2, DLATRS, DRSCL, XERBLA, DSYTRS
 184 *     ..
 185 *     .. Intrinsic Functions ..
 186       INTRINSIC          ABS, MAX
 187 *     ..
 188 *     .. Executable Statements ..
 189 *
 190       DLA_SYRCOND = 0.0D+0
 191 *
 192       INFO = 0
 193       IF( N.LT.0 ) THEN
 194          INFO = -2
 195       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
 196          INFO = -4
 197       ELSE IF( LDAF.LT.MAX( 1, N ) ) THEN
 198          INFO = -6
 199       END IF
 200       IF( INFO.NE.0 ) THEN
 201          CALL XERBLA( 'DLA_SYRCOND', -INFO )
 202          RETURN
 203       END IF
 204       IF( N.EQ.0 ) THEN
 205          DLA_SYRCOND = 1.0D+0
 206          RETURN
 207       END IF
 208       UP = .FALSE.
 209       IF ( LSAME( UPLO, 'U' ) ) UP = .TRUE.
 210 *
 211 *     Compute the equilibration matrix R such that
 212 *     inv(R)*A*C has unit 1-norm.
 213 *
 214       IF ( UP ) THEN
 215          DO I = 1, N
 216             TMP = 0.0D+0
 217             IF ( CMODE .EQ. 1 ) THEN
 218                DO J = 1, I
 219                   TMP = TMP + ABS( A( J, I ) * C( J ) )
 220                END DO
 221                DO J = I+1, N
 222                   TMP = TMP + ABS( A( I, J ) * C( J ) )
 223                END DO
 224             ELSE IF ( CMODE .EQ. 0 ) THEN
 225                DO J = 1, I
 226                   TMP = TMP + ABS( A( J, I ) )
 227                END DO
 228                DO J = I+1, N
 229                   TMP = TMP + ABS( A( I, J ) )
 230                END DO
 231             ELSE
 232                DO J = 1, I
 233                   TMP = TMP + ABS( A( J, I ) / C( J ) )
 234                END DO
 235                DO J = I+1, N
 236                   TMP = TMP + ABS( A( I, J ) / C( J ) )
 237                END DO
 238             END IF
 239             WORK( 2*N+I ) = TMP
 240          END DO
 241       ELSE
 242          DO I = 1, N
 243             TMP = 0.0D+0
 244             IF ( CMODE .EQ. 1 ) THEN
 245                DO J = 1, I
 246                   TMP = TMP + ABS( A( I, J ) * C( J ) )
 247                END DO
 248                DO J = I+1, N
 249                   TMP = TMP + ABS( A( J, I ) * C( J ) )
 250                END DO
 251             ELSE IF ( CMODE .EQ. 0 ) THEN
 252                DO J = 1, I
 253                   TMP = TMP + ABS( A( I, J ) )
 254                END DO
 255                DO J = I+1, N
 256                   TMP = TMP + ABS( A( J, I ) )
 257                END DO
 258             ELSE
 259                DO J = 1, I
 260                   TMP = TMP + ABS( A( I, J) / C( J ) )
 261                END DO
 262                DO J = I+1, N
 263                   TMP = TMP + ABS( A( J, I) / C( J ) )
 264                END DO
 265             END IF
 266             WORK( 2*N+I ) = TMP
 267          END DO
 268       ENDIF
 269 *
 270 *     Estimate the norm of inv(op(A)).
 271 *
 272       SMLNUM = DLAMCH( 'Safe minimum' )
 273       AINVNM = 0.0D+0
 274       NORMIN = 'N'
 275
 276       KASE = 0
 277    10 CONTINUE
 278       CALL DLACN2( N, WORK( N+1 ), WORK, IWORK, AINVNM, KASE, ISAVE )
 279       IF( KASE.NE.0 ) THEN
 280          IF( KASE.EQ.2 ) THEN
 281 *
 282 *           Multiply by R.
 283 *
 284             DO I = 1, N
 285                WORK( I ) = WORK( I ) * WORK( 2*N+I )
 286             END DO
 287
 288             IF ( UP ) THEN
 289                CALL DSYTRS( 'U', N, 1, AF, LDAF, IPIV, WORK, N, INFO )
 290             ELSE
 291                CALL DSYTRS( 'L', N, 1, AF, LDAF, IPIV, WORK, N, INFO )
 292             ENDIF
 293 *
 294 *           Multiply by inv(C).
 295 *
 296             IF ( CMODE .EQ. 1 ) THEN
 297                DO I = 1, N
 298                   WORK( I ) = WORK( I ) / C( I )
 299                END DO
 300             ELSE IF ( CMODE .EQ. -1 ) THEN
 301                DO I = 1, N
 302                   WORK( I ) = WORK( I ) * C( I )
 303                END DO
 304             END IF
 305          ELSE
 306 *
 307 *           Multiply by inv(C**T).
 308 *
 309             IF ( CMODE .EQ. 1 ) THEN
 310                DO I = 1, N
 311                   WORK( I ) = WORK( I ) / C( I )
 312                END DO
 313             ELSE IF ( CMODE .EQ. -1 ) THEN
 314                DO I = 1, N
 315                   WORK( I ) = WORK( I ) * C( I )
 316                END DO
 317             END IF
 318
 319             IF ( UP ) THEN
 320                CALL DSYTRS( 'U', N, 1, AF, LDAF, IPIV, WORK, N, INFO )
 321             ELSE
 322                CALL DSYTRS( 'L', N, 1, AF, LDAF, IPIV, WORK, N, INFO )
 323             ENDIF
 324 *
 325 *           Multiply by R.
 326 *
 327             DO I = 1, N
 328                WORK( I ) = WORK( I ) * WORK( 2*N+I )
 329             END DO
 330          END IF
 331 *
 332          GO TO 10
 333       END IF
 334 *
 335 *     Compute the estimate of the reciprocal condition number.
 336 *
 337       IF( AINVNM .NE. 0.0D+0 )
 338      $   DLA_SYRCOND = ( 1.0D+0 / AINVNM )
 339 *
 340       RETURN
 341 *
 342       END