Annotation of rpl/lapack/lapack/zgesc2.f, revision 1.7
1.1 bertrand 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
CVSweb interface <joel.bertrand@systella.fr>