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