|
1
2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 |
SUBROUTINE SGETC2( N, A, LDA, IPIV, JPIV, INFO )
* * -- LAPACK auxiliary routine (version 3.2) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * November 2006 * * .. Scalar Arguments .. INTEGER INFO, LDA, N * .. * .. Array Arguments .. INTEGER IPIV( * ), JPIV( * ) REAL A( LDA, * ) * .. * * Purpose * ======= * * SGETC2 computes an LU factorization with complete pivoting of the * n-by-n matrix A. The factorization has the form A = P * L * U * Q, * where P and Q are permutation matrices, L is lower triangular with * unit diagonal elements and U is upper triangular. * * This is the Level 2 BLAS algorithm. * * Arguments * ========= * * N (input) INTEGER * The order of the matrix A. N >= 0. * * A (input/output) REAL array, dimension (LDA, N) * On entry, the n-by-n matrix A to be factored. * On exit, the factors L and U from the factorization * A = P*L*U*Q; the unit diagonal elements of L are not stored. * If U(k, k) appears to be less than SMIN, U(k, k) is given the * value of SMIN, i.e., giving a nonsingular perturbed system. * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,N). * * IPIV (output) INTEGER array, dimension(N). * The pivot indices; for 1 <= i <= N, row i of the * matrix has been interchanged with row IPIV(i). * * JPIV (output) INTEGER array, dimension(N). * The pivot indices; for 1 <= j <= N, column j of the * matrix has been interchanged with column JPIV(j). * * INFO (output) INTEGER * = 0: successful exit * > 0: if INFO = k, U(k, k) is likely to produce owerflow if * we try to solve for x in Ax = b. So U is perturbed to * avoid the overflow. * * Further Details * =============== * * Based on contributions by * Bo Kagstrom and Peter Poromaa, Department of Computing Science, * Umea University, S-901 87 Umea, Sweden. * * ===================================================================== * * .. Parameters .. REAL ZERO, ONE PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) * .. * .. Local Scalars .. INTEGER I, IP, IPV, J, JP, JPV REAL BIGNUM, EPS, SMIN, SMLNUM, XMAX * .. * .. External Subroutines .. EXTERNAL SGER, SLABAD, SSWAP * .. * .. External Functions .. REAL SLAMCH EXTERNAL SLAMCH * .. * .. Intrinsic Functions .. INTRINSIC ABS, MAX * .. * .. Executable Statements .. * * Set constants to control overflow * INFO = 0 EPS = SLAMCH( 'P' ) SMLNUM = SLAMCH( 'S' ) / EPS BIGNUM = ONE / SMLNUM CALL SLABAD( SMLNUM, BIGNUM ) * * Factorize A using complete pivoting. * Set pivots less than SMIN to SMIN. * DO 40 I = 1, N - 1 * * Find max element in matrix A * XMAX = ZERO DO 20 IP = I, N DO 10 JP = I, N IF( ABS( A( IP, JP ) ).GE.XMAX ) THEN XMAX = ABS( A( IP, JP ) ) IPV = IP JPV = JP END IF 10 CONTINUE 20 CONTINUE IF( I.EQ.1 ) $ SMIN = MAX( EPS*XMAX, SMLNUM ) * * Swap rows * IF( IPV.NE.I ) $ CALL SSWAP( N, A( IPV, 1 ), LDA, A( I, 1 ), LDA ) IPIV( I ) = IPV * * Swap columns * IF( JPV.NE.I ) $ CALL SSWAP( N, A( 1, JPV ), 1, A( 1, I ), 1 ) JPIV( I ) = JPV * * Check for singularity * IF( ABS( A( I, I ) ).LT.SMIN ) THEN INFO = I A( I, I ) = SMIN END IF DO 30 J = I + 1, N A( J, I ) = A( J, I ) / A( I, I ) 30 CONTINUE CALL SGER( N-I, N-I, -ONE, A( I+1, I ), 1, A( I, I+1 ), LDA, $ A( I+1, I+1 ), LDA ) 40 CONTINUE * IF( ABS( A( N, N ) ).LT.SMIN ) THEN INFO = N A( N, N ) = SMIN END IF * RETURN * * End of SGETC2 * END |