1: SUBROUTINE DGETC2( N, A, LDA, IPIV, JPIV, INFO )
2: *
3: * -- LAPACK auxiliary routine (version 3.2) --
4: * -- LAPACK is a software package provided by Univ. of Tennessee, --
5: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6: * November 2006
7: *
8: * .. Scalar Arguments ..
9: INTEGER INFO, LDA, N
10: * ..
11: * .. Array Arguments ..
12: INTEGER IPIV( * ), JPIV( * )
13: DOUBLE PRECISION A( LDA, * )
14: * ..
15: *
16: * Purpose
17: * =======
18: *
19: * DGETC2 computes an LU factorization with complete pivoting of the
20: * n-by-n matrix A. The factorization has the form A = P * L * U * Q,
21: * where P and Q are permutation matrices, L is lower triangular with
22: * unit diagonal elements and U is upper triangular.
23: *
24: * This is the Level 2 BLAS algorithm.
25: *
26: * Arguments
27: * =========
28: *
29: * N (input) INTEGER
30: * The order of the matrix A. N >= 0.
31: *
32: * A (input/output) DOUBLE PRECISION array, dimension (LDA, N)
33: * On entry, the n-by-n matrix A to be factored.
34: * On exit, the factors L and U from the factorization
35: * A = P*L*U*Q; the unit diagonal elements of L are not stored.
36: * If U(k, k) appears to be less than SMIN, U(k, k) is given the
37: * value of SMIN, i.e., giving a nonsingular perturbed system.
38: *
39: * LDA (input) INTEGER
40: * The leading dimension of the array A. LDA >= max(1,N).
41: *
42: * IPIV (output) INTEGER array, dimension(N).
43: * The pivot indices; for 1 <= i <= N, row i of the
44: * matrix has been interchanged with row IPIV(i).
45: *
46: * JPIV (output) INTEGER array, dimension(N).
47: * The pivot indices; for 1 <= j <= N, column j of the
48: * matrix has been interchanged with column JPIV(j).
49: *
50: * INFO (output) INTEGER
51: * = 0: successful exit
52: * > 0: if INFO = k, U(k, k) is likely to produce owerflow if
53: * we try to solve for x in Ax = b. So U is perturbed to
54: * avoid the overflow.
55: *
56: * Further Details
57: * ===============
58: *
59: * Based on contributions by
60: * Bo Kagstrom and Peter Poromaa, Department of Computing Science,
61: * Umea University, S-901 87 Umea, Sweden.
62: *
63: * =====================================================================
64: *
65: * .. Parameters ..
66: DOUBLE PRECISION ZERO, ONE
67: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
68: * ..
69: * .. Local Scalars ..
70: INTEGER I, IP, IPV, J, JP, JPV
71: DOUBLE PRECISION BIGNUM, EPS, SMIN, SMLNUM, XMAX
72: * ..
73: * .. External Subroutines ..
74: EXTERNAL DGER, DSWAP
75: * ..
76: * .. External Functions ..
77: DOUBLE PRECISION DLAMCH
78: EXTERNAL DLAMCH
79: * ..
80: * .. Intrinsic Functions ..
81: INTRINSIC ABS, MAX
82: * ..
83: * .. Executable Statements ..
84: *
85: * Set constants to control overflow
86: *
87: INFO = 0
88: EPS = DLAMCH( 'P' )
89: SMLNUM = DLAMCH( 'S' ) / EPS
90: BIGNUM = ONE / SMLNUM
91: CALL DLABAD( SMLNUM, BIGNUM )
92: *
93: * Factorize A using complete pivoting.
94: * Set pivots less than SMIN to SMIN.
95: *
96: DO 40 I = 1, N - 1
97: *
98: * Find max element in matrix A
99: *
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 DSWAP( 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 DSWAP( 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 DGER( 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 DGETC2
146: *
147: END
CVSweb interface <joel.bertrand@systella.fr>