SRC/cgetrf2.f

   1 *> \brief \b CGETRF2
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *  Definition:
   9 *  ===========
  10 *
  11 *       RECURSIVE SUBROUTINE CGETRF2( M, N, A, LDA, IPIV, INFO )
  12 *
  13 *       .. Scalar Arguments ..
  14 *       INTEGER            INFO, LDA, M, N
  15 *       ..
  16 *       .. Array Arguments ..
  17 *       INTEGER            IPIV( * )
  18 *       COMPLEX            A( LDA, * )
  19 *       ..
  20 *
  21 *
  22 *> \par Purpose:
  23 *  =============
  24 *>
  25 *> \verbatim
  26 *>
  27 *> CGETRF2 computes an LU factorization of a general M-by-N matrix A
  28 *> using partial pivoting with row interchanges.
  29 *>
  30 *> The factorization has the form
  31 *>    A = P * L * U
  32 *> where P is a permutation matrix, L is lower triangular with unit
  33 *> diagonal elements (lower trapezoidal if m > n), and U is upper
  34 *> triangular (upper trapezoidal if m < n).
  35 *>
  36 *> This is the recursive version of the algorithm. It divides
  37 *> the matrix into four submatrices:
  38 *>
  39 *>        [  A11 | A12  ]  where A11 is n1 by n1 and A22 is n2 by n2
  40 *>    A = [ -----|----- ]  with n1 = min(m,n)/2
  41 *>        [  A21 | A22  ]       n2 = n-n1
  42 *>
  43 *>                                       [ A11 ]
  44 *> The subroutine calls itself to factor [ --- ],
  45 *>                                       [ A12 ]
  46 *>                 [ A12 ]
  47 *> do the swaps on [ --- ], solve A12, update A22,
  48 *>                 [ A22 ]
  49 *>
  50 *> then calls itself to factor A22 and do the swaps on A21.
  51 *>
  52 *> \endverbatim
  53 *
  54 *  Arguments:
  55 *  ==========
  56 *
  57 *> \param[in] M
  58 *> \verbatim
  59 *>          M is INTEGER
  60 *>          The number of rows of the matrix A.  M >= 0.
  61 *> \endverbatim
  62 *>
  63 *> \param[in] N
  64 *> \verbatim
  65 *>          N is INTEGER
  66 *>          The number of columns of the matrix A.  N >= 0.
  67 *> \endverbatim
  68 *>
  69 *> \param[in,out] A
  70 *> \verbatim
  71 *>          A is COMPLEX array, dimension (LDA,N)
  72 *>          On entry, the M-by-N matrix to be factored.
  73 *>          On exit, the factors L and U from the factorization
  74 *>          A = P*L*U; the unit diagonal elements of L are not stored.
  75 *> \endverbatim
  76 *>
  77 *> \param[in] LDA
  78 *> \verbatim
  79 *>          LDA is INTEGER
  80 *>          The leading dimension of the array A.  LDA >= max(1,M).
  81 *> \endverbatim
  82 *>
  83 *> \param[out] IPIV
  84 *> \verbatim
  85 *>          IPIV is INTEGER array, dimension (min(M,N))
  86 *>          The pivot indices; for 1 <= i <= min(M,N), row i of the
  87 *>          matrix was interchanged with row IPIV(i).
  88 *> \endverbatim
  89 *>
  90 *> \param[out] INFO
  91 *> \verbatim
  92 *>          INFO is INTEGER
  93 *>          = 0:  successful exit
  94 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
  95 *>          > 0:  if INFO = i, U(i,i) is exactly zero. The factorization
  96 *>                has been completed, but the factor U is exactly
  97 *>                singular, and division by zero will occur if it is used
  98 *>                to solve a system of equations.
  99 *> \endverbatim
 100 *
 101 *  Authors:
 102 *  ========
 103 *
 104 *> \author Univ. of Tennessee
 105 *> \author Univ. of California Berkeley
 106 *> \author Univ. of Colorado Denver
 107 *> \author NAG Ltd.
 108 *
 109 *> \date June 2016
 110 *
 111 *> \ingroup complexGEcomputational
 112 *
 113 *  =====================================================================
 114       RECURSIVE SUBROUTINE CGETRF2( M, N, A, LDA, IPIV, INFO )
 115 *
 116 *  -- LAPACK computational routine (version 3.6.1) --
 117 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 118 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 119 *     June 2016
 120 *
 121 *     .. Scalar Arguments ..
 122       INTEGER            INFO, LDA, M, N
 123 *     ..
 124 *     .. Array Arguments ..
 125       INTEGER            IPIV( * )
 126       COMPLEX            A( LDA, * )
 127 *     ..
 128 *
 129 *  =====================================================================
 130 *
 131 *     .. Parameters ..
 132       COMPLEX            ONE, ZERO
 133       PARAMETER          ( ONE = ( 1.0E+0, 0.0E+0 ),
 134      $                     ZERO = ( 0.0E+0, 0.0E+0 ) )
 135 *     ..
 136 *     .. Local Scalars ..
 137       REAL               SFMIN
 138       COMPLEX            TEMP
 139       INTEGER            I, IINFO, N1, N2
 140 *     ..
 141 *     .. External Functions ..
 142       REAL               SLAMCH
 143       INTEGER            ICAMAX
 144       EXTERNAL           SLAMCH, ICAMAX
 145 *     ..
 146 *     .. External Subroutines ..
 147       EXTERNAL           CGEMM, CSCAL, CLASWP, CTRSM, XERBLA
 148 *     ..
 149 *     .. Intrinsic Functions ..
 150       INTRINSIC          MAX, MIN
 151 *     ..
 152 *     .. Executable Statements ..
 153 *
 154 *     Test the input parameters
 155 *
 156       INFO = 0
 157       IF( M.LT.0 ) THEN
 158          INFO = -1
 159       ELSE IF( N.LT.0 ) THEN
 160          INFO = -2
 161       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
 162          INFO = -4
 163       END IF
 164       IF( INFO.NE.0 ) THEN
 165          CALL XERBLA( 'CGETRF2', -INFO )
 166          RETURN
 167       END IF
 168 *
 169 *     Quick return if possible
 170 *
 171       IF( M.EQ.0 .OR. N.EQ.0 )
 172      $   RETURN
 173
 174       IF ( M.EQ.1 ) THEN
 175 *
 176 *        Use unblocked code for one row case
 177 *        Just need to handle IPIV and INFO
 178 *
 179          IPIV( 1 ) = 1
 180          IF ( A(1,1).EQ.ZERO )
 181      $      INFO = 1
 182 *
 183       ELSE IF( N.EQ.1 ) THEN
 184 *
 185 *        Use unblocked code for one column case
 186 *
 187 *
 188 *        Compute machine safe minimum
 189 *
 190          SFMIN = SLAMCH('S')
 191 *
 192 *        Find pivot and test for singularity
 193 *
 194          I = ICAMAX( M, A( 1, 1 ), 1 )
 195          IPIV( 1 ) = I
 196          IF( A( I, 1 ).NE.ZERO ) THEN
 197 *
 198 *           Apply the interchange
 199 *
 200             IF( I.NE.1 ) THEN
 201                TEMP = A( 1, 1 )
 202                A( 1, 1 ) = A( I, 1 )
 203                A( I, 1 ) = TEMP
 204             END IF
 205 *
 206 *           Compute elements 2:M of the column
 207 *
 208             IF( ABS(A( 1, 1 )) .GE. SFMIN ) THEN
 209                CALL CSCAL( M-1, ONE / A( 1, 1 ), A( 2, 1 ), 1 )
 210             ELSE
 211                DO 10 I = 1, M-1
 212                   A( 1+I, 1 ) = A( 1+I, 1 ) / A( 1, 1 )
 213    10          CONTINUE
 214             END IF
 215 *
 216          ELSE
 217             INFO = 1
 218          END IF
 219 *
 220       ELSE
 221 *
 222 *        Use recursive code
 223 *
 224          N1 = MIN( M, N ) / 2
 225          N2 = N-N1
 226 *
 227 *               [ A11 ]
 228 *        Factor [ --- ]
 229 *               [ A21 ]
 230 *
 231          CALL CGETRF2( M, N1, A, LDA, IPIV, IINFO )
 232
 233          IF ( INFO.EQ.0 .AND. IINFO.GT.0 )
 234      $      INFO = IINFO
 235 *
 236 *                              [ A12 ]
 237 *        Apply interchanges to [ --- ]
 238 *                              [ A22 ]
 239 *
 240          CALL CLASWP( N2, A( 1, N1+1 ), LDA, 1, N1, IPIV, 1 )
 241 *
 242 *        Solve A12
 243 *
 244          CALL CTRSM( 'L', 'L', 'N', 'U', N1, N2, ONE, A, LDA,
 245      $               A( 1, N1+1 ), LDA )
 246 *
 247 *        Update A22
 248 *
 249          CALL CGEMM( 'N', 'N', M-N1, N2, N1, -ONE, A( N1+1, 1 ), LDA,
 250      $               A( 1, N1+1 ), LDA, ONE, A( N1+1, N1+1 ), LDA )
 251 *
 252 *        Factor A22
 253 *
 254          CALL CGETRF2( M-N1, N2, A( N1+1, N1+1 ), LDA, IPIV( N1+1 ),
 255      $                 IINFO )
 256 *
 257 *        Adjust INFO and the pivot indices
 258 *
 259          IF ( INFO.EQ.0 .AND. IINFO.GT.0 )
 260      $      INFO = IINFO + N1
 261          DO 20 I = N1+1, MIN( M, N )
 262             IPIV( I ) = IPIV( I ) + N1
 263    20    CONTINUE
 264 *
 265 *        Apply interchanges to A21
 266 *
 267          CALL CLASWP( N1, A( 1, 1 ), LDA, N1+1, MIN( M, N), IPIV, 1 )
 268 *
 269       END IF
 270       RETURN
 271 *
 272 *     End of CGETRF2
 273 *
 274       END