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