00001 SUBROUTINE CGESC2( N, A, LDA, RHS, IPIV, JPIV, SCALE ) 00002 * 00003 * -- LAPACK auxiliary routine (version 3.2) -- 00004 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00005 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00006 * November 2006 00007 * 00008 * .. Scalar Arguments .. 00009 INTEGER LDA, N 00010 REAL SCALE 00011 * .. 00012 * .. Array Arguments .. 00013 INTEGER IPIV( * ), JPIV( * ) 00014 COMPLEX A( LDA, * ), RHS( * ) 00015 * .. 00016 * 00017 * Purpose 00018 * ======= 00019 * 00020 * CGESC2 solves a system of linear equations 00021 * 00022 * A * X = scale* RHS 00023 * 00024 * with a general N-by-N matrix A using the LU factorization with 00025 * complete pivoting computed by CGETC2. 00026 * 00027 * 00028 * Arguments 00029 * ========= 00030 * 00031 * N (input) INTEGER 00032 * The number of columns of the matrix A. 00033 * 00034 * A (input) COMPLEX array, dimension (LDA, N) 00035 * On entry, the LU part of the factorization of the n-by-n 00036 * matrix A computed by CGETC2: A = P * L * U * Q 00037 * 00038 * LDA (input) INTEGER 00039 * The leading dimension of the array A. LDA >= max(1, N). 00040 * 00041 * RHS (input/output) COMPLEX array, dimension N. 00042 * On entry, the right hand side vector b. 00043 * On exit, the solution vector X. 00044 * 00045 * IPIV (input) INTEGER array, dimension (N). 00046 * The pivot indices; for 1 <= i <= N, row i of the 00047 * matrix has been interchanged with row IPIV(i). 00048 * 00049 * JPIV (input) INTEGER array, dimension (N). 00050 * The pivot indices; for 1 <= j <= N, column j of the 00051 * matrix has been interchanged with column JPIV(j). 00052 * 00053 * SCALE (output) REAL 00054 * On exit, SCALE contains the scale factor. SCALE is chosen 00055 * 0 <= SCALE <= 1 to prevent owerflow in the solution. 00056 * 00057 * Further Details 00058 * =============== 00059 * 00060 * Based on contributions by 00061 * Bo Kagstrom and Peter Poromaa, Department of Computing Science, 00062 * Umea University, S-901 87 Umea, Sweden. 00063 * 00064 * ===================================================================== 00065 * 00066 * .. Parameters .. 00067 REAL ZERO, ONE, TWO 00068 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0, TWO = 2.0E+0 ) 00069 * .. 00070 * .. Local Scalars .. 00071 INTEGER I, J 00072 REAL BIGNUM, EPS, SMLNUM 00073 COMPLEX TEMP 00074 * .. 00075 * .. External Subroutines .. 00076 EXTERNAL CLASWP, CSCAL, SLABAD 00077 * .. 00078 * .. External Functions .. 00079 INTEGER ICAMAX 00080 REAL SLAMCH 00081 EXTERNAL ICAMAX, SLAMCH 00082 * .. 00083 * .. Intrinsic Functions .. 00084 INTRINSIC ABS, CMPLX, REAL 00085 * .. 00086 * .. Executable Statements .. 00087 * 00088 * Set constant to control overflow 00089 * 00090 EPS = SLAMCH( 'P' ) 00091 SMLNUM = SLAMCH( 'S' ) / EPS 00092 BIGNUM = ONE / SMLNUM 00093 CALL SLABAD( SMLNUM, BIGNUM ) 00094 * 00095 * Apply permutations IPIV to RHS 00096 * 00097 CALL CLASWP( 1, RHS, LDA, 1, N-1, IPIV, 1 ) 00098 * 00099 * Solve for L part 00100 * 00101 DO 20 I = 1, N - 1 00102 DO 10 J = I + 1, N 00103 RHS( J ) = RHS( J ) - A( J, I )*RHS( I ) 00104 10 CONTINUE 00105 20 CONTINUE 00106 * 00107 * Solve for U part 00108 * 00109 SCALE = ONE 00110 * 00111 * Check for scaling 00112 * 00113 I = ICAMAX( N, RHS, 1 ) 00114 IF( TWO*SMLNUM*ABS( RHS( I ) ).GT.ABS( A( N, N ) ) ) THEN 00115 TEMP = CMPLX( ONE / TWO, ZERO ) / ABS( RHS( I ) ) 00116 CALL CSCAL( N, TEMP, RHS( 1 ), 1 ) 00117 SCALE = SCALE*REAL( TEMP ) 00118 END IF 00119 DO 40 I = N, 1, -1 00120 TEMP = CMPLX( ONE, ZERO ) / A( I, I ) 00121 RHS( I ) = RHS( I )*TEMP 00122 DO 30 J = I + 1, N 00123 RHS( I ) = RHS( I ) - RHS( J )*( A( I, J )*TEMP ) 00124 30 CONTINUE 00125 40 CONTINUE 00126 * 00127 * Apply permutations JPIV to the solution (RHS) 00128 * 00129 CALL CLASWP( 1, RHS, LDA, 1, N-1, JPIV, -1 ) 00130 RETURN 00131 * 00132 * End of CGESC2 00133 * 00134 END
1.7.2