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