src/LAPACK/zgetrf.f

   1       SUBROUTINE ZGETRF( M, N, A, LDA, IPIV, INFO )
   2 *
   3 *  -- LAPACK routine (version 3.1) --
   4 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
   5 *     November 2006
   6 *
   7 *     .. Scalar Arguments ..
   8       INTEGER            INFO, LDA, M, N
   9 *     ..
  10 *     .. Array Arguments ..
  11       INTEGER            IPIV( * )
  12       COMPLEX*16         A( LDA, * )
  13 *     ..
  14 *
  15 *  Purpose
  16 *  =======
  17 *
  18 *  ZGETRF computes an LU factorization of a general M-by-N matrix A
  19 *  using partial pivoting with row interchanges.
  20 *
  21 *  The factorization has the form
  22 *     A = P * L * U
  23 *  where P is a permutation matrix, L is lower triangular with unit
  24 *  diagonal elements (lower trapezoidal if m > n), and U is upper
  25 *  triangular (upper trapezoidal if m < n).
  26 *
  27 *  This is the right-looking Level 3 BLAS version of the algorithm.
  28 *
  29 *  Arguments
  30 *  =========
  31 *
  32 *  M       (input) INTEGER
  33 *          The number of rows of the matrix A.  M >= 0.
  34 *
  35 *  N       (input) INTEGER
  36 *          The number of columns of the matrix A.  N >= 0.
  37 *
  38 *  A       (input/output) COMPLEX*16 array, dimension (LDA,N)
  39 *          On entry, the M-by-N matrix to be factored.
  40 *          On exit, the factors L and U from the factorization
  41 *          A = P*L*U; the unit diagonal elements of L are not stored.
  42 *
  43 *  LDA     (input) INTEGER
  44 *          The leading dimension of the array A.  LDA >= max(1,M).
  45 *
  46 *  IPIV    (output) INTEGER array, dimension (min(M,N))
  47 *          The pivot indices; for 1 <= i <= min(M,N), row i of the
  48 *          matrix was interchanged with row IPIV(i).
  49 *
  50 *  INFO    (output) INTEGER
  51 *          = 0:  successful exit
  52 *          < 0:  if INFO = -i, the i-th argument had an illegal value
  53 *          > 0:  if INFO = i, U(i,i) is exactly zero. The factorization
  54 *                has been completed, but the factor U is exactly
  55 *                singular, and division by zero will occur if it is used
  56 *                to solve a system of equations.
  57 *
  58 *  =====================================================================
  59 *
  60 *     .. Parameters ..
  61       COMPLEX*16         ONE
  62       PARAMETER          ( ONE = ( 1.0D+0, 0.0D+0 ) )
  63 *     ..
  64 *     .. Local Scalars ..
  65       INTEGER            I, IINFO, J, JB, NB
  66 *     ..
  67 *     .. External Subroutines ..
  68       EXTERNAL           XERBLA, ZGEMM, ZGETF2, ZLASWP, ZTRSM
  69 *     ..
  70 *     .. External Functions ..
  71       INTEGER            ILAENV
  72       EXTERNAL           ILAENV
  73 *     ..
  74 *     .. Intrinsic Functions ..
  75       INTRINSIC          MAX, MIN
  76 *     ..
  77 *     .. Executable Statements ..
  78 *
  79 *     Test the input parameters.
  80 *
  81       INFO = 0
  82       IF( M.LT.0 ) THEN
  83          INFO = -1
  84       ELSE IF( N.LT.0 ) THEN
  85          INFO = -2
  86       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
  87          INFO = -4
  88       END IF
  89       IF( INFO.NE.0 ) THEN
  90          CALL XERBLA( 'ZGETRF', -INFO )
  91          RETURN
  92       END IF
  93 *
  94 *     Quick return if possible
  95 *
  96       IF( M.EQ.0 .OR. N.EQ.0 )
  97      $   RETURN
  98 *
  99 *     Determine the block size for this environment.
 100 *
 101       NB = ILAENV( 1, 'ZGETRF', ' ', M, N, -1, -1 )
 102       IF( NB.LE.1 .OR. NB.GE.MIN( M, N ) ) THEN
 103 *
 104 *        Use unblocked code.
 105 *
 106          CALL ZGETF2( M, N, A, LDA, IPIV, INFO )
 107       ELSE
 108 *
 109 *        Use blocked code.
 110 *
 111          DO 20 J = 1, MIN( M, N ), NB
 112             JB = MIN( MIN( M, N )-J+1, NB )
 113 *
 114 *           Factor diagonal and subdiagonal blocks and test for exact
 115 *           singularity.
 116 *
 117             CALL ZGETF2( M-J+1, JB, A( J, J ), LDA, IPIV( J ), IINFO )
 118 *
 119 *           Adjust INFO and the pivot indices.
 120 *
 121             IF( INFO.EQ.0 .AND. IINFO.GT.0 )
 122      $         INFO = IINFO + J - 1
 123             DO 10 I = J, MIN( M, J+JB-1 )
 124                IPIV( I ) = J - 1 + IPIV( I )
 125    10       CONTINUE
 126 *
 127 *           Apply interchanges to columns 1:J-1.
 128 *
 129             CALL ZLASWP( J-1, A, LDA, J, J+JB-1, IPIV, 1 )
 130 *
 131             IF( J+JB.LE.N ) THEN
 132 *
 133 *              Apply interchanges to columns J+JB:N.
 134 *
 135                CALL ZLASWP( N-J-JB+1, A( 1, J+JB ), LDA, J, J+JB-1,
 136      $                      IPIV, 1 )
 137 *
 138 *              Compute block row of U.
 139 *
 140                CALL ZTRSM( 'Left', 'Lower', 'No transpose', 'Unit', JB,
 141      $                     N-J-JB+1, ONE, A( J, J ), LDA, A( J, J+JB ),
 142      $                     LDA )
 143                IF( J+JB.LE.M ) THEN
 144 *
 145 *                 Update trailing submatrix.
 146 *
 147                   CALL ZGEMM( 'No transpose', 'No transpose', M-J-JB+1,
 148      $                        N-J-JB+1, JB, -ONE, A( J+JB, J ), LDA,
 149      $                        A( J, J+JB ), LDA, ONE, A( J+JB, J+JB ),
 150      $                        LDA )
 151                END IF
 152             END IF
 153    20    CONTINUE
 154       END IF
 155       RETURN
 156 *
 157 *     End of ZGETRF
 158 *
 159       END