SRC/cgesc2.f

   1 *> \brief \b CGESC2 solves a system of linear equations using the LU factorization with complete pivoting computed by sgetc2.
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download CGESC2 + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cgesc2.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cgesc2.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cgesc2.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE CGESC2( N, A, LDA, RHS, IPIV, JPIV, SCALE )
  22 *
  23 *       .. Scalar Arguments ..
  24 *       INTEGER            LDA, N
  25 *       REAL               SCALE
  26 *       ..
  27 *       .. Array Arguments ..
  28 *       INTEGER            IPIV( * ), JPIV( * )
  29 *       COMPLEX            A( LDA, * ), RHS( * )
  30 *       ..
  31 *
  32 *
  33 *> \par Purpose:
  34 *  =============
  35 *>
  36 *> \verbatim
  37 *>
  38 *> CGESC2 solves a system of linear equations
  39 *>
  40 *>           A * X = scale* RHS
  41 *>
  42 *> with a general N-by-N matrix A using the LU factorization with
  43 *> complete pivoting computed by CGETC2.
  44 *>
  45 *> \endverbatim
  46 *
  47 *  Arguments:
  48 *  ==========
  49 *
  50 *> \param[in] N
  51 *> \verbatim
  52 *>          N is INTEGER
  53 *>          The number of columns of the matrix A.
  54 *> \endverbatim
  55 *>
  56 *> \param[in] A
  57 *> \verbatim
  58 *>          A is COMPLEX array, dimension (LDA, N)
  59 *>          On entry, the  LU part of the factorization of the n-by-n
  60 *>          matrix A computed by CGETC2:  A = P * L * U * Q
  61 *> \endverbatim
  62 *>
  63 *> \param[in] LDA
  64 *> \verbatim
  65 *>          LDA is INTEGER
  66 *>          The leading dimension of the array A.  LDA >= max(1, N).
  67 *> \endverbatim
  68 *>
  69 *> \param[in,out] RHS
  70 *> \verbatim
  71 *>          RHS is COMPLEX array, dimension N.
  72 *>          On entry, the right hand side vector b.
  73 *>          On exit, the solution vector X.
  74 *> \endverbatim
  75 *>
  76 *> \param[in] IPIV
  77 *> \verbatim
  78 *>          IPIV is INTEGER array, dimension (N).
  79 *>          The pivot indices; for 1 <= i <= N, row i of the
  80 *>          matrix has been interchanged with row IPIV(i).
  81 *> \endverbatim
  82 *>
  83 *> \param[in] JPIV
  84 *> \verbatim
  85 *>          JPIV is INTEGER array, dimension (N).
  86 *>          The pivot indices; for 1 <= j <= N, column j of the
  87 *>          matrix has been interchanged with column JPIV(j).
  88 *> \endverbatim
  89 *>
  90 *> \param[out] SCALE
  91 *> \verbatim
  92 *>          SCALE is REAL
  93 *>           On exit, SCALE contains the scale factor. SCALE is chosen
  94 *>           0 <= SCALE <= 1 to prevent owerflow in the solution.
  95 *> \endverbatim
  96 *
  97 *  Authors:
  98 *  ========
  99 *
 100 *> \author Univ. of Tennessee
 101 *> \author Univ. of California Berkeley
 102 *> \author Univ. of Colorado Denver
 103 *> \author NAG Ltd.
 104 *
 105 *> \date September 2012
 106 *
 107 *> \ingroup complexGEauxiliary
 108 *
 109 *> \par Contributors:
 110 *  ==================
 111 *>
 112 *>     Bo Kagstrom and Peter Poromaa, Department of Computing Science,
 113 *>     Umea University, S-901 87 Umea, Sweden.
 114 *
 115 *  =====================================================================
 116       SUBROUTINE CGESC2( N, A, LDA, RHS, IPIV, JPIV, SCALE )
 117 *
 118 *  -- LAPACK auxiliary routine (version 3.4.2) --
 119 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 120 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 121 *     September 2012
 122 *
 123 *     .. Scalar Arguments ..
 124       INTEGER            LDA, N
 125       REAL               SCALE
 126 *     ..
 127 *     .. Array Arguments ..
 128       INTEGER            IPIV( * ), JPIV( * )
 129       COMPLEX            A( LDA, * ), RHS( * )
 130 *     ..
 131 *
 132 *  =====================================================================
 133 *
 134 *     .. Parameters ..
 135       REAL               ZERO, ONE, TWO
 136       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0, TWO = 2.0E+0 )
 137 *     ..
 138 *     .. Local Scalars ..
 139       INTEGER            I, J
 140       REAL               BIGNUM, EPS, SMLNUM
 141       COMPLEX            TEMP
 142 *     ..
 143 *     .. External Subroutines ..
 144       EXTERNAL           CLASWP, CSCAL, SLABAD
 145 *     ..
 146 *     .. External Functions ..
 147       INTEGER            ICAMAX
 148       REAL               SLAMCH
 149       EXTERNAL           ICAMAX, SLAMCH
 150 *     ..
 151 *     .. Intrinsic Functions ..
 152       INTRINSIC          ABS, CMPLX, REAL
 153 *     ..
 154 *     .. Executable Statements ..
 155 *
 156 *     Set constant to control overflow
 157 *
 158       EPS = SLAMCH( 'P' )
 159       SMLNUM = SLAMCH( 'S' ) / EPS
 160       BIGNUM = ONE / SMLNUM
 161       CALL SLABAD( SMLNUM, BIGNUM )
 162 *
 163 *     Apply permutations IPIV to RHS
 164 *
 165       CALL CLASWP( 1, RHS, LDA, 1, N-1, IPIV, 1 )
 166 *
 167 *     Solve for L part
 168 *
 169       DO 20 I = 1, N - 1
 170          DO 10 J = I + 1, N
 171             RHS( J ) = RHS( J ) - A( J, I )*RHS( I )
 172    10    CONTINUE
 173    20 CONTINUE
 174 *
 175 *     Solve for U part
 176 *
 177       SCALE = ONE
 178 *
 179 *     Check for scaling
 180 *
 181       I = ICAMAX( N, RHS, 1 )
 182       IF( TWO*SMLNUM*ABS( RHS( I ) ).GT.ABS( A( N, N ) ) ) THEN
 183          TEMP = CMPLX( ONE / TWO, ZERO ) / ABS( RHS( I ) )
 184          CALL CSCAL( N, TEMP, RHS( 1 ), 1 )
 185          SCALE = SCALE*REAL( TEMP )
 186       END IF
 187       DO 40 I = N, 1, -1
 188          TEMP = CMPLX( ONE, ZERO ) / A( I, I )
 189          RHS( I ) = RHS( I )*TEMP
 190          DO 30 J = I + 1, N
 191             RHS( I ) = RHS( I ) - RHS( J )*( A( I, J )*TEMP )
 192    30    CONTINUE
 193    40 CONTINUE
 194 *
 195 *     Apply permutations JPIV to the solution (RHS)
 196 *
 197       CALL CLASWP( 1, RHS, LDA, 1, N-1, JPIV, -1 )
 198       RETURN
 199 *
 200 *     End of CGESC2
 201 *
 202       END