001:       REAL FUNCTION CLA_PORCOND_X( UPLO, N, A, LDA, AF, LDAF, X, INFO,
002:      $                             WORK, RWORK )
003: *
004: *     -- LAPACK routine (version 3.2.1)                                 --
005: *     -- Contributed by James Demmel, Deaglan Halligan, Yozo Hida and --
006: *     -- Jason Riedy of Univ. of California Berkeley.                 --
007: *     -- April 2009                                                   --
008: *
009: *     -- LAPACK is a software package provided by Univ. of Tennessee, --
010: *     -- Univ. of California Berkeley and NAG Ltd.                    --
011: *
012:       IMPLICIT NONE
013: *     ..
014: *     .. Scalar Arguments ..
015:       CHARACTER          UPLO
016:       INTEGER            N, LDA, LDAF, INFO
017: *     ..
018: *     .. Array Arguments ..
019:       COMPLEX            A( LDA, * ), AF( LDAF, * ), WORK( * ), X( * )
020:       REAL               RWORK( * )
021: *     ..
022: *
023: *  Purpose
024: *  =======
025: *
026: *     CLA_PORCOND_X Computes the infinity norm condition number of
027: *     op(A) * diag(X) where X is a COMPLEX vector.
028: *
029: *  Arguments
030: *  =========
031: *
032: *     UPLO    (input) CHARACTER*1
033: *       = 'U':  Upper triangle of A is stored;
034: *       = 'L':  Lower triangle of A is stored.
035: *
036: *     N       (input) INTEGER
037: *     The number of linear equations, i.e., the order of the
038: *     matrix A.  N >= 0.
039: *
040: *     A       (input) COMPLEX array, dimension (LDA,N)
041: *     On entry, the N-by-N matrix A.
042: *
043: *     LDA     (input) INTEGER
044: *     The leading dimension of the array A.  LDA >= max(1,N).
045: *
046: *     AF      (input) COMPLEX array, dimension (LDAF,N)
047: *     The triangular factor U or L from the Cholesky factorization
048: *     A = U**T*U or A = L*L**T, as computed by CPOTRF.
049: *
050: *     LDAF    (input) INTEGER
051: *     The leading dimension of the array AF.  LDAF >= max(1,N).
052: *
053: *     X       (input) COMPLEX array, dimension (N)
054: *     The vector X in the formula op(A) * diag(X).
055: *
056: *     INFO    (output) INTEGER
057: *       = 0:  Successful exit.
058: *     i > 0:  The ith argument is invalid.
059: *
060: *     WORK    (input) COMPLEX array, dimension (2*N).
061: *     Workspace.
062: *
063: *     RWORK   (input) REAL array, dimension (N).
064: *     Workspace.
065: *
066: *  =====================================================================
067: *
068: *     .. Local Scalars ..
069:       INTEGER            KASE, I, J
070:       REAL               AINVNM, ANORM, TMP
071:       LOGICAL            UP
072:       COMPLEX            ZDUM
073: *     ..
074: *     .. Local Arrays ..
075:       INTEGER            ISAVE( 3 )
076: *     ..
077: *     .. External Functions ..
078:       LOGICAL            LSAME
079:       EXTERNAL           LSAME
080: *     ..
081: *     .. External Subroutines ..
082:       EXTERNAL           CLACN2, CPOTRS, XERBLA
083: *     ..
084: *     .. Intrinsic Functions ..
085:       INTRINSIC          ABS, MAX, REAL, AIMAG
086: *     ..
087: *     .. Statement Functions ..
088:       REAL CABS1
089: *     ..
090: *     .. Statement Function Definitions ..
091:       CABS1( ZDUM ) = ABS( REAL( ZDUM ) ) + ABS( AIMAG( ZDUM ) )
092: *     ..
093: *     .. Executable Statements ..
094: *
095:       CLA_PORCOND_X = 0.0E+0
096: *
097:       INFO = 0
098:       IF( N.LT.0 ) THEN
099:          INFO = -2
100:       END IF
101:       IF( INFO.NE.0 ) THEN
102:          CALL XERBLA( 'CLA_PORCOND_X', -INFO )
103:          RETURN
104:       END IF
105:       UP = .FALSE.
106:       IF ( LSAME( UPLO, 'U' ) ) UP = .TRUE.
107: *
108: *     Compute norm of op(A)*op2(C).
109: *
110:       ANORM = 0.0
111:       IF ( UP ) THEN
112:          DO I = 1, N
113:             TMP = 0.0E+0
114:             DO J = 1, I
115:                TMP = TMP + CABS1( A( J, I ) * X( J ) )
116:             END DO
117:             DO J = I+1, N
118:                TMP = TMP + CABS1( A( I, J ) * X( J ) )
119:             END DO
120:             RWORK( I ) = TMP
121:             ANORM = MAX( ANORM, TMP )
122:          END DO
123:       ELSE
124:          DO I = 1, N
125:             TMP = 0.0E+0
126:             DO J = 1, I
127:                TMP = TMP + CABS1( A( I, J ) * X( J ) )
128:             END DO
129:             DO J = I+1, N
130:                TMP = TMP + CABS1( A( J, I ) * X( J ) )
131:             END DO
132:             RWORK( I ) = TMP
133:             ANORM = MAX( ANORM, TMP )
134:          END DO
135:       END IF
136: *
137: *     Quick return if possible.
138: *
139:       IF( N.EQ.0 ) THEN
140:          CLA_PORCOND_X = 1.0E+0
141:          RETURN
142:       ELSE IF( ANORM .EQ. 0.0E+0 ) THEN
143:          RETURN
144:       END IF
145: *
146: *     Estimate the norm of inv(op(A)).
147: *
148:       AINVNM = 0.0E+0
149: *
150:       KASE = 0
151:    10 CONTINUE
152:       CALL CLACN2( N, WORK( N+1 ), WORK, AINVNM, KASE, ISAVE )
153:       IF( KASE.NE.0 ) THEN
154:          IF( KASE.EQ.2 ) THEN
155: *
156: *           Multiply by R.
157: *
158:             DO I = 1, N
159:                WORK( I ) = WORK( I ) * RWORK( I )
160:             END DO
161: *
162:             IF ( UP ) THEN
163:                CALL CPOTRS( 'U', N, 1, AF, LDAF,
164:      $            WORK, N, INFO )
165:             ELSE
166:                CALL CPOTRS( 'L', N, 1, AF, LDAF,
167:      $            WORK, N, INFO )
168:             ENDIF
169: *
170: *           Multiply by inv(X).
171: *
172:             DO I = 1, N
173:                WORK( I ) = WORK( I ) / X( I )
174:             END DO
175:          ELSE
176: *
177: *           Multiply by inv(X').
178: *
179:             DO I = 1, N
180:                WORK( I ) = WORK( I ) / X( I )
181:             END DO
182: *
183:             IF ( UP ) THEN
184:                CALL CPOTRS( 'U', N, 1, AF, LDAF,
185:      $            WORK, N, INFO )
186:             ELSE
187:                CALL CPOTRS( 'L', N, 1, AF, LDAF,
188:      $            WORK, N, INFO )
189:             END IF
190: *
191: *           Multiply by R.
192: *
193:             DO I = 1, N
194:                WORK( I ) = WORK( I ) * RWORK( I )
195:             END DO
196:          END IF
197:          GO TO 10
198:       END IF
199: *
200: *     Compute the estimate of the reciprocal condition number.
201: *
202:       IF( AINVNM .NE. 0.0E+0 )
203:      $   CLA_PORCOND_X = 1.0E+0 / AINVNM
204: *
205:       RETURN
206: *
207:       END
208: