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