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.

cppt02.f 5.4 kB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206
  1. *> \brief \b CPPT02
  2. *
  3. * =========== DOCUMENTATION ===========
  4. *
  5. * Online html documentation available at
  6. * http://www.netlib.org/lapack/explore-html/
  7. *
  8. * Definition:
  9. * ===========
  10. *
  11. * SUBROUTINE CPPT02( UPLO, N, NRHS, A, X, LDX, B, LDB, RWORK,
  12. * RESID )
  13. *
  14. * .. Scalar Arguments ..
  15. * CHARACTER UPLO
  16. * INTEGER LDB, LDX, N, NRHS
  17. * REAL RESID
  18. * ..
  19. * .. Array Arguments ..
  20. * REAL RWORK( * )
  21. * COMPLEX A( * ), B( LDB, * ), X( LDX, * )
  22. * ..
  23. *
  24. *
  25. *> \par Purpose:
  26. * =============
  27. *>
  28. *> \verbatim
  29. *>
  30. *> CPPT02 computes the residual in the solution of a Hermitian system
  31. *> of linear equations A*x = b when packed storage is used for the
  32. *> coefficient matrix. The ratio computed is
  33. *>
  34. *> RESID = norm(B - A*X) / ( norm(A) * norm(X) * EPS),
  35. *>
  36. *> where EPS is the machine precision.
  37. *> \endverbatim
  38. *
  39. * Arguments:
  40. * ==========
  41. *
  42. *> \param[in] UPLO
  43. *> \verbatim
  44. *> UPLO is CHARACTER*1
  45. *> Specifies whether the upper or lower triangular part of the
  46. *> Hermitian matrix A is stored:
  47. *> = 'U': Upper triangular
  48. *> = 'L': Lower triangular
  49. *> \endverbatim
  50. *>
  51. *> \param[in] N
  52. *> \verbatim
  53. *> N is INTEGER
  54. *> The number of rows and columns of the matrix A. N >= 0.
  55. *> \endverbatim
  56. *>
  57. *> \param[in] NRHS
  58. *> \verbatim
  59. *> NRHS is INTEGER
  60. *> The number of columns of B, the matrix of right hand sides.
  61. *> NRHS >= 0.
  62. *> \endverbatim
  63. *>
  64. *> \param[in] A
  65. *> \verbatim
  66. *> A is COMPLEX array, dimension (N*(N+1)/2)
  67. *> The original Hermitian matrix A, stored as a packed
  68. *> triangular matrix.
  69. *> \endverbatim
  70. *>
  71. *> \param[in] X
  72. *> \verbatim
  73. *> X is COMPLEX array, dimension (LDX,NRHS)
  74. *> The computed solution vectors for the system of linear
  75. *> equations.
  76. *> \endverbatim
  77. *>
  78. *> \param[in] LDX
  79. *> \verbatim
  80. *> LDX is INTEGER
  81. *> The leading dimension of the array X. LDX >= max(1,N).
  82. *> \endverbatim
  83. *>
  84. *> \param[in,out] B
  85. *> \verbatim
  86. *> B is COMPLEX array, dimension (LDB,NRHS)
  87. *> On entry, the right hand side vectors for the system of
  88. *> linear equations.
  89. *> On exit, B is overwritten with the difference B - A*X.
  90. *> \endverbatim
  91. *>
  92. *> \param[in] LDB
  93. *> \verbatim
  94. *> LDB is INTEGER
  95. *> The leading dimension of the array B. LDB >= max(1,N).
  96. *> \endverbatim
  97. *>
  98. *> \param[out] RWORK
  99. *> \verbatim
  100. *> RWORK is REAL array, dimension (N)
  101. *> \endverbatim
  102. *>
  103. *> \param[out] RESID
  104. *> \verbatim
  105. *> RESID is REAL
  106. *> The maximum over the number of right hand sides of
  107. *> norm(B - A*X) / ( norm(A) * norm(X) * EPS ).
  108. *> \endverbatim
  109. *
  110. * Authors:
  111. * ========
  112. *
  113. *> \author Univ. of Tennessee
  114. *> \author Univ. of California Berkeley
  115. *> \author Univ. of Colorado Denver
  116. *> \author NAG Ltd.
  117. *
  118. *> \date November 2011
  119. *
  120. *> \ingroup complex_lin
  121. *
  122. * =====================================================================
  123. SUBROUTINE CPPT02( UPLO, N, NRHS, A, X, LDX, B, LDB, RWORK,
  124. $ RESID )
  125. *
  126. * -- LAPACK test routine (version 3.4.0) --
  127. * -- LAPACK is a software package provided by Univ. of Tennessee, --
  128. * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  129. * November 2011
  130. *
  131. * .. Scalar Arguments ..
  132. CHARACTER UPLO
  133. INTEGER LDB, LDX, N, NRHS
  134. REAL RESID
  135. * ..
  136. * .. Array Arguments ..
  137. REAL RWORK( * )
  138. COMPLEX A( * ), B( LDB, * ), X( LDX, * )
  139. * ..
  140. *
  141. * =====================================================================
  142. *
  143. * .. Parameters ..
  144. REAL ZERO, ONE
  145. PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 )
  146. COMPLEX CONE
  147. PARAMETER ( CONE = ( 1.0E+0, 0.0E+0 ) )
  148. * ..
  149. * .. Local Scalars ..
  150. INTEGER J
  151. REAL ANORM, BNORM, EPS, XNORM
  152. * ..
  153. * .. External Functions ..
  154. REAL CLANHP, SCASUM, SLAMCH
  155. EXTERNAL CLANHP, SCASUM, SLAMCH
  156. * ..
  157. * .. External Subroutines ..
  158. EXTERNAL CHPMV
  159. * ..
  160. * .. Intrinsic Functions ..
  161. INTRINSIC MAX
  162. * ..
  163. * .. Executable Statements ..
  164. *
  165. * Quick exit if N = 0 or NRHS = 0.
  166. *
  167. IF( N.LE.0 .OR. NRHS.LE.0 ) THEN
  168. RESID = ZERO
  169. RETURN
  170. END IF
  171. *
  172. * Exit with RESID = 1/EPS if ANORM = 0.
  173. *
  174. EPS = SLAMCH( 'Epsilon' )
  175. ANORM = CLANHP( '1', UPLO, N, A, RWORK )
  176. IF( ANORM.LE.ZERO ) THEN
  177. RESID = ONE / EPS
  178. RETURN
  179. END IF
  180. *
  181. * Compute B - A*X for the matrix of right hand sides B.
  182. *
  183. DO 10 J = 1, NRHS
  184. CALL CHPMV( UPLO, N, -CONE, A, X( 1, J ), 1, CONE, B( 1, J ),
  185. $ 1 )
  186. 10 CONTINUE
  187. *
  188. * Compute the maximum over the number of right hand sides of
  189. * norm( B - A*X ) / ( norm(A) * norm(X) * EPS ) .
  190. *
  191. RESID = ZERO
  192. DO 20 J = 1, NRHS
  193. BNORM = SCASUM( N, B( 1, J ), 1 )
  194. XNORM = SCASUM( N, X( 1, J ), 1 )
  195. IF( XNORM.LE.ZERO ) THEN
  196. RESID = ONE / EPS
  197. ELSE
  198. RESID = MAX( RESID, ( ( BNORM/ANORM )/XNORM )/EPS )
  199. END IF
  200. 20 CONTINUE
  201. *
  202. RETURN
  203. *
  204. * End of CPPT02
  205. *
  206. END