dgetc2.f

*> \brief \b DGETC2 computes the LU factorization with complete pivoting of the general n-by-n matrix.
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*> \htmlonly
*> Download DGETC2 + dependencies
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgetc2.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgetc2.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgetc2.f">
*> [TXT]</a>
*> \endhtmlonly
*
*  Definition:
*  ===========
*
*       SUBROUTINE DGETC2( N, A, LDA, IPIV, JPIV, INFO )
*
*       .. Scalar Arguments ..
*       INTEGER            INFO, LDA, N
*       ..
*       .. Array Arguments ..
*       INTEGER            IPIV( * ), JPIV( * )
*       DOUBLE PRECISION   A( LDA, * )
*       ..
*
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> DGETC2 computes an LU factorization with complete pivoting of the
*> n-by-n matrix A. The factorization has the form A = P * L * U * Q,
*> where P and Q are permutation matrices, L is lower triangular with
*> unit diagonal elements and U is upper triangular.
*>
*> This is the Level 2 BLAS algorithm.
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] N
*> \verbatim
*>          N is INTEGER
*>          The order of the matrix A. N >= 0.
*> \endverbatim
*>
*> \param[in,out] A
*> \verbatim
*>          A is DOUBLE PRECISION array, dimension (LDA, N)
*>          On entry, the n-by-n matrix A to be factored.
*>          On exit, the factors L and U from the factorization
*>          A = P*L*U*Q; the unit diagonal elements of L are not stored.
*>          If U(k, k) appears to be less than SMIN, U(k, k) is given the
*>          value of SMIN, i.e., giving a nonsingular perturbed system.
*> \endverbatim
*>
*> \param[in] LDA
*> \verbatim
*>          LDA is INTEGER
*>          The leading dimension of the array A.  LDA >= max(1,N).
*> \endverbatim
*>
*> \param[out] IPIV
*> \verbatim
*>          IPIV is INTEGER array, dimension(N).
*>          The pivot indices; for 1 <= i <= N, row i of the
*>          matrix has been interchanged with row IPIV(i).
*> \endverbatim
*>
*> \param[out] JPIV
*> \verbatim
*>          JPIV is INTEGER array, dimension(N).
*>          The pivot indices; for 1 <= j <= N, column j of the
*>          matrix has been interchanged with column JPIV(j).
*> \endverbatim
*>
*> \param[out] INFO
*> \verbatim
*>          INFO is INTEGER
*>           = 0: successful exit
*>           > 0: if INFO = k, U(k, k) is likely to produce overflow if
*>                we try to solve for x in Ax = b. So U is perturbed to
*>                avoid the overflow.
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup getc2
*
*> \par Contributors:
*  ==================
*>
*>     Bo Kagstrom and Peter Poromaa, Department of Computing Science,
*>     Umea University, S-901 87 Umea, Sweden.
*
*  =====================================================================
      SUBROUTINE dgetc2( N, A, LDA, IPIV, JPIV, INFO )
*
*  -- LAPACK auxiliary routine --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
*     .. Scalar Arguments ..
      INTEGER            INFO, LDA, N
*     ..
*     .. Array Arguments ..
      INTEGER            IPIV( * ), JPIV( * )
      DOUBLE PRECISION   A( lda, * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   ZERO, ONE
      parameter( zero = 0.0d+0, one = 1.0d+0 )
*     ..
*     .. Local Scalars ..
      INTEGER            I, IP, IPV, J, JP, JPV
      DOUBLE PRECISION   BIGNUM, EPS, SMIN, SMLNUM, XMAX
*     ..
*     .. External Subroutines ..
      EXTERNAL           dger, dswap, dlabad
*     ..
*     .. External Functions ..
      DOUBLE PRECISION   DLAMCH
      EXTERNAL           dlamch
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, max
*     ..
*     .. Executable Statements ..
*
      info = 0
*
*     Quick return if possible
*
      IF( n.EQ.0 )
     $   RETURN
*
*     Set constants to control overflow
*
      eps = dlamch( 'P' )
      smlnum = dlamch( 'S' ) / eps
      bignum = one / smlnum
      CALL dlabad( smlnum, bignum )
*
*     Handle the case N=1 by itself
*
      IF( n.EQ.1 ) THEN
         ipiv( 1 ) = 1
         jpiv( 1 ) = 1
         IF( abs( a( 1, 1 ) ).LT.smlnum ) THEN
            info = 1
            a( 1, 1 ) = smlnum
         END IF
         RETURN
      END IF
*
*     Factorize A using complete pivoting.
*     Set pivots less than SMIN to SMIN.
*
      DO 40 i = 1, n - 1
*
*        Find max element in matrix A
*
         xmax = zero
         DO 20 ip = i, n
            DO 10 jp = i, n
               IF( abs( a( ip, jp ) ).GE.xmax ) THEN
                  xmax = abs( a( ip, jp ) )
                  ipv = ip
                  jpv = jp
               END IF
   10       CONTINUE
   20    CONTINUE
         IF( i.EQ.1 )
     $      smin = max( eps*xmax, smlnum )
*
*        Swap rows
*
         IF( ipv.NE.i )
     $      CALL dswap( n, a( ipv, 1 ), lda, a( i, 1 ), lda )
         ipiv( i ) = ipv
*
*        Swap columns
*
         IF( jpv.NE.i )
     $      CALL dswap( n, a( 1, jpv ), 1, a( 1, i ), 1 )
         jpiv( i ) = jpv
*
*        Check for singularity
*
         IF( abs( a( i, i ) ).LT.smin ) THEN
            info = i
            a( i, i ) = smin
         END IF
         DO 30 j = i + 1, n
            a( j, i ) = a( j, i ) / a( i, i )
   30    CONTINUE
         CALL dger( n-i, n-i, -one, a( i+1, i ), 1, a( i, i+1 ), lda,
     $              a( i+1, i+1 ), lda )
   40 CONTINUE
*
      IF( abs( a( n, n ) ).LT.smin ) THEN
         info = n
         a( n, n ) = smin
      END IF
*
*     Set last pivots to N
*
      ipiv( n ) = n
      jpiv( n ) = n
*
      RETURN
*
*     End of DGETC2
*
      END