LAPACK: SRC/dla_porcond.f Source File

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