1 SUBROUTINE ZGETC2( 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 COMPLEX*16 A( LDA, * )
14 * ..
15 *
16 * Purpose
17 * =======
18 *
19 * ZGETC2 computes an LU factorization, using 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 a level 1 BLAS version of the algorithm.
25 *
26 * Arguments
27 * =========
28 *
29 * N (input) INTEGER
30 * The order of the matrix A. N >= 0.
31 *
32 * A (input/output) COMPLEX*16 array, dimension (LDA, N)
33 * On entry, the n-by-n matrix 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, 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 overflow if
53 * one tries to solve for x in Ax = b. So U is perturbed
54 * to 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 ZGERU, ZSWAP
75 * ..
76 * .. External Functions ..
77 DOUBLE PRECISION DLAMCH
78 EXTERNAL DLAMCH
79 * ..
80 * .. Intrinsic Functions ..
81 INTRINSIC ABS, DCMPLX, 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 ZSWAP( 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 ZSWAP( 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 ) = DCMPLX( SMIN, ZERO )
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 ZGERU( N-I, N-I, -DCMPLX( ONE ), A( I+1, I ), 1,
135 $ A( I, I+1 ), LDA, 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 ) = DCMPLX( SMIN, ZERO )
141 END IF
142 RETURN
143 *
144 * End of ZGETC2
145 *
146 END
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 COMPLEX*16 A( LDA, * )
14 * ..
15 *
16 * Purpose
17 * =======
18 *
19 * ZGETC2 computes an LU factorization, using 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 a level 1 BLAS version of the algorithm.
25 *
26 * Arguments
27 * =========
28 *
29 * N (input) INTEGER
30 * The order of the matrix A. N >= 0.
31 *
32 * A (input/output) COMPLEX*16 array, dimension (LDA, N)
33 * On entry, the n-by-n matrix 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, 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 overflow if
53 * one tries to solve for x in Ax = b. So U is perturbed
54 * to 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 ZGERU, ZSWAP
75 * ..
76 * .. External Functions ..
77 DOUBLE PRECISION DLAMCH
78 EXTERNAL DLAMCH
79 * ..
80 * .. Intrinsic Functions ..
81 INTRINSIC ABS, DCMPLX, 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 ZSWAP( 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 ZSWAP( 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 ) = DCMPLX( SMIN, ZERO )
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 ZGERU( N-I, N-I, -DCMPLX( ONE ), A( I+1, I ), 1,
135 $ A( I, I+1 ), LDA, 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 ) = DCMPLX( SMIN, ZERO )
141 END IF
142 RETURN
143 *
144 * End of ZGETC2
145 *
146 END