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