LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
dla_porcond.f
Go to the documentation of this file.
00001 *> \brief \b DLA_PORCOND
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download DLA_PORCOND + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dla_porcond.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dla_porcond.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dla_porcond.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       DOUBLE PRECISION FUNCTION DLA_PORCOND( UPLO, N, A, LDA, AF, LDAF,
00022 *                                              CMODE, C, INFO, WORK,
00023 *                                              IWORK )
00024 * 
00025 *       .. Scalar Arguments ..
00026 *       CHARACTER          UPLO
00027 *       INTEGER            N, LDA, LDAF, INFO, CMODE
00028 *       DOUBLE PRECISION   A( LDA, * ), AF( LDAF, * ), WORK( * ),
00029 *      $                   C( * )
00030 *       ..
00031 *       .. Array Arguments ..
00032 *       INTEGER            IWORK( * )
00033 *       ..
00034 *  
00035 *
00036 *> \par Purpose:
00037 *  =============
00038 *>
00039 *> \verbatim
00040 *>
00041 *>    DLA_PORCOND Estimates the Skeel condition number of  op(A) * op2(C)
00042 *>    where op2 is determined by CMODE as follows
00043 *>    CMODE =  1    op2(C) = C
00044 *>    CMODE =  0    op2(C) = I
00045 *>    CMODE = -1    op2(C) = inv(C)
00046 *>    The Skeel condition number  cond(A) = norminf( |inv(A)||A| )
00047 *>    is computed by computing scaling factors R such that
00048 *>    diag(R)*A*op2(C) is row equilibrated and computing the standard
00049 *>    infinity-norm condition number.
00050 *> \endverbatim
00051 *
00052 *  Arguments:
00053 *  ==========
00054 *
00055 *> \param[in] UPLO
00056 *> \verbatim
00057 *>          UPLO is CHARACTER*1
00058 *>       = 'U':  Upper triangle of A is stored;
00059 *>       = 'L':  Lower triangle of A is stored.
00060 *> \endverbatim
00061 *>
00062 *> \param[in] N
00063 *> \verbatim
00064 *>          N is INTEGER
00065 *>     The number of linear equations, i.e., the order of the
00066 *>     matrix A.  N >= 0.
00067 *> \endverbatim
00068 *>
00069 *> \param[in] A
00070 *> \verbatim
00071 *>          A is DOUBLE PRECISION array, dimension (LDA,N)
00072 *>     On entry, the N-by-N matrix A.
00073 *> \endverbatim
00074 *>
00075 *> \param[in] LDA
00076 *> \verbatim
00077 *>          LDA is INTEGER
00078 *>     The leading dimension of the array A.  LDA >= max(1,N).
00079 *> \endverbatim
00080 *>
00081 *> \param[in] AF
00082 *> \verbatim
00083 *>          AF is DOUBLE PRECISION array, dimension (LDAF,N)
00084 *>     The triangular factor U or L from the Cholesky factorization
00085 *>     A = U**T*U or A = L*L**T, as computed by DPOTRF.
00086 *> \endverbatim
00087 *>
00088 *> \param[in] LDAF
00089 *> \verbatim
00090 *>          LDAF is INTEGER
00091 *>     The leading dimension of the array AF.  LDAF >= max(1,N).
00092 *> \endverbatim
00093 *>
00094 *> \param[in] CMODE
00095 *> \verbatim
00096 *>          CMODE is INTEGER
00097 *>     Determines op2(C) in the formula op(A) * op2(C) as follows:
00098 *>     CMODE =  1    op2(C) = C
00099 *>     CMODE =  0    op2(C) = I
00100 *>     CMODE = -1    op2(C) = inv(C)
00101 *> \endverbatim
00102 *>
00103 *> \param[in] C
00104 *> \verbatim
00105 *>          C is DOUBLE PRECISION array, dimension (N)
00106 *>     The vector C in the formula op(A) * op2(C).
00107 *> \endverbatim
00108 *>
00109 *> \param[out] INFO
00110 *> \verbatim
00111 *>          INFO is INTEGER
00112 *>       = 0:  Successful exit.
00113 *>     i > 0:  The ith argument is invalid.
00114 *> \endverbatim
00115 *>
00116 *> \param[in] WORK
00117 *> \verbatim
00118 *>          WORK is DOUBLE PRECISION array, dimension (3*N).
00119 *>     Workspace.
00120 *> \endverbatim
00121 *>
00122 *> \param[in] IWORK
00123 *> \verbatim
00124 *>          IWORK is INTEGER array, dimension (N).
00125 *>     Workspace.
00126 *> \endverbatim
00127 *
00128 *  Authors:
00129 *  ========
00130 *
00131 *> \author Univ. of Tennessee 
00132 *> \author Univ. of California Berkeley 
00133 *> \author Univ. of Colorado Denver 
00134 *> \author NAG Ltd. 
00135 *
00136 *> \date November 2011
00137 *
00138 *> \ingroup doublePOcomputational
00139 *
00140 *  =====================================================================
00141       DOUBLE PRECISION FUNCTION DLA_PORCOND( UPLO, N, A, LDA, AF, LDAF,
00142      $                                       CMODE, C, INFO, WORK,
00143      $                                       IWORK )
00144 *
00145 *  -- LAPACK computational routine (version 3.4.0) --
00146 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00147 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00148 *     November 2011
00149 *
00150 *     .. Scalar Arguments ..
00151       CHARACTER          UPLO
00152       INTEGER            N, LDA, LDAF, INFO, CMODE
00153       DOUBLE PRECISION   A( LDA, * ), AF( LDAF, * ), WORK( * ),
00154      $                   C( * )
00155 *     ..
00156 *     .. Array Arguments ..
00157       INTEGER            IWORK( * )
00158 *     ..
00159 *
00160 *  =====================================================================
00161 *
00162 *     .. Local Scalars ..
00163       INTEGER            KASE, I, J
00164       DOUBLE PRECISION   AINVNM, TMP
00165       LOGICAL            UP
00166 *     ..
00167 *     .. Array Arguments ..
00168       INTEGER            ISAVE( 3 )
00169 *     ..
00170 *     .. External Functions ..
00171       LOGICAL            LSAME
00172       INTEGER            IDAMAX
00173       EXTERNAL           LSAME, IDAMAX
00174 *     ..
00175 *     .. External Subroutines ..
00176       EXTERNAL           DLACN2, DPOTRS, XERBLA
00177 *     ..
00178 *     .. Intrinsic Functions ..
00179       INTRINSIC          ABS, MAX
00180 *     ..
00181 *     .. Executable Statements ..
00182 *
00183       DLA_PORCOND = 0.0D+0
00184 *
00185       INFO = 0
00186       IF( N.LT.0 ) THEN
00187          INFO = -2
00188       END IF
00189       IF( INFO.NE.0 ) THEN
00190          CALL XERBLA( 'DLA_PORCOND', -INFO )
00191          RETURN
00192       END IF
00193 
00194       IF( N.EQ.0 ) THEN
00195          DLA_PORCOND = 1.0D+0
00196          RETURN
00197       END IF
00198       UP = .FALSE.
00199       IF ( LSAME( UPLO, 'U' ) ) UP = .TRUE.
00200 *
00201 *     Compute the equilibration matrix R such that
00202 *     inv(R)*A*C has unit 1-norm.
00203 *
00204       IF ( UP ) THEN
00205          DO I = 1, N
00206             TMP = 0.0D+0
00207             IF ( CMODE .EQ. 1 ) THEN
00208                DO J = 1, I
00209                   TMP = TMP + ABS( A( J, I ) * C( J ) )
00210                END DO
00211                DO J = I+1, N
00212                   TMP = TMP + ABS( A( I, J ) * C( J ) )
00213                END DO
00214             ELSE IF ( CMODE .EQ. 0 ) THEN
00215                DO J = 1, I
00216                   TMP = TMP + ABS( A( J, I ) )
00217                END DO
00218                DO J = I+1, N
00219                   TMP = TMP + ABS( A( I, J ) )
00220                END DO
00221             ELSE
00222                DO J = 1, I
00223                   TMP = TMP + ABS( A( J ,I ) / C( J ) )
00224                END DO
00225                DO J = I+1, N
00226                   TMP = TMP + ABS( A( I, J ) / C( J ) )
00227                END DO
00228             END IF
00229             WORK( 2*N+I ) = TMP
00230          END DO
00231       ELSE
00232          DO I = 1, N
00233             TMP = 0.0D+0
00234             IF ( CMODE .EQ. 1 ) THEN
00235                DO J = 1, I
00236                   TMP = TMP + ABS( A( I, J ) * C( J ) )
00237                END DO
00238                DO J = I+1, N
00239                   TMP = TMP + ABS( A( J, I ) * C( J ) )
00240                END DO
00241             ELSE IF ( CMODE .EQ. 0 ) THEN
00242                DO J = 1, I
00243                   TMP = TMP + ABS( A( I, J ) )
00244                END DO
00245                DO J = I+1, N
00246                   TMP = TMP + ABS( A( J, I ) )
00247                END DO
00248             ELSE
00249                DO J = 1, I
00250                   TMP = TMP + ABS( A( I, J ) / C( J ) )
00251                END DO
00252                DO J = I+1, N
00253                   TMP = TMP + ABS( A( J, I ) / C( J ) )
00254                END DO
00255             END IF
00256             WORK( 2*N+I ) = TMP
00257          END DO
00258       ENDIF
00259 *
00260 *     Estimate the norm of inv(op(A)).
00261 *
00262       AINVNM = 0.0D+0
00263 
00264       KASE = 0
00265    10 CONTINUE
00266       CALL DLACN2( N, WORK( N+1 ), WORK, IWORK, AINVNM, KASE, ISAVE )
00267       IF( KASE.NE.0 ) THEN
00268          IF( KASE.EQ.2 ) THEN
00269 *
00270 *           Multiply by R.
00271 *
00272             DO I = 1, N
00273                WORK( I ) = WORK( I ) * WORK( 2*N+I )
00274             END DO
00275 
00276             IF (UP) THEN
00277                CALL DPOTRS( 'Upper', N, 1, AF, LDAF, WORK, N, INFO )
00278             ELSE
00279                CALL DPOTRS( 'Lower', N, 1, AF, LDAF, WORK, N, INFO )
00280             ENDIF
00281 *
00282 *           Multiply by inv(C).
00283 *
00284             IF ( CMODE .EQ. 1 ) THEN
00285                DO I = 1, N
00286                   WORK( I ) = WORK( I ) / C( I )
00287                END DO
00288             ELSE IF ( CMODE .EQ. -1 ) THEN
00289                DO I = 1, N
00290                   WORK( I ) = WORK( I ) * C( I )
00291                END DO
00292             END IF
00293          ELSE
00294 *
00295 *           Multiply by inv(C**T).
00296 *
00297             IF ( CMODE .EQ. 1 ) THEN
00298                DO I = 1, N
00299                   WORK( I ) = WORK( I ) / C( I )
00300                END DO
00301             ELSE IF ( CMODE .EQ. -1 ) THEN
00302                DO I = 1, N
00303                   WORK( I ) = WORK( I ) * C( I )
00304                END DO
00305             END IF
00306 
00307             IF ( UP ) THEN
00308                CALL DPOTRS( 'Upper', N, 1, AF, LDAF, WORK, N, INFO )
00309             ELSE
00310                CALL DPOTRS( 'Lower', N, 1, AF, LDAF, WORK, N, INFO )
00311             ENDIF
00312 *
00313 *           Multiply by R.
00314 *
00315             DO I = 1, N
00316                WORK( I ) = WORK( I ) * WORK( 2*N+I )
00317             END DO
00318          END IF
00319          GO TO 10
00320       END IF
00321 *
00322 *     Compute the estimate of the reciprocal condition number.
00323 *
00324       IF( AINVNM .NE. 0.0D+0 )
00325      $   DLA_PORCOND = ( 1.0D+0 / AINVNM )
00326 *
00327       RETURN
00328 *
00329       END
 All Files Functions