LAPACK: SRC/dla_syrcond.f Source File

Go to the documentation of this file.
00001       DOUBLE PRECISION FUNCTION DLA_SYRCOND( UPLO, N, A, LDA, AF, LDAF, 
00002      $                                       IPIV, CMODE, C, INFO, WORK,
00003      $                                       IWORK )
00004 *
00005 *     -- LAPACK routine (version 3.2.1)                                 --
00006 *     -- Contributed by James Demmel, Deaglan Halligan, Yozo Hida and --
00007 *     -- Jason Riedy of Univ. of California Berkeley.                 --
00008 *     -- April 2009                                                   --
00009 *
00010 *     -- LAPACK is a software package provided by Univ. of Tennessee, --
00011 *     -- Univ. of California Berkeley and NAG Ltd.                    --
00012 *
00013       IMPLICIT NONE
00014 *     ..
00015 *     .. Scalar Arguments ..
00016       CHARACTER          UPLO
00017       INTEGER            N, LDA, LDAF, INFO, CMODE
00018 *     ..
00019 *     .. Array Arguments
00020       INTEGER            IWORK( * ), IPIV( * )
00021       DOUBLE PRECISION   A( LDA, * ), AF( LDAF, * ), WORK( * ), C( * )
00022 *     ..
00023 *
00024 *  Purpose
00025 *  =======
00026 *
00027 *     DLA_SYRCOND estimates the Skeel condition number of  op(A) * op2(C)
00028 *     where op2 is determined by CMODE as follows
00029 *     CMODE =  1    op2(C) = C
00030 *     CMODE =  0    op2(C) = I
00031 *     CMODE = -1    op2(C) = inv(C)
00032 *     The Skeel condition number cond(A) = norminf( |inv(A)||A| )
00033 *     is computed by computing scaling factors R such that
00034 *     diag(R)*A*op2(C) is row equilibrated and computing the standard
00035 *     infinity-norm condition number.
00036 *
00037 *  Arguments
00038 *  ==========
00039 *
00040 *     UPLO    (input) CHARACTER*1
00041 *       = 'U':  Upper triangle of A is stored;
00042 *       = 'L':  Lower triangle of A is stored.
00043 *
00044 *     N       (input) INTEGER
00045 *     The number of linear equations, i.e., the order of the
00046 *     matrix A.  N >= 0.
00047 *
00048 *     A       (input) DOUBLE PRECISION array, dimension (LDA,N)
00049 *     On entry, the N-by-N matrix A.
00050 *
00051 *     LDA     (input) INTEGER
00052 *     The leading dimension of the array A.  LDA >= max(1,N).
00053 *
00054 *     AF      (input) DOUBLE PRECISION array, dimension (LDAF,N)
00055 *     The block diagonal matrix D and the multipliers used to
00056 *     obtain the factor U or L as computed by DSYTRF.
00057 *
00058 *     LDAF    (input) INTEGER
00059 *     The leading dimension of the array AF.  LDAF >= max(1,N).
00060 *
00061 *     IPIV    (input) INTEGER array, dimension (N)
00062 *     Details of the interchanges and the block structure of D
00063 *     as determined by DSYTRF.
00064 *
00065 *     CMODE   (input) INTEGER
00066 *     Determines op2(C) in the formula op(A) * op2(C) as follows:
00067 *     CMODE =  1    op2(C) = C
00068 *     CMODE =  0    op2(C) = I
00069 *     CMODE = -1    op2(C) = inv(C)
00070 *
00071 *     C       (input) DOUBLE PRECISION array, dimension (N)
00072 *     The vector C in the formula op(A) * op2(C).
00073 *
00074 *     INFO    (output) INTEGER
00075 *       = 0:  Successful exit.
00076 *     i > 0:  The ith argument is invalid.
00077 *
00078 *     WORK    (input) DOUBLE PRECISION array, dimension (3*N).
00079 *     Workspace.
00080 *
00081 *     IWORK   (input) INTEGER array, dimension (N).
00082 *     Workspace.
00083 *
00084 *  =====================================================================
00085 *
00086 *     .. Local Scalars ..
00087       CHARACTER          NORMIN
00088       INTEGER            KASE, I, J
00089       DOUBLE PRECISION   AINVNM, SMLNUM, TMP
00090       LOGICAL            UP
00091 *     ..
00092 *     .. Local Arrays ..
00093       INTEGER            ISAVE( 3 )
00094 *     ..
00095 *     .. External Functions ..
00096       LOGICAL            LSAME
00097       INTEGER            IDAMAX
00098       DOUBLE PRECISION   DLAMCH
00099       EXTERNAL           LSAME, IDAMAX, DLAMCH
00100 *     ..
00101 *     .. External Subroutines ..
00102       EXTERNAL           DLACN2, DLATRS, DRSCL, XERBLA, DSYTRS
00103 *     ..
00104 *     .. Intrinsic Functions ..
00105       INTRINSIC          ABS, MAX
00106 *     ..
00107 *     .. Executable Statements ..
00108 *
00109       DLA_SYRCOND = 0.0D+0
00110 *
00111       INFO = 0
00112       IF( N.LT.0 ) THEN
00113          INFO = -2
00114       END IF
00115       IF( INFO.NE.0 ) THEN
00116          CALL XERBLA( 'DLA_SYRCOND', -INFO )
00117          RETURN
00118       END IF
00119       IF( N.EQ.0 ) THEN
00120          DLA_SYRCOND = 1.0D+0
00121          RETURN
00122       END IF
00123       UP = .FALSE.
00124       IF ( LSAME( UPLO, 'U' ) ) UP = .TRUE.
00125 *
00126 *     Compute the equilibration matrix R such that
00127 *     inv(R)*A*C has unit 1-norm.
00128 *
00129       IF ( UP ) THEN
00130          DO I = 1, N
00131             TMP = 0.0D+0
00132             IF ( CMODE .EQ. 1 ) THEN
00133                DO J = 1, I
00134                   TMP = TMP + ABS( A( J, I ) * C( J ) )
00135                END DO
00136                DO J = I+1, N
00137                   TMP = TMP + ABS( A( I, J ) * C( J ) )
00138                END DO
00139             ELSE IF ( CMODE .EQ. 0 ) THEN
00140                DO J = 1, I
00141                   TMP = TMP + ABS( A( J, I ) )
00142                END DO
00143                DO J = I+1, N
00144                   TMP = TMP + ABS( A( I, J ) )
00145                END DO
00146             ELSE
00147                DO J = 1, I
00148                   TMP = TMP + ABS( A( J, I ) / C( J ) )
00149                END DO
00150                DO J = I+1, N
00151                   TMP = TMP + ABS( A( I, J ) / C( J ) )
00152                END DO
00153             END IF
00154             WORK( 2*N+I ) = TMP
00155          END DO
00156       ELSE
00157          DO I = 1, N
00158             TMP = 0.0D+0
00159             IF ( CMODE .EQ. 1 ) THEN
00160                DO J = 1, I
00161                   TMP = TMP + ABS( A( I, J ) * C( J ) )
00162                END DO
00163                DO J = I+1, N
00164                   TMP = TMP + ABS( A( J, I ) * C( J ) )
00165                END DO
00166             ELSE IF ( CMODE .EQ. 0 ) THEN
00167                DO J = 1, I
00168                   TMP = TMP + ABS( A( I, J ) )
00169                END DO
00170                DO J = I+1, N
00171                   TMP = TMP + ABS( A( J, I ) )
00172                END DO
00173             ELSE
00174                DO J = 1, I
00175                   TMP = TMP + ABS( A( I, J) / C( J ) )
00176                END DO
00177                DO J = I+1, N
00178                   TMP = TMP + ABS( A( J, I) / C( J ) )
00179                END DO
00180             END IF
00181             WORK( 2*N+I ) = TMP
00182          END DO
00183       ENDIF
00184 *
00185 *     Estimate the norm of inv(op(A)).
00186 *
00187       SMLNUM = DLAMCH( 'Safe minimum' )
00188       AINVNM = 0.0D+0
00189       NORMIN = 'N'
00190 
00191       KASE = 0
00192    10 CONTINUE
00193       CALL DLACN2( N, WORK( N+1 ), WORK, IWORK, AINVNM, KASE, ISAVE )
00194       IF( KASE.NE.0 ) THEN
00195          IF( KASE.EQ.2 ) THEN
00196 *
00197 *           Multiply by R.
00198 *
00199             DO I = 1, N
00200                WORK( I ) = WORK( I ) * WORK( 2*N+I )
00201             END DO
00202 
00203             IF ( UP ) THEN
00204                CALL DSYTRS( 'U', N, 1, AF, LDAF, IPIV, WORK, N, INFO )
00205             ELSE
00206                CALL DSYTRS( 'L', N, 1, AF, LDAF, IPIV, WORK, N, INFO )
00207             ENDIF
00208 *
00209 *           Multiply by inv(C).
00210 *
00211             IF ( CMODE .EQ. 1 ) THEN
00212                DO I = 1, N
00213                   WORK( I ) = WORK( I ) / C( I )
00214                END DO
00215             ELSE IF ( CMODE .EQ. -1 ) THEN
00216                DO I = 1, N
00217                   WORK( I ) = WORK( I ) * C( I )
00218                END DO
00219             END IF
00220          ELSE
00221 *
00222 *           Multiply by inv(C').
00223 *
00224             IF ( CMODE .EQ. 1 ) THEN
00225                DO I = 1, N
00226                   WORK( I ) = WORK( I ) / C( I )
00227                END DO
00228             ELSE IF ( CMODE .EQ. -1 ) THEN
00229                DO I = 1, N
00230                   WORK( I ) = WORK( I ) * C( I )
00231                END DO
00232             END IF
00233 
00234             IF ( UP ) THEN
00235                CALL DSYTRS( 'U', N, 1, AF, LDAF, IPIV, WORK, N, INFO )
00236             ELSE
00237                CALL DSYTRS( 'L', N, 1, AF, LDAF, IPIV, WORK, N, INFO )
00238             ENDIF
00239 *
00240 *           Multiply by R.
00241 *
00242             DO I = 1, N
00243                WORK( I ) = WORK( I ) * WORK( 2*N+I )
00244             END DO
00245          END IF
00246 *
00247          GO TO 10
00248       END IF
00249 *
00250 *     Compute the estimate of the reciprocal condition number.
00251 *
00252       IF( AINVNM .NE. 0.0D+0 )
00253      $   DLA_SYRCOND = ( 1.0D+0 / AINVNM )
00254 *
00255       RETURN
00256 *
00257       END