001:       DOUBLE PRECISION FUNCTION DLA_SYRCOND( UPLO, N, A, LDA, AF, LDAF, 
002:      $                                       IPIV, CMODE, C, INFO, WORK,
003:      $                                       IWORK )
004: *
005: *     -- LAPACK routine (version 3.2.1)                                 --
006: *     -- Contributed by James Demmel, Deaglan Halligan, Yozo Hida and --
007: *     -- Jason Riedy of Univ. of California Berkeley.                 --
008: *     -- April 2009                                                   --
009: *
010: *     -- LAPACK is a software package provided by Univ. of Tennessee, --
011: *     -- Univ. of California Berkeley and NAG Ltd.                    --
012: *
013:       IMPLICIT NONE
014: *     ..
015: *     .. Scalar Arguments ..
016:       CHARACTER          UPLO
017:       INTEGER            N, LDA, LDAF, INFO, CMODE
018: *     ..
019: *     .. Array Arguments
020:       INTEGER            IWORK( * ), IPIV( * )
021:       DOUBLE PRECISION   A( LDA, * ), AF( LDAF, * ), WORK( * ), C( * )
022: *     ..
023: *
024: *  Purpose
025: *  =======
026: *
027: *     DLA_SYRCOND estimates the Skeel condition number of  op(A) * op2(C)
028: *     where op2 is determined by CMODE as follows
029: *     CMODE =  1    op2(C) = C
030: *     CMODE =  0    op2(C) = I
031: *     CMODE = -1    op2(C) = inv(C)
032: *     The Skeel condition number cond(A) = norminf( |inv(A)||A| )
033: *     is computed by computing scaling factors R such that
034: *     diag(R)*A*op2(C) is row equilibrated and computing the standard
035: *     infinity-norm condition number.
036: *
037: *  Arguments
038: *  ==========
039: *
040: *     UPLO    (input) CHARACTER*1
041: *       = 'U':  Upper triangle of A is stored;
042: *       = 'L':  Lower triangle of A is stored.
043: *
044: *     N       (input) INTEGER
045: *     The number of linear equations, i.e., the order of the
046: *     matrix A.  N >= 0.
047: *
048: *     A       (input) DOUBLE PRECISION array, dimension (LDA,N)
049: *     On entry, the N-by-N matrix A.
050: *
051: *     LDA     (input) INTEGER
052: *     The leading dimension of the array A.  LDA >= max(1,N).
053: *
054: *     AF      (input) DOUBLE PRECISION array, dimension (LDAF,N)
055: *     The block diagonal matrix D and the multipliers used to
056: *     obtain the factor U or L as computed by DSYTRF.
057: *
058: *     LDAF    (input) INTEGER
059: *     The leading dimension of the array AF.  LDAF >= max(1,N).
060: *
061: *     IPIV    (input) INTEGER array, dimension (N)
062: *     Details of the interchanges and the block structure of D
063: *     as determined by DSYTRF.
064: *
065: *     CMODE   (input) INTEGER
066: *     Determines op2(C) in the formula op(A) * op2(C) as follows:
067: *     CMODE =  1    op2(C) = C
068: *     CMODE =  0    op2(C) = I
069: *     CMODE = -1    op2(C) = inv(C)
070: *
071: *     C       (input) DOUBLE PRECISION array, dimension (N)
072: *     The vector C in the formula op(A) * op2(C).
073: *
074: *     INFO    (output) INTEGER
075: *       = 0:  Successful exit.
076: *     i > 0:  The ith argument is invalid.
077: *
078: *     WORK    (input) DOUBLE PRECISION array, dimension (3*N).
079: *     Workspace.
080: *
081: *     IWORK   (input) INTEGER array, dimension (N).
082: *     Workspace.
083: *
084: *  =====================================================================
085: *
086: *     .. Local Scalars ..
087:       CHARACTER          NORMIN
088:       INTEGER            KASE, I, J
089:       DOUBLE PRECISION   AINVNM, SMLNUM, TMP
090:       LOGICAL            UP
091: *     ..
092: *     .. Local Arrays ..
093:       INTEGER            ISAVE( 3 )
094: *     ..
095: *     .. External Functions ..
096:       LOGICAL            LSAME
097:       INTEGER            IDAMAX
098:       DOUBLE PRECISION   DLAMCH
099:       EXTERNAL           LSAME, IDAMAX, DLAMCH
100: *     ..
101: *     .. External Subroutines ..
102:       EXTERNAL           DLACN2, DLATRS, DRSCL, XERBLA, DSYTRS
103: *     ..
104: *     .. Intrinsic Functions ..
105:       INTRINSIC          ABS, MAX
106: *     ..
107: *     .. Executable Statements ..
108: *
109:       DLA_SYRCOND = 0.0D+0
110: *
111:       INFO = 0
112:       IF( N.LT.0 ) THEN
113:          INFO = -2
114:       END IF
115:       IF( INFO.NE.0 ) THEN
116:          CALL XERBLA( 'DLA_SYRCOND', -INFO )
117:          RETURN
118:       END IF
119:       IF( N.EQ.0 ) THEN
120:          DLA_SYRCOND = 1.0D+0
121:          RETURN
122:       END IF
123:       UP = .FALSE.
124:       IF ( LSAME( UPLO, 'U' ) ) UP = .TRUE.
125: *
126: *     Compute the equilibration matrix R such that
127: *     inv(R)*A*C has unit 1-norm.
128: *
129:       IF ( UP ) THEN
130:          DO I = 1, N
131:             TMP = 0.0D+0
132:             IF ( CMODE .EQ. 1 ) THEN
133:                DO J = 1, I
134:                   TMP = TMP + ABS( A( J, I ) * C( J ) )
135:                END DO
136:                DO J = I+1, N
137:                   TMP = TMP + ABS( A( I, J ) * C( J ) )
138:                END DO
139:             ELSE IF ( CMODE .EQ. 0 ) THEN
140:                DO J = 1, I
141:                   TMP = TMP + ABS( A( J, I ) )
142:                END DO
143:                DO J = I+1, N
144:                   TMP = TMP + ABS( A( I, J ) )
145:                END DO
146:             ELSE
147:                DO J = 1, I
148:                   TMP = TMP + ABS( A( J, I ) / C( J ) )
149:                END DO
150:                DO J = I+1, N
151:                   TMP = TMP + ABS( A( I, J ) / C( J ) )
152:                END DO
153:             END IF
154:             WORK( 2*N+I ) = TMP
155:          END DO
156:       ELSE
157:          DO I = 1, N
158:             TMP = 0.0D+0
159:             IF ( CMODE .EQ. 1 ) THEN
160:                DO J = 1, I
161:                   TMP = TMP + ABS( A( I, J ) * C( J ) )
162:                END DO
163:                DO J = I+1, N
164:                   TMP = TMP + ABS( A( J, I ) * C( J ) )
165:                END DO
166:             ELSE IF ( CMODE .EQ. 0 ) THEN
167:                DO J = 1, I
168:                   TMP = TMP + ABS( A( I, J ) )
169:                END DO
170:                DO J = I+1, N
171:                   TMP = TMP + ABS( A( J, I ) )
172:                END DO
173:             ELSE
174:                DO J = 1, I
175:                   TMP = TMP + ABS( A( I, J) / C( J ) )
176:                END DO
177:                DO J = I+1, N
178:                   TMP = TMP + ABS( A( J, I) / C( J ) )
179:                END DO
180:             END IF
181:             WORK( 2*N+I ) = TMP
182:          END DO
183:       ENDIF
184: *
185: *     Estimate the norm of inv(op(A)).
186: *
187:       SMLNUM = DLAMCH( 'Safe minimum' )
188:       AINVNM = 0.0D+0
189:       NORMIN = 'N'
190: 
191:       KASE = 0
192:    10 CONTINUE
193:       CALL DLACN2( N, WORK( N+1 ), WORK, IWORK, AINVNM, KASE, ISAVE )
194:       IF( KASE.NE.0 ) THEN
195:          IF( KASE.EQ.2 ) THEN
196: *
197: *           Multiply by R.
198: *
199:             DO I = 1, N
200:                WORK( I ) = WORK( I ) * WORK( 2*N+I )
201:             END DO
202: 
203:             IF ( UP ) THEN
204:                CALL DSYTRS( 'U', N, 1, AF, LDAF, IPIV, WORK, N, INFO )
205:             ELSE
206:                CALL DSYTRS( 'L', N, 1, AF, LDAF, IPIV, WORK, N, INFO )
207:             ENDIF
208: *
209: *           Multiply by inv(C).
210: *
211:             IF ( CMODE .EQ. 1 ) THEN
212:                DO I = 1, N
213:                   WORK( I ) = WORK( I ) / C( I )
214:                END DO
215:             ELSE IF ( CMODE .EQ. -1 ) THEN
216:                DO I = 1, N
217:                   WORK( I ) = WORK( I ) * C( I )
218:                END DO
219:             END IF
220:          ELSE
221: *
222: *           Multiply by inv(C').
223: *
224:             IF ( CMODE .EQ. 1 ) THEN
225:                DO I = 1, N
226:                   WORK( I ) = WORK( I ) / C( I )
227:                END DO
228:             ELSE IF ( CMODE .EQ. -1 ) THEN
229:                DO I = 1, N
230:                   WORK( I ) = WORK( I ) * C( I )
231:                END DO
232:             END IF
233: 
234:             IF ( UP ) THEN
235:                CALL DSYTRS( 'U', N, 1, AF, LDAF, IPIV, WORK, N, INFO )
236:             ELSE
237:                CALL DSYTRS( 'L', N, 1, AF, LDAF, IPIV, WORK, N, INFO )
238:             ENDIF
239: *
240: *           Multiply by R.
241: *
242:             DO I = 1, N
243:                WORK( I ) = WORK( I ) * WORK( 2*N+I )
244:             END DO
245:          END IF
246: *
247:          GO TO 10
248:       END IF
249: *
250: *     Compute the estimate of the reciprocal condition number.
251: *
252:       IF( AINVNM .NE. 0.0D+0 )
253:      $   DLA_SYRCOND = ( 1.0D+0 / AINVNM )
254: *
255:       RETURN
256: *
257:       END
258: