diff options
Diffstat (limited to 'SRC/dgetrf2.f')
-rw-r--r-- | SRC/dgetrf2.f | 272 |
1 files changed, 272 insertions, 0 deletions
diff --git a/SRC/dgetrf2.f b/SRC/dgetrf2.f new file mode 100644 index 00000000..b1871b5d --- /dev/null +++ b/SRC/dgetrf2.f @@ -0,0 +1,272 @@ +*> \brief \b DGETRF2 +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +* Definition: +* =========== +* +* RECURSIVE SUBROUTINE DGETRF2( M, N, A, LDA, IPIV, INFO ) +* +* .. Scalar Arguments .. +* INTEGER INFO, LDA, M, N +* .. +* .. Array Arguments .. +* INTEGER IPIV( * ) +* DOUBLE PRECISION A( LDA, * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> DGETRF2 computes an LU factorization of a general M-by-N matrix A +*> using partial pivoting with row interchanges. +*> +*> The factorization has the form +*> A = P * L * U +*> where P is a permutation matrix, L is lower triangular with unit +*> diagonal elements (lower trapezoidal if m > n), and U is upper +*> triangular (upper trapezoidal if m < n). +*> +*> This is the recursive version of the algorithm. It divides +*> the matrix into four submatrices: +*> +*> [ A11 | A12 ] where A11 is n1 by n1 and A22 is n2 by n2 +*> A = [ -----|----- ] with n1 = min(m,n) +*> [ A21 | A22 ] n2 = n-n1 +*> +*> [ A11 ] +*> The subroutine calls itself to factor [ --- ], +*> [ A12 ] +*> [ A12 ] +*> do the swaps on [ --- ], solve A12, update A22, +*> [ A22 ] +*> +*> then calls itself to factor A22 and do the swaps on A21. +*> +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] M +*> \verbatim +*> M is INTEGER +*> The number of rows of the matrix A. M >= 0. +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The number of columns of the matrix A. N >= 0. +*> \endverbatim +*> +*> \param[in,out] A +*> \verbatim +*> A is DOUBLE PRECISION array, dimension (LDA,N) +*> On entry, the M-by-N matrix to be factored. +*> On exit, the factors L and U from the factorization +*> A = P*L*U; the unit diagonal elements of L are not stored. +*> \endverbatim +*> +*> \param[in] LDA +*> \verbatim +*> LDA is INTEGER +*> The leading dimension of the array A. LDA >= max(1,M). +*> \endverbatim +*> +*> \param[out] IPIV +*> \verbatim +*> IPIV is INTEGER array, dimension (min(M,N)) +*> The pivot indices; for 1 <= i <= min(M,N), row i of the +*> matrix was interchanged with row IPIV(i). +*> \endverbatim +*> +*> \param[out] INFO +*> \verbatim +*> INFO is INTEGER +*> = 0: successful exit +*> < 0: if INFO = -i, the i-th argument had an illegal value +*> > 0: if INFO = i, U(i,i) is exactly zero. The factorization +*> has been completed, but the factor U is exactly +*> singular, and division by zero will occur if it is used +*> to solve a system of equations. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2015 +* +*> \ingroup doubleGEcomputational +* +* ===================================================================== + RECURSIVE SUBROUTINE DGETRF2( M, N, A, LDA, IPIV, INFO ) +* +* -- LAPACK computational routine (version 3.6.0) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2015 +* +* .. Scalar Arguments .. + INTEGER INFO, LDA, M, N +* .. +* .. Array Arguments .. + INTEGER IPIV( * ) + DOUBLE PRECISION A( LDA, * ) +* .. +* +* ===================================================================== +* +* .. Parameters .. + DOUBLE PRECISION ONE, ZERO + PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) +* .. +* .. Local Scalars .. + DOUBLE PRECISION SFMIN, TEMP + INTEGER I, IINFO, N1, N2 +* .. +* .. External Functions .. + DOUBLE PRECISION DLAMCH + INTEGER IDAMAX + EXTERNAL DLAMCH, IDAMAX +* .. +* .. External Subroutines .. + EXTERNAL DGEMM, DSCAL, DLASWP, DTRSM, XERBLA +* .. +* .. Intrinsic Functions .. + INTRINSIC MAX, MIN +* .. +* .. Executable Statements .. +* +* Test the input parameters +* + INFO = 0 + IF( M.LT.0 ) THEN + INFO = -1 + ELSE IF( N.LT.0 ) THEN + INFO = -2 + ELSE IF( LDA.LT.MAX( 1, M ) ) THEN + INFO = -4 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'DGETRF2', -INFO ) + RETURN + END IF +* +* Quick return if possible +* + IF( M.EQ.0 .OR. N.EQ.0 ) + $ RETURN + + IF ( M.EQ.1 ) THEN +* +* Use unblocked code for one row case +* Just need to handle IPIV and INFO +* + IPIV( 1 ) = 1 + IF ( A(1,1).EQ.ZERO ) + $ INFO = 1 +* + ELSE IF( N.EQ.1 ) THEN +* +* Use unblocked code for one column case +* +* +* Compute machine safe minimum +* + SFMIN = DLAMCH('S') +* +* Find pivot and test for singularity +* + I = IDAMAX( M, A( 1, 1 ), 1 ) + IPIV( 1 ) = I + IF( A( I, 1 ).NE.ZERO ) THEN +* +* Apply the interchange +* + IF( I.NE.1 ) THEN + TEMP = A( 1, 1 ) + A( 1, 1 ) = A( I, 1 ) + A( I, 1 ) = TEMP + END IF +* +* Compute elements 2:M of the column +* + IF( ABS(A( 1, 1 )) .GE. SFMIN ) THEN + CALL DSCAL( M-1, ONE / A( 1, 1 ), A( 2, 1 ), 1 ) + ELSE + DO 10 I = 1, M-1 + A( 1+I, 1 ) = A( 1+I, 1 ) / A( 1, 1 ) + 10 CONTINUE + END IF +* + ELSE + INFO = 1 + END IF +* + ELSE +* +* Use recursive code +* + N1 = MIN( M, N ) / 2 + N2 = N-N1 +* +* [ A11 ] +* Factor [ --- ] +* [ A21 ] +* + CALL DGETRF2( M, N1, A, LDA, IPIV, IINFO ) + + IF ( INFO.EQ.0 .AND. IINFO.GT.0 ) + $ INFO = IINFO +* +* [ A12 ] +* Apply interchanges to [ --- ] +* [ A22 ] +* + CALL DLASWP( N2, A( 1, N1+1 ), LDA, 1, N1, IPIV, 1 ) +* +* Solve A12 +* + CALL DTRSM( 'L', 'L', 'N', 'U', N1, N2, ONE, A, LDA, + $ A( 1, N1+1 ), LDA ) +* +* Update A22 +* + CALL DGEMM( 'N', 'N', M-N1, N2, N1, -ONE, A( N1+1, 1 ), LDA, + $ A( 1, N1+1 ), LDA, ONE, A( N1+1, N1+1 ), LDA ) +* +* Factor A22 +* + CALL DGETRF2( M-N1, N2, A( N1+1, N1+1 ), LDA, IPIV( N1+1 ), + $ IINFO ) +* +* Adjust INFO and the pivot indices +* + IF ( INFO.EQ.0 .AND. IINFO.GT.0 ) + $ INFO = IINFO + N1 + DO 20 I = N1+1, MIN( M, N ) + IPIV( I ) = IPIV( I ) + N1 + 20 CONTINUE +* +* Apply interchanges to A21 +* + CALL DLASWP( N1, A( 1, 1 ), LDA, N1+1, MIN( M, N), IPIV, 1 ) +* + END IF + RETURN +* +* End of DGETRF2 +* + END |