1 SUBROUTINE ZGESC2( N, A, LDA, RHS, IPIV, JPIV, SCALE )
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 LDA, N
10 DOUBLE PRECISION SCALE
11 * ..
12 * .. Array Arguments ..
13 INTEGER IPIV( * ), JPIV( * )
14 COMPLEX*16 A( LDA, * ), RHS( * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * ZGESC2 solves a system of linear equations
21 *
22 * A * X = scale* RHS
23 *
24 * with a general N-by-N matrix A using the LU factorization with
25 * complete pivoting computed by ZGETC2.
26 *
27 *
28 * Arguments
29 * =========
30 *
31 * N (input) INTEGER
32 * The number of columns of the matrix A.
33 *
34 * A (input) COMPLEX*16 array, dimension (LDA, N)
35 * On entry, the LU part of the factorization of the n-by-n
36 * matrix A computed by ZGETC2: A = P * L * U * Q
37 *
38 * LDA (input) INTEGER
39 * The leading dimension of the array A. LDA >= max(1, N).
40 *
41 * RHS (input/output) COMPLEX*16 array, dimension N.
42 * On entry, the right hand side vector b.
43 * On exit, the solution vector X.
44 *
45 * IPIV (input) INTEGER array, dimension (N).
46 * The pivot indices; for 1 <= i <= N, row i of the
47 * matrix has been interchanged with row IPIV(i).
48 *
49 * JPIV (input) INTEGER array, dimension (N).
50 * The pivot indices; for 1 <= j <= N, column j of the
51 * matrix has been interchanged with column JPIV(j).
52 *
53 * SCALE (output) DOUBLE PRECISION
54 * On exit, SCALE contains the scale factor. SCALE is chosen
55 * 0 <= SCALE <= 1 to prevent owerflow in the solution.
56 *
57 * Further Details
58 * ===============
59 *
60 * Based on contributions by
61 * Bo Kagstrom and Peter Poromaa, Department of Computing Science,
62 * Umea University, S-901 87 Umea, Sweden.
63 *
64 * =====================================================================
65 *
66 * .. Parameters ..
67 DOUBLE PRECISION ZERO, ONE, TWO
68 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0 )
69 * ..
70 * .. Local Scalars ..
71 INTEGER I, J
72 DOUBLE PRECISION BIGNUM, EPS, SMLNUM
73 COMPLEX*16 TEMP
74 * ..
75 * .. External Subroutines ..
76 EXTERNAL ZLASWP, ZSCAL
77 * ..
78 * .. External Functions ..
79 INTEGER IZAMAX
80 DOUBLE PRECISION DLAMCH
81 EXTERNAL IZAMAX, DLAMCH
82 * ..
83 * .. Intrinsic Functions ..
84 INTRINSIC ABS, DBLE, DCMPLX
85 * ..
86 * .. Executable Statements ..
87 *
88 * Set constant to control overflow
89 *
90 EPS = DLAMCH( 'P' )
91 SMLNUM = DLAMCH( 'S' ) / EPS
92 BIGNUM = ONE / SMLNUM
93 CALL DLABAD( SMLNUM, BIGNUM )
94 *
95 * Apply permutations IPIV to RHS
96 *
97 CALL ZLASWP( 1, RHS, LDA, 1, N-1, IPIV, 1 )
98 *
99 * Solve for L part
100 *
101 DO 20 I = 1, N - 1
102 DO 10 J = I + 1, N
103 RHS( J ) = RHS( J ) - A( J, I )*RHS( I )
104 10 CONTINUE
105 20 CONTINUE
106 *
107 * Solve for U part
108 *
109 SCALE = ONE
110 *
111 * Check for scaling
112 *
113 I = IZAMAX( N, RHS, 1 )
114 IF( TWO*SMLNUM*ABS( RHS( I ) ).GT.ABS( A( N, N ) ) ) THEN
115 TEMP = DCMPLX( ONE / TWO, ZERO ) / ABS( RHS( I ) )
116 CALL ZSCAL( N, TEMP, RHS( 1 ), 1 )
117 SCALE = SCALE*DBLE( TEMP )
118 END IF
119 DO 40 I = N, 1, -1
120 TEMP = DCMPLX( ONE, ZERO ) / A( I, I )
121 RHS( I ) = RHS( I )*TEMP
122 DO 30 J = I + 1, N
123 RHS( I ) = RHS( I ) - RHS( J )*( A( I, J )*TEMP )
124 30 CONTINUE
125 40 CONTINUE
126 *
127 * Apply permutations JPIV to the solution (RHS)
128 *
129 CALL ZLASWP( 1, RHS, LDA, 1, N-1, JPIV, -1 )
130 RETURN
131 *
132 * End of ZGESC2
133 *
134 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 LDA, N
10 DOUBLE PRECISION SCALE
11 * ..
12 * .. Array Arguments ..
13 INTEGER IPIV( * ), JPIV( * )
14 COMPLEX*16 A( LDA, * ), RHS( * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * ZGESC2 solves a system of linear equations
21 *
22 * A * X = scale* RHS
23 *
24 * with a general N-by-N matrix A using the LU factorization with
25 * complete pivoting computed by ZGETC2.
26 *
27 *
28 * Arguments
29 * =========
30 *
31 * N (input) INTEGER
32 * The number of columns of the matrix A.
33 *
34 * A (input) COMPLEX*16 array, dimension (LDA, N)
35 * On entry, the LU part of the factorization of the n-by-n
36 * matrix A computed by ZGETC2: A = P * L * U * Q
37 *
38 * LDA (input) INTEGER
39 * The leading dimension of the array A. LDA >= max(1, N).
40 *
41 * RHS (input/output) COMPLEX*16 array, dimension N.
42 * On entry, the right hand side vector b.
43 * On exit, the solution vector X.
44 *
45 * IPIV (input) INTEGER array, dimension (N).
46 * The pivot indices; for 1 <= i <= N, row i of the
47 * matrix has been interchanged with row IPIV(i).
48 *
49 * JPIV (input) INTEGER array, dimension (N).
50 * The pivot indices; for 1 <= j <= N, column j of the
51 * matrix has been interchanged with column JPIV(j).
52 *
53 * SCALE (output) DOUBLE PRECISION
54 * On exit, SCALE contains the scale factor. SCALE is chosen
55 * 0 <= SCALE <= 1 to prevent owerflow in the solution.
56 *
57 * Further Details
58 * ===============
59 *
60 * Based on contributions by
61 * Bo Kagstrom and Peter Poromaa, Department of Computing Science,
62 * Umea University, S-901 87 Umea, Sweden.
63 *
64 * =====================================================================
65 *
66 * .. Parameters ..
67 DOUBLE PRECISION ZERO, ONE, TWO
68 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0 )
69 * ..
70 * .. Local Scalars ..
71 INTEGER I, J
72 DOUBLE PRECISION BIGNUM, EPS, SMLNUM
73 COMPLEX*16 TEMP
74 * ..
75 * .. External Subroutines ..
76 EXTERNAL ZLASWP, ZSCAL
77 * ..
78 * .. External Functions ..
79 INTEGER IZAMAX
80 DOUBLE PRECISION DLAMCH
81 EXTERNAL IZAMAX, DLAMCH
82 * ..
83 * .. Intrinsic Functions ..
84 INTRINSIC ABS, DBLE, DCMPLX
85 * ..
86 * .. Executable Statements ..
87 *
88 * Set constant to control overflow
89 *
90 EPS = DLAMCH( 'P' )
91 SMLNUM = DLAMCH( 'S' ) / EPS
92 BIGNUM = ONE / SMLNUM
93 CALL DLABAD( SMLNUM, BIGNUM )
94 *
95 * Apply permutations IPIV to RHS
96 *
97 CALL ZLASWP( 1, RHS, LDA, 1, N-1, IPIV, 1 )
98 *
99 * Solve for L part
100 *
101 DO 20 I = 1, N - 1
102 DO 10 J = I + 1, N
103 RHS( J ) = RHS( J ) - A( J, I )*RHS( I )
104 10 CONTINUE
105 20 CONTINUE
106 *
107 * Solve for U part
108 *
109 SCALE = ONE
110 *
111 * Check for scaling
112 *
113 I = IZAMAX( N, RHS, 1 )
114 IF( TWO*SMLNUM*ABS( RHS( I ) ).GT.ABS( A( N, N ) ) ) THEN
115 TEMP = DCMPLX( ONE / TWO, ZERO ) / ABS( RHS( I ) )
116 CALL ZSCAL( N, TEMP, RHS( 1 ), 1 )
117 SCALE = SCALE*DBLE( TEMP )
118 END IF
119 DO 40 I = N, 1, -1
120 TEMP = DCMPLX( ONE, ZERO ) / A( I, I )
121 RHS( I ) = RHS( I )*TEMP
122 DO 30 J = I + 1, N
123 RHS( I ) = RHS( I ) - RHS( J )*( A( I, J )*TEMP )
124 30 CONTINUE
125 40 CONTINUE
126 *
127 * Apply permutations JPIV to the solution (RHS)
128 *
129 CALL ZLASWP( 1, RHS, LDA, 1, N-1, JPIV, -1 )
130 RETURN
131 *
132 * End of ZGESC2
133 *
134 END