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.

sgesc2.f 5.4 kB

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