You can not select more than 25 topics Topics must start with a chinese character,a letter or number, can include dashes ('-') and can be up to 35 characters long.

zgesc2.f 5.5 kB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202
  1. *> \brief \b ZGESC2 solves a system of linear equations using the LU factorization with complete pivoting computed by sgetc2.
  2. *
  3. * =========== DOCUMENTATION ===========
  4. *
  5. * Online html documentation available at
  6. * http://www.netlib.org/lapack/explore-html/
  7. *
  8. *> \htmlonly
  9. *> Download ZGESC2 + dependencies
  10. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgesc2.f">
  11. *> [TGZ]</a>
  12. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgesc2.f">
  13. *> [ZIP]</a>
  14. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgesc2.f">
  15. *> [TXT]</a>
  16. *> \endhtmlonly
  17. *
  18. * Definition:
  19. * ===========
  20. *
  21. * SUBROUTINE ZGESC2( N, A, LDA, RHS, IPIV, JPIV, SCALE )
  22. *
  23. * .. Scalar Arguments ..
  24. * INTEGER LDA, N
  25. * DOUBLE PRECISION SCALE
  26. * ..
  27. * .. Array Arguments ..
  28. * INTEGER IPIV( * ), JPIV( * )
  29. * COMPLEX*16 A( LDA, * ), RHS( * )
  30. * ..
  31. *
  32. *
  33. *> \par Purpose:
  34. * =============
  35. *>
  36. *> \verbatim
  37. *>
  38. *> ZGESC2 solves a system of linear equations
  39. *>
  40. *> A * X = scale* RHS
  41. *>
  42. *> with a general N-by-N matrix A using the LU factorization with
  43. *> complete pivoting computed by ZGETC2.
  44. *>
  45. *> \endverbatim
  46. *
  47. * Arguments:
  48. * ==========
  49. *
  50. *> \param[in] N
  51. *> \verbatim
  52. *> N is INTEGER
  53. *> The number of columns of the matrix A.
  54. *> \endverbatim
  55. *>
  56. *> \param[in] A
  57. *> \verbatim
  58. *> A is COMPLEX*16 array, dimension (LDA, N)
  59. *> On entry, the LU part of the factorization of the n-by-n
  60. *> matrix A computed by ZGETC2: A = P * L * U * Q
  61. *> \endverbatim
  62. *>
  63. *> \param[in] LDA
  64. *> \verbatim
  65. *> LDA is INTEGER
  66. *> The leading dimension of the array A. LDA >= max(1, N).
  67. *> \endverbatim
  68. *>
  69. *> \param[in,out] RHS
  70. *> \verbatim
  71. *> RHS is COMPLEX*16 array, dimension N.
  72. *> On entry, the right hand side vector b.
  73. *> On exit, the solution vector X.
  74. *> \endverbatim
  75. *>
  76. *> \param[in] IPIV
  77. *> \verbatim
  78. *> IPIV is INTEGER array, dimension (N).
  79. *> The pivot indices; for 1 <= i <= N, row i of the
  80. *> matrix has been interchanged with row IPIV(i).
  81. *> \endverbatim
  82. *>
  83. *> \param[in] JPIV
  84. *> \verbatim
  85. *> JPIV is INTEGER array, dimension (N).
  86. *> The pivot indices; for 1 <= j <= N, column j of the
  87. *> matrix has been interchanged with column JPIV(j).
  88. *> \endverbatim
  89. *>
  90. *> \param[out] SCALE
  91. *> \verbatim
  92. *> SCALE is DOUBLE PRECISION
  93. *> On exit, SCALE contains the scale factor. SCALE is chosen
  94. *> 0 <= SCALE <= 1 to prevent overflow in the solution.
  95. *> \endverbatim
  96. *
  97. * Authors:
  98. * ========
  99. *
  100. *> \author Univ. of Tennessee
  101. *> \author Univ. of California Berkeley
  102. *> \author Univ. of Colorado Denver
  103. *> \author NAG Ltd.
  104. *
  105. *> \date November 2017
  106. *
  107. *> \ingroup complex16GEauxiliary
  108. *
  109. *> \par Contributors:
  110. * ==================
  111. *>
  112. *> Bo Kagstrom and Peter Poromaa, Department of Computing Science,
  113. *> Umea University, S-901 87 Umea, Sweden.
  114. *
  115. * =====================================================================
  116. SUBROUTINE ZGESC2( N, A, LDA, RHS, IPIV, JPIV, SCALE )
  117. *
  118. * -- LAPACK auxiliary routine (version 3.8.0) --
  119. * -- LAPACK is a software package provided by Univ. of Tennessee, --
  120. * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  121. * November 2017
  122. *
  123. * .. Scalar Arguments ..
  124. INTEGER LDA, N
  125. DOUBLE PRECISION SCALE
  126. * ..
  127. * .. Array Arguments ..
  128. INTEGER IPIV( * ), JPIV( * )
  129. COMPLEX*16 A( LDA, * ), RHS( * )
  130. * ..
  131. *
  132. * =====================================================================
  133. *
  134. * .. Parameters ..
  135. DOUBLE PRECISION ZERO, ONE, TWO
  136. PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0 )
  137. * ..
  138. * .. Local Scalars ..
  139. INTEGER I, J
  140. DOUBLE PRECISION BIGNUM, EPS, SMLNUM
  141. COMPLEX*16 TEMP
  142. * ..
  143. * .. External Subroutines ..
  144. EXTERNAL ZLASWP, ZSCAL, DLABAD
  145. * ..
  146. * .. External Functions ..
  147. INTEGER IZAMAX
  148. DOUBLE PRECISION DLAMCH
  149. EXTERNAL IZAMAX, DLAMCH
  150. * ..
  151. * .. Intrinsic Functions ..
  152. INTRINSIC ABS, DBLE, DCMPLX
  153. * ..
  154. * .. Executable Statements ..
  155. *
  156. * Set constant to control overflow
  157. *
  158. EPS = DLAMCH( 'P' )
  159. SMLNUM = DLAMCH( 'S' ) / EPS
  160. BIGNUM = ONE / SMLNUM
  161. CALL DLABAD( SMLNUM, BIGNUM )
  162. *
  163. * Apply permutations IPIV to RHS
  164. *
  165. CALL ZLASWP( 1, RHS, LDA, 1, N-1, IPIV, 1 )
  166. *
  167. * Solve for L part
  168. *
  169. DO 20 I = 1, N - 1
  170. DO 10 J = I + 1, N
  171. RHS( J ) = RHS( J ) - A( J, I )*RHS( I )
  172. 10 CONTINUE
  173. 20 CONTINUE
  174. *
  175. * Solve for U part
  176. *
  177. SCALE = ONE
  178. *
  179. * Check for scaling
  180. *
  181. I = IZAMAX( N, RHS, 1 )
  182. IF( TWO*SMLNUM*ABS( RHS( I ) ).GT.ABS( A( N, N ) ) ) THEN
  183. TEMP = DCMPLX( ONE / TWO, ZERO ) / ABS( RHS( I ) )
  184. CALL ZSCAL( N, TEMP, RHS( 1 ), 1 )
  185. SCALE = SCALE*DBLE( TEMP )
  186. END IF
  187. DO 40 I = N, 1, -1
  188. TEMP = DCMPLX( ONE, ZERO ) / A( I, I )
  189. RHS( I ) = RHS( I )*TEMP
  190. DO 30 J = I + 1, N
  191. RHS( I ) = RHS( I ) - RHS( J )*( A( I, J )*TEMP )
  192. 30 CONTINUE
  193. 40 CONTINUE
  194. *
  195. * Apply permutations JPIV to the solution (RHS)
  196. *
  197. CALL ZLASWP( 1, RHS, LDA, 1, N-1, JPIV, -1 )
  198. RETURN
  199. *
  200. * End of ZGESC2
  201. *
  202. END