001:       SUBROUTINE SGETC2( N, A, LDA, IPIV, JPIV, INFO )
002: *
003: *  -- LAPACK auxiliary routine (version 3.2) --
004: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
005: *     November 2006
006: *
007: *     .. Scalar Arguments ..
008:       INTEGER            INFO, LDA, N
009: *     ..
010: *     .. Array Arguments ..
011:       INTEGER            IPIV( * ), JPIV( * )
012:       REAL               A( LDA, * )
013: *     ..
014: *
015: *  Purpose
016: *  =======
017: *
018: *  SGETC2 computes an LU factorization with complete pivoting of the
019: *  n-by-n matrix A. The factorization has the form A = P * L * U * Q,
020: *  where P and Q are permutation matrices, L is lower triangular with
021: *  unit diagonal elements and U is upper triangular.
022: *
023: *  This is the Level 2 BLAS algorithm.
024: *
025: *  Arguments
026: *  =========
027: *
028: *  N       (input) INTEGER
029: *          The order of the matrix A. N >= 0.
030: *
031: *  A       (input/output) REAL array, dimension (LDA, N)
032: *          On entry, the n-by-n matrix A to be factored.
033: *          On exit, the factors L and U from the factorization
034: *          A = P*L*U*Q; the unit diagonal elements of L are not stored.
035: *          If U(k, k) appears to be less than SMIN, U(k, k) is given the
036: *          value of SMIN, i.e., giving a nonsingular perturbed system.
037: *
038: *  LDA     (input) INTEGER
039: *          The leading dimension of the array A.  LDA >= max(1,N).
040: *
041: *  IPIV    (output) INTEGER array, dimension(N).
042: *          The pivot indices; for 1 <= i <= N, row i of the
043: *          matrix has been interchanged with row IPIV(i).
044: *
045: *  JPIV    (output) INTEGER array, dimension(N).
046: *          The pivot indices; for 1 <= j <= N, column j of the
047: *          matrix has been interchanged with column JPIV(j).
048: *
049: *  INFO    (output) INTEGER
050: *           = 0: successful exit
051: *           > 0: if INFO = k, U(k, k) is likely to produce owerflow if
052: *                we try to solve for x in Ax = b. So U is perturbed to
053: *                avoid the overflow.
054: *
055: *  Further Details
056: *  ===============
057: *
058: *  Based on contributions by
059: *     Bo Kagstrom and Peter Poromaa, Department of Computing Science,
060: *     Umea University, S-901 87 Umea, Sweden.
061: *
062: *  =====================================================================
063: *
064: *     .. Parameters ..
065:       REAL               ZERO, ONE
066:       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
067: *     ..
068: *     .. Local Scalars ..
069:       INTEGER            I, IP, IPV, J, JP, JPV
070:       REAL               BIGNUM, EPS, SMIN, SMLNUM, XMAX
071: *     ..
072: *     .. External Subroutines ..
073:       EXTERNAL           SGER, SLABAD, SSWAP
074: *     ..
075: *     .. External Functions ..
076:       REAL               SLAMCH
077:       EXTERNAL           SLAMCH
078: *     ..
079: *     .. Intrinsic Functions ..
080:       INTRINSIC          ABS, MAX
081: *     ..
082: *     .. Executable Statements ..
083: *
084: *     Set constants to control overflow
085: *
086:       INFO = 0
087:       EPS = SLAMCH( 'P' )
088:       SMLNUM = SLAMCH( 'S' ) / EPS
089:       BIGNUM = ONE / SMLNUM
090:       CALL SLABAD( SMLNUM, BIGNUM )
091: *
092: *     Factorize A using complete pivoting.
093: *     Set pivots less than SMIN to SMIN.
094: *
095:       DO 40 I = 1, N - 1
096: *
097: *        Find max element in matrix A
098: *
099:          XMAX = ZERO
100:          DO 20 IP = I, N
101:             DO 10 JP = I, N
102:                IF( ABS( A( IP, JP ) ).GE.XMAX ) THEN
103:                   XMAX = ABS( A( IP, JP ) )
104:                   IPV = IP
105:                   JPV = JP
106:                END IF
107:    10       CONTINUE
108:    20    CONTINUE
109:          IF( I.EQ.1 )
110:      $      SMIN = MAX( EPS*XMAX, SMLNUM )
111: *
112: *        Swap rows
113: *
114:          IF( IPV.NE.I )
115:      $      CALL SSWAP( N, A( IPV, 1 ), LDA, A( I, 1 ), LDA )
116:          IPIV( I ) = IPV
117: *
118: *        Swap columns
119: *
120:          IF( JPV.NE.I )
121:      $      CALL SSWAP( N, A( 1, JPV ), 1, A( 1, I ), 1 )
122:          JPIV( I ) = JPV
123: *
124: *        Check for singularity
125: *
126:          IF( ABS( A( I, I ) ).LT.SMIN ) THEN
127:             INFO = I
128:             A( I, I ) = SMIN
129:          END IF
130:          DO 30 J = I + 1, N
131:             A( J, I ) = A( J, I ) / A( I, I )
132:    30    CONTINUE
133:          CALL SGER( N-I, N-I, -ONE, A( I+1, I ), 1, A( I, I+1 ), LDA,
134:      $              A( I+1, I+1 ), LDA )
135:    40 CONTINUE
136: *
137:       IF( ABS( A( N, N ) ).LT.SMIN ) THEN
138:          INFO = N
139:          A( N, N ) = SMIN
140:       END IF
141: *
142:       RETURN
143: *
144: *     End of SGETC2
145: *
146:       END
147: