Télécharger dgetrf2.eso

Retour à la liste

Numérotation des lignes :

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

© Cast3M 2003 - Tous droits réservés.
Mentions légales