SRC/dla_porcond.f

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