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

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199
  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. *> \ingroup complex16GEauxiliary
  106. *
  107. *> \par Contributors:
  108. * ==================
  109. *>
  110. *> Bo Kagstrom and Peter Poromaa, Department of Computing Science,
  111. *> Umea University, S-901 87 Umea, Sweden.
  112. *
  113. * =====================================================================
  114. SUBROUTINE ZGESC2( N, A, LDA, RHS, IPIV, JPIV, SCALE )
  115. *
  116. * -- LAPACK auxiliary routine --
  117. * -- LAPACK is a software package provided by Univ. of Tennessee, --
  118. * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  119. *
  120. * .. Scalar Arguments ..
  121. INTEGER LDA, N
  122. DOUBLE PRECISION SCALE
  123. * ..
  124. * .. Array Arguments ..
  125. INTEGER IPIV( * ), JPIV( * )
  126. COMPLEX*16 A( LDA, * ), RHS( * )
  127. * ..
  128. *
  129. * =====================================================================
  130. *
  131. * .. Parameters ..
  132. DOUBLE PRECISION ZERO, ONE, TWO
  133. PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0 )
  134. * ..
  135. * .. Local Scalars ..
  136. INTEGER I, J
  137. DOUBLE PRECISION BIGNUM, EPS, SMLNUM
  138. COMPLEX*16 TEMP
  139. * ..
  140. * .. External Subroutines ..
  141. EXTERNAL ZLASWP, ZSCAL, DLABAD
  142. * ..
  143. * .. External Functions ..
  144. INTEGER IZAMAX
  145. DOUBLE PRECISION DLAMCH
  146. EXTERNAL IZAMAX, DLAMCH
  147. * ..
  148. * .. Intrinsic Functions ..
  149. INTRINSIC ABS, DBLE, DCMPLX
  150. * ..
  151. * .. Executable Statements ..
  152. *
  153. * Set constant to control overflow
  154. *
  155. EPS = DLAMCH( 'P' )
  156. SMLNUM = DLAMCH( 'S' ) / EPS
  157. BIGNUM = ONE / SMLNUM
  158. CALL DLABAD( SMLNUM, BIGNUM )
  159. *
  160. * Apply permutations IPIV to RHS
  161. *
  162. CALL ZLASWP( 1, RHS, LDA, 1, N-1, IPIV, 1 )
  163. *
  164. * Solve for L part
  165. *
  166. DO 20 I = 1, N - 1
  167. DO 10 J = I + 1, N
  168. RHS( J ) = RHS( J ) - A( J, I )*RHS( I )
  169. 10 CONTINUE
  170. 20 CONTINUE
  171. *
  172. * Solve for U part
  173. *
  174. SCALE = ONE
  175. *
  176. * Check for scaling
  177. *
  178. I = IZAMAX( N, RHS, 1 )
  179. IF( TWO*SMLNUM*ABS( RHS( I ) ).GT.ABS( A( N, N ) ) ) THEN
  180. TEMP = DCMPLX( ONE / TWO, ZERO ) / ABS( RHS( I ) )
  181. CALL ZSCAL( N, TEMP, RHS( 1 ), 1 )
  182. SCALE = SCALE*DBLE( TEMP )
  183. END IF
  184. DO 40 I = N, 1, -1
  185. TEMP = DCMPLX( ONE, ZERO ) / A( I, I )
  186. RHS( I ) = RHS( I )*TEMP
  187. DO 30 J = I + 1, N
  188. RHS( I ) = RHS( I ) - RHS( J )*( A( I, J )*TEMP )
  189. 30 CONTINUE
  190. 40 CONTINUE
  191. *
  192. * Apply permutations JPIV to the solution (RHS)
  193. *
  194. CALL ZLASWP( 1, RHS, LDA, 1, N-1, JPIV, -1 )
  195. RETURN
  196. *
  197. * End of ZGESC2
  198. *
  199. END