LAPACK: SRC/sla_porcond.f Source File

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