SRC/sgetc2.f

   1 *> \brief \b SGETC2 computes the LU factorization with complete pivoting of the general n-by-n matrix.
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download SGETC2 + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sgetc2.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sgetc2.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sgetc2.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE SGETC2( N, A, LDA, IPIV, JPIV, INFO )
  22 *
  23 *       .. Scalar Arguments ..
  24 *       INTEGER            INFO, LDA, N
  25 *       ..
  26 *       .. Array Arguments ..
  27 *       INTEGER            IPIV( * ), JPIV( * )
  28 *       REAL               A( LDA, * )
  29 *       ..
  30 *
  31 *
  32 *> \par Purpose:
  33 *  =============
  34 *>
  35 *> \verbatim
  36 *>
  37 *> SGETC2 computes an LU factorization with complete pivoting of the
  38 *> n-by-n matrix A. The factorization has the form A = P * L * U * Q,
  39 *> where P and Q are permutation matrices, L is lower triangular with
  40 *> unit diagonal elements and U is upper triangular.
  41 *>
  42 *> This is the Level 2 BLAS algorithm.
  43 *> \endverbatim
  44 *
  45 *  Arguments:
  46 *  ==========
  47 *
  48 *> \param[in] N
  49 *> \verbatim
  50 *>          N is INTEGER
  51 *>          The order of the matrix A. N >= 0.
  52 *> \endverbatim
  53 *>
  54 *> \param[in,out] A
  55 *> \verbatim
  56 *>          A is REAL array, dimension (LDA, N)
  57 *>          On entry, the n-by-n matrix A to be factored.
  58 *>          On exit, the factors L and U from the factorization
  59 *>          A = P*L*U*Q; the unit diagonal elements of L are not stored.
  60 *>          If U(k, k) appears to be less than SMIN, U(k, k) is given the
  61 *>          value of SMIN, i.e., giving a nonsingular perturbed system.
  62 *> \endverbatim
  63 *>
  64 *> \param[in] LDA
  65 *> \verbatim
  66 *>          LDA is INTEGER
  67 *>          The leading dimension of the array A.  LDA >= max(1,N).
  68 *> \endverbatim
  69 *>
  70 *> \param[out] IPIV
  71 *> \verbatim
  72 *>          IPIV is INTEGER array, dimension(N).
  73 *>          The pivot indices; for 1 <= i <= N, row i of the
  74 *>          matrix has been interchanged with row IPIV(i).
  75 *> \endverbatim
  76 *>
  77 *> \param[out] JPIV
  78 *> \verbatim
  79 *>          JPIV is INTEGER array, dimension(N).
  80 *>          The pivot indices; for 1 <= j <= N, column j of the
  81 *>          matrix has been interchanged with column JPIV(j).
  82 *> \endverbatim
  83 *>
  84 *> \param[out] INFO
  85 *> \verbatim
  86 *>          INFO is INTEGER
  87 *>           = 0: successful exit
  88 *>           > 0: if INFO = k, U(k, k) is likely to produce owerflow if
  89 *>                we try to solve for x in Ax = b. So U is perturbed to
  90 *>                avoid the overflow.
  91 *> \endverbatim
  92 *
  93 *  Authors:
  94 *  ========
  95 *
  96 *> \author Univ. of Tennessee
  97 *> \author Univ. of California Berkeley
  98 *> \author Univ. of Colorado Denver
  99 *> \author NAG Ltd.
 100 *
 101 *> \date June 2016
 102 *
 103 *> \ingroup realGEauxiliary
 104 *
 105 *> \par Contributors:
 106 *  ==================
 107 *>
 108 *>     Bo Kagstrom and Peter Poromaa, Department of Computing Science,
 109 *>     Umea University, S-901 87 Umea, Sweden.
 110 *
 111 *  =====================================================================
 112       SUBROUTINE SGETC2( N, A, LDA, IPIV, JPIV, INFO )
 113 *
 114 *  -- LAPACK auxiliary routine (version 3.6.1) --
 115 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 116 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 117 *     June 2016
 118 *
 119 *     .. Scalar Arguments ..
 120       INTEGER            INFO, LDA, N
 121 *     ..
 122 *     .. Array Arguments ..
 123       INTEGER            IPIV( * ), JPIV( * )
 124       REAL               A( LDA, * )
 125 *     ..
 126 *
 127 *  =====================================================================
 128 *
 129 *     .. Parameters ..
 130       REAL               ZERO, ONE
 131       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
 132 *     ..
 133 *     .. Local Scalars ..
 134       INTEGER            I, IP, IPV, J, JP, JPV
 135       REAL               BIGNUM, EPS, SMIN, SMLNUM, XMAX
 136 *     ..
 137 *     .. External Subroutines ..
 138       EXTERNAL           SGER, SLABAD, SSWAP
 139 *     ..
 140 *     .. External Functions ..
 141       REAL               SLAMCH
 142       EXTERNAL           SLAMCH
 143 *     ..
 144 *     .. Intrinsic Functions ..
 145       INTRINSIC          ABS, MAX
 146 *     ..
 147 *     .. Executable Statements ..
 148 *
 149       INFO = 0
 150 *
 151 *     Quick return if possible
 152 *
 153       IF( N.EQ.0 )
 154      $   RETURN
 155 *
 156 *     Set constants to control overflow
 157 *
 158       EPS = SLAMCH( 'P' )
 159       SMLNUM = SLAMCH( 'S' ) / EPS
 160       BIGNUM = ONE / SMLNUM
 161       CALL SLABAD( SMLNUM, BIGNUM )
 162 *
 163 *     Handle the case N=1 by itself
 164 *
 165       IF( N.EQ.1 ) THEN
 166          IPIV( 1 ) = 1
 167          JPIV( 1 ) = 1
 168          IF( ABS( A( 1, 1 ) ).LT.SMLNUM ) THEN
 169             INFO = 1
 170             A( 1, 1 ) = SMLNUM
 171          END IF
 172          RETURN
 173       END IF
 174 *
 175 *     Factorize A using complete pivoting.
 176 *     Set pivots less than SMIN to SMIN.
 177 *
 178       DO 40 I = 1, N - 1
 179 *
 180 *        Find max element in matrix A
 181 *
 182          XMAX = ZERO
 183          DO 20 IP = I, N
 184             DO 10 JP = I, N
 185                IF( ABS( A( IP, JP ) ).GE.XMAX ) THEN
 186                   XMAX = ABS( A( IP, JP ) )
 187                   IPV = IP
 188                   JPV = JP
 189                END IF
 190    10       CONTINUE
 191    20    CONTINUE
 192          IF( I.EQ.1 )
 193      $      SMIN = MAX( EPS*XMAX, SMLNUM )
 194 *
 195 *        Swap rows
 196 *
 197          IF( IPV.NE.I )
 198      $      CALL SSWAP( N, A( IPV, 1 ), LDA, A( I, 1 ), LDA )
 199          IPIV( I ) = IPV
 200 *
 201 *        Swap columns
 202 *
 203          IF( JPV.NE.I )
 204      $      CALL SSWAP( N, A( 1, JPV ), 1, A( 1, I ), 1 )
 205          JPIV( I ) = JPV
 206 *
 207 *        Check for singularity
 208 *
 209          IF( ABS( A( I, I ) ).LT.SMIN ) THEN
 210             INFO = I
 211             A( I, I ) = SMIN
 212          END IF
 213          DO 30 J = I + 1, N
 214             A( J, I ) = A( J, I ) / A( I, I )
 215    30    CONTINUE
 216          CALL SGER( N-I, N-I, -ONE, A( I+1, I ), 1, A( I, I+1 ), LDA,
 217      $              A( I+1, I+1 ), LDA )
 218    40 CONTINUE
 219 *
 220       IF( ABS( A( N, N ) ).LT.SMIN ) THEN
 221          INFO = N
 222          A( N, N ) = SMIN
 223       END IF
 224 *
 225 *     Set last pivots to N
 226 *
 227       IPIV( N ) = N
 228       JPIV( N ) = N
 229 *
 230       RETURN
 231 *
 232 *     End of SGETC2
 233 *
 234       END