SRC/sgesc2.f

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