LAPACK 3.3.1 Linear Algebra PACKage

zgetc2.f

Go to the documentation of this file.
```00001       SUBROUTINE ZGETC2( N, A, LDA, IPIV, JPIV, INFO )
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            INFO, LDA, N
00010 *     ..
00011 *     .. Array Arguments ..
00012       INTEGER            IPIV( * ), JPIV( * )
00013       COMPLEX*16         A( LDA, * )
00014 *     ..
00015 *
00016 *  Purpose
00017 *  =======
00018 *
00019 *  ZGETC2 computes an LU factorization, using complete pivoting, of the
00020 *  n-by-n matrix A. The factorization has the form A = P * L * U * Q,
00021 *  where P and Q are permutation matrices, L is lower triangular with
00022 *  unit diagonal elements and U is upper triangular.
00023 *
00024 *  This is a level 1 BLAS version of the algorithm.
00025 *
00026 *  Arguments
00027 *  =========
00028 *
00029 *  N       (input) INTEGER
00030 *          The order of the matrix A. N >= 0.
00031 *
00032 *  A       (input/output) COMPLEX*16 array, dimension (LDA, N)
00033 *          On entry, the n-by-n matrix to be factored.
00034 *          On exit, the factors L and U from the factorization
00035 *          A = P*L*U*Q; the unit diagonal elements of L are not stored.
00036 *          If U(k, k) appears to be less than SMIN, U(k, k) is given the
00037 *          value of SMIN, giving a nonsingular perturbed system.
00038 *
00039 *  LDA     (input) INTEGER
00040 *          The leading dimension of the array A.  LDA >= max(1, N).
00041 *
00042 *  IPIV    (output) INTEGER array, dimension (N).
00043 *          The pivot indices; for 1 <= i <= N, row i of the
00044 *          matrix has been interchanged with row IPIV(i).
00045 *
00046 *  JPIV    (output) INTEGER array, dimension (N).
00047 *          The pivot indices; for 1 <= j <= N, column j of the
00048 *          matrix has been interchanged with column JPIV(j).
00049 *
00050 *  INFO    (output) INTEGER
00051 *           = 0: successful exit
00052 *           > 0: if INFO = k, U(k, k) is likely to produce overflow if
00053 *                one tries to solve for x in Ax = b. So U is perturbed
00054 *                to avoid the overflow.
00055 *
00056 *  Further Details
00057 *  ===============
00058 *
00059 *  Based on contributions by
00060 *     Bo Kagstrom and Peter Poromaa, Department of Computing Science,
00061 *     Umea University, S-901 87 Umea, Sweden.
00062 *
00063 *  =====================================================================
00064 *
00065 *     .. Parameters ..
00066       DOUBLE PRECISION   ZERO, ONE
00067       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
00068 *     ..
00069 *     .. Local Scalars ..
00070       INTEGER            I, IP, IPV, J, JP, JPV
00071       DOUBLE PRECISION   BIGNUM, EPS, SMIN, SMLNUM, XMAX
00072 *     ..
00073 *     .. External Subroutines ..
00074       EXTERNAL           ZGERU, ZSWAP
00075 *     ..
00076 *     .. External Functions ..
00077       DOUBLE PRECISION   DLAMCH
00078       EXTERNAL           DLAMCH
00079 *     ..
00080 *     .. Intrinsic Functions ..
00081       INTRINSIC          ABS, DCMPLX, MAX
00082 *     ..
00083 *     .. Executable Statements ..
00084 *
00085 *     Set constants to control overflow
00086 *
00087       INFO = 0
00088       EPS = DLAMCH( 'P' )
00089       SMLNUM = DLAMCH( 'S' ) / EPS
00090       BIGNUM = ONE / SMLNUM
00091       CALL DLABAD( SMLNUM, BIGNUM )
00092 *
00093 *     Factorize A using complete pivoting.
00094 *     Set pivots less than SMIN to SMIN
00095 *
00096       DO 40 I = 1, N - 1
00097 *
00098 *        Find max element in matrix A
00099 *
00100          XMAX = ZERO
00101          DO 20 IP = I, N
00102             DO 10 JP = I, N
00103                IF( ABS( A( IP, JP ) ).GE.XMAX ) THEN
00104                   XMAX = ABS( A( IP, JP ) )
00105                   IPV = IP
00106                   JPV = JP
00107                END IF
00108    10       CONTINUE
00109    20    CONTINUE
00110          IF( I.EQ.1 )
00111      \$      SMIN = MAX( EPS*XMAX, SMLNUM )
00112 *
00113 *        Swap rows
00114 *
00115          IF( IPV.NE.I )
00116      \$      CALL ZSWAP( N, A( IPV, 1 ), LDA, A( I, 1 ), LDA )
00117          IPIV( I ) = IPV
00118 *
00119 *        Swap columns
00120 *
00121          IF( JPV.NE.I )
00122      \$      CALL ZSWAP( N, A( 1, JPV ), 1, A( 1, I ), 1 )
00123          JPIV( I ) = JPV
00124 *
00125 *        Check for singularity
00126 *
00127          IF( ABS( A( I, I ) ).LT.SMIN ) THEN
00128             INFO = I
00129             A( I, I ) = DCMPLX( SMIN, ZERO )
00130          END IF
00131          DO 30 J = I + 1, N
00132             A( J, I ) = A( J, I ) / A( I, I )
00133    30    CONTINUE
00134          CALL ZGERU( N-I, N-I, -DCMPLX( ONE ), A( I+1, I ), 1,
00135      \$               A( I, I+1 ), LDA, A( I+1, I+1 ), LDA )
00136    40 CONTINUE
00137 *
00138       IF( ABS( A( N, N ) ).LT.SMIN ) THEN
00139          INFO = N
00140          A( N, N ) = DCMPLX( SMIN, ZERO )
00141       END IF
00142       RETURN
00143 *
00144 *     End of ZGETC2
00145 *
00146       END
```