LAPACK: SRC/cla_gercond_x.f Source File

Go to the documentation of this file.
00001       REAL FUNCTION CLA_GERCOND_X( TRANS, N, A, LDA, AF, LDAF, IPIV, X,
00002      $                             INFO, WORK, RWORK )
00003 *
00004 *     -- LAPACK routine (version 3.2.1)                                 --
00005 *     -- Contributed by James Demmel, Deaglan Halligan, Yozo Hida and --
00006 *     -- Jason Riedy of Univ. of California Berkeley.                 --
00007 *     -- April 2009                                                   --
00008 *
00009 *     -- LAPACK is a software package provided by Univ. of Tennessee, --
00010 *     -- Univ. of California Berkeley and NAG Ltd.                    --
00011 *
00012       IMPLICIT NONE
00013 *     ..
00014 *     .. Scalar Arguments ..
00015       CHARACTER          TRANS
00016       INTEGER            N, LDA, LDAF, INFO
00017 *     ..
00018 *     .. Array Arguments ..
00019       INTEGER            IPIV( * )
00020       COMPLEX            A( LDA, * ), AF( LDAF, * ), WORK( * ), X( * )
00021       REAL               RWORK( * )
00022 *     ..
00023 *
00024 *  Purpose
00025 *  =======
00026 * 
00027 *     CLA_GERCOND_X computes the infinity norm condition number of
00028 *     op(A) * diag(X) where X is a COMPLEX vector.
00029 *
00030 *  Arguments
00031 *  =========
00032 *
00033 *     TRANS   (input) CHARACTER*1
00034 *     Specifies the form of the system of equations:
00035 *       = 'N':  A * X = B     (No transpose)
00036 *       = 'T':  A**T * X = B  (Transpose)
00037 *       = 'C':  A**H * X = B  (Conjugate Transpose = Transpose)
00038 *
00039 *     N       (input) INTEGER
00040 *     The number of linear equations, i.e., the order of the
00041 *     matrix A.  N >= 0.
00042 *
00043 *     A       (input) COMPLEX array, dimension (LDA,N)
00044 *     On entry, the N-by-N matrix A.
00045 *
00046 *     LDA     (input) INTEGER
00047 *     The leading dimension of the array A.  LDA >= max(1,N).
00048 *
00049 *     AF      (input) COMPLEX array, dimension (LDAF,N)
00050 *     The factors L and U from the factorization
00051 *     A = P*L*U as computed by CGETRF.
00052 *
00053 *     LDAF    (input) INTEGER
00054 *     The leading dimension of the array AF.  LDAF >= max(1,N).
00055 *
00056 *     IPIV    (input) INTEGER array, dimension (N)
00057 *     The pivot indices from the factorization A = P*L*U
00058 *     as computed by CGETRF; row i of the matrix was interchanged
00059 *     with row IPIV(i).
00060 *
00061 *     X       (input) COMPLEX array, dimension (N)
00062 *     The vector X in the formula op(A) * diag(X).
00063 *
00064 *     INFO    (output) INTEGER
00065 *       = 0:  Successful exit.
00066 *     i > 0:  The ith argument is invalid.
00067 *
00068 *     WORK    (input) COMPLEX array, dimension (2*N).
00069 *     Workspace.
00070 *
00071 *     RWORK   (input) REAL array, dimension (N).
00072 *     Workspace.
00073 *
00074 *  =====================================================================
00075 *
00076 *     .. Local Scalars ..
00077       LOGICAL            NOTRANS
00078       INTEGER            KASE
00079       REAL               AINVNM, ANORM, TMP
00080       INTEGER            I, J
00081       COMPLEX            ZDUM
00082 *     ..
00083 *     .. Local Arrays ..
00084       INTEGER            ISAVE( 3 )
00085 *     ..
00086 *     .. External Functions ..
00087       LOGICAL            LSAME
00088       EXTERNAL           LSAME
00089 *     ..
00090 *     .. External Subroutines ..
00091       EXTERNAL           CLACN2, CGETRS, XERBLA
00092 *     ..
00093 *     .. Intrinsic Functions ..
00094       INTRINSIC          ABS, MAX, REAL, AIMAG
00095 *     ..
00096 *     .. Statement Functions ..
00097       REAL               CABS1
00098 *     ..
00099 *     .. Statement Function Definitions ..
00100       CABS1( ZDUM ) = ABS( REAL( ZDUM ) ) + ABS( AIMAG( ZDUM ) )
00101 *     ..
00102 *     .. Executable Statements ..
00103 *
00104       CLA_GERCOND_X = 0.0E+0
00105 *
00106       INFO = 0
00107       NOTRANS = LSAME( TRANS, 'N' )
00108       IF ( .NOT. NOTRANS .AND. .NOT. LSAME( TRANS, 'T' ) .AND. .NOT.
00109      $     LSAME( TRANS, 'C' ) ) THEN
00110          INFO = -1
00111       ELSE IF( N.LT.0 ) THEN
00112          INFO = -2
00113       END IF
00114       IF( INFO.NE.0 ) THEN
00115          CALL XERBLA( 'CLA_GERCOND_X', -INFO )
00116          RETURN
00117       END IF
00118 *
00119 *     Compute norm of op(A)*op2(C).
00120 *
00121       ANORM = 0.0
00122       IF ( NOTRANS ) THEN
00123          DO I = 1, N
00124             TMP = 0.0E+0
00125             DO J = 1, N
00126                TMP = TMP + CABS1( A( I, J ) * X( J ) )
00127             END DO
00128             RWORK( I ) = TMP
00129             ANORM = MAX( ANORM, TMP )
00130          END DO
00131       ELSE
00132          DO I = 1, N
00133             TMP = 0.0E+0
00134             DO J = 1, N
00135                TMP = TMP + CABS1( A( J, I ) * X( J ) )
00136             END DO
00137             RWORK( I ) = TMP
00138             ANORM = MAX( ANORM, TMP )
00139          END DO
00140       END IF
00141 *
00142 *     Quick return if possible.
00143 *
00144       IF( N.EQ.0 ) THEN
00145          CLA_GERCOND_X = 1.0E+0
00146          RETURN
00147       ELSE IF( ANORM .EQ. 0.0E+0 ) THEN
00148          RETURN
00149       END IF
00150 *
00151 *     Estimate the norm of inv(op(A)).
00152 *
00153       AINVNM = 0.0E+0
00154 *
00155       KASE = 0
00156    10 CONTINUE
00157       CALL CLACN2( N, WORK( N+1 ), WORK, AINVNM, KASE, ISAVE )
00158       IF( KASE.NE.0 ) THEN
00159          IF( KASE.EQ.2 ) THEN
00160 *           Multiply by R.
00161             DO I = 1, N
00162                WORK( I ) = WORK( I ) * RWORK( I )
00163             END DO
00164 *
00165             IF ( NOTRANS ) THEN
00166                CALL CGETRS( 'No transpose', N, 1, AF, LDAF, IPIV,
00167      $            WORK, N, INFO )
00168             ELSE
00169                CALL CGETRS( 'Conjugate transpose', N, 1, AF, LDAF, IPIV,
00170      $            WORK, N, INFO )
00171             ENDIF
00172 *
00173 *           Multiply by inv(X).
00174 *
00175             DO I = 1, N
00176                WORK( I ) = WORK( I ) / X( I )
00177             END DO
00178          ELSE
00179 *
00180 *           Multiply by inv(X').
00181 *
00182             DO I = 1, N
00183                WORK( I ) = WORK( I ) / X( I )
00184             END DO
00185 *
00186             IF ( NOTRANS ) THEN
00187                CALL CGETRS( 'Conjugate transpose', N, 1, AF, LDAF, IPIV,
00188      $            WORK, N, INFO )
00189             ELSE
00190                CALL CGETRS( 'No transpose', N, 1, AF, LDAF, IPIV,
00191      $            WORK, N, INFO )
00192             END IF
00193 *
00194 *           Multiply by R.
00195 *
00196             DO I = 1, N
00197                WORK( I ) = WORK( I ) * RWORK( I )
00198             END DO
00199          END IF
00200          GO TO 10
00201       END IF
00202 *
00203 *     Compute the estimate of the reciprocal condition number.
00204 *
00205       IF( AINVNM .NE. 0.0E+0 )
00206      $   CLA_GERCOND_X = 1.0E+0 / AINVNM
00207 *
00208       RETURN
00209 *
00210       END