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.

spst01.f 7.9 kB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303
  1. *> \brief \b SPST01
  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 SPST01( UPLO, N, A, LDA, AFAC, LDAFAC, PERM, LDPERM,
  12. * PIV, RWORK, RESID, RANK )
  13. *
  14. * .. Scalar Arguments ..
  15. * REAL RESID
  16. * INTEGER LDA, LDAFAC, LDPERM, N, RANK
  17. * CHARACTER UPLO
  18. * ..
  19. * .. Array Arguments ..
  20. * REAL A( LDA, * ), AFAC( LDAFAC, * ),
  21. * $ PERM( LDPERM, * ), RWORK( * )
  22. * INTEGER PIV( * )
  23. * ..
  24. *
  25. *
  26. *> \par Purpose:
  27. * =============
  28. *>
  29. *> \verbatim
  30. *>
  31. *> SPST01 reconstructs a symmetric positive semidefinite matrix A
  32. *> from its L or U factors and the permutation matrix P and computes
  33. *> the residual
  34. *> norm( P*L*L'*P' - A ) / ( N * norm(A) * EPS ) or
  35. *> norm( P*U'*U*P' - A ) / ( N * norm(A) * EPS ),
  36. *> where EPS is the machine epsilon.
  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. *> symmetric 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] A
  58. *> \verbatim
  59. *> A is REAL array, dimension (LDA,N)
  60. *> The original symmetric matrix A.
  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] AFAC
  70. *> \verbatim
  71. *> AFAC is REAL array, dimension (LDAFAC,N)
  72. *> The factor L or U from the L*L' or U'*U
  73. *> factorization of A.
  74. *> \endverbatim
  75. *>
  76. *> \param[in] LDAFAC
  77. *> \verbatim
  78. *> LDAFAC is INTEGER
  79. *> The leading dimension of the array AFAC. LDAFAC >= max(1,N).
  80. *> \endverbatim
  81. *>
  82. *> \param[out] PERM
  83. *> \verbatim
  84. *> PERM is REAL array, dimension (LDPERM,N)
  85. *> Overwritten with the reconstructed matrix, and then with the
  86. *> difference P*L*L'*P' - A (or P*U'*U*P' - A)
  87. *> \endverbatim
  88. *>
  89. *> \param[in] LDPERM
  90. *> \verbatim
  91. *> LDPERM is INTEGER
  92. *> The leading dimension of the array PERM.
  93. *> LDAPERM >= max(1,N).
  94. *> \endverbatim
  95. *>
  96. *> \param[in] PIV
  97. *> \verbatim
  98. *> PIV is INTEGER array, dimension (N)
  99. *> PIV is such that the nonzero entries are
  100. *> P( PIV( K ), K ) = 1.
  101. *> \endverbatim
  102. *>
  103. *> \param[out] RWORK
  104. *> \verbatim
  105. *> RWORK is REAL array, dimension (N)
  106. *> \endverbatim
  107. *>
  108. *> \param[out] RESID
  109. *> \verbatim
  110. *> RESID is REAL
  111. *> If UPLO = 'L', norm(L*L' - A) / ( N * norm(A) * EPS )
  112. *> If UPLO = 'U', norm(U'*U - A) / ( N * norm(A) * EPS )
  113. *> \endverbatim
  114. *>
  115. *> \param[in] RANK
  116. *> \verbatim
  117. *> RANK is INTEGER
  118. *> number of nonzero singular values of A.
  119. *> \endverbatim
  120. *
  121. * Authors:
  122. * ========
  123. *
  124. *> \author Univ. of Tennessee
  125. *> \author Univ. of California Berkeley
  126. *> \author Univ. of Colorado Denver
  127. *> \author NAG Ltd.
  128. *
  129. *> \ingroup single_lin
  130. *
  131. * =====================================================================
  132. SUBROUTINE SPST01( UPLO, N, A, LDA, AFAC, LDAFAC, PERM, LDPERM,
  133. $ PIV, RWORK, RESID, RANK )
  134. *
  135. * -- LAPACK test routine --
  136. * -- LAPACK is a software package provided by Univ. of Tennessee, --
  137. * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  138. *
  139. * .. Scalar Arguments ..
  140. REAL RESID
  141. INTEGER LDA, LDAFAC, LDPERM, N, RANK
  142. CHARACTER UPLO
  143. * ..
  144. * .. Array Arguments ..
  145. REAL A( LDA, * ), AFAC( LDAFAC, * ),
  146. $ PERM( LDPERM, * ), RWORK( * )
  147. INTEGER PIV( * )
  148. * ..
  149. *
  150. * =====================================================================
  151. *
  152. * .. Parameters ..
  153. REAL ZERO, ONE
  154. PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 )
  155. * ..
  156. * .. Local Scalars ..
  157. REAL ANORM, EPS, T
  158. INTEGER I, J, K
  159. * ..
  160. * .. External Functions ..
  161. REAL SDOT, SLAMCH, SLANSY
  162. LOGICAL LSAME
  163. EXTERNAL SDOT, SLAMCH, SLANSY, LSAME
  164. * ..
  165. * .. External Subroutines ..
  166. EXTERNAL SSCAL, SSYR, STRMV
  167. * ..
  168. * .. Intrinsic Functions ..
  169. INTRINSIC REAL
  170. * ..
  171. * .. Executable Statements ..
  172. *
  173. * Quick exit if N = 0.
  174. *
  175. IF( N.LE.0 ) THEN
  176. RESID = ZERO
  177. RETURN
  178. END IF
  179. *
  180. * Exit with RESID = 1/EPS if ANORM = 0.
  181. *
  182. EPS = SLAMCH( 'Epsilon' )
  183. ANORM = SLANSY( '1', UPLO, N, A, LDA, RWORK )
  184. IF( ANORM.LE.ZERO ) THEN
  185. RESID = ONE / EPS
  186. RETURN
  187. END IF
  188. *
  189. * Compute the product U'*U, overwriting U.
  190. *
  191. IF( LSAME( UPLO, 'U' ) ) THEN
  192. *
  193. IF( RANK.LT.N ) THEN
  194. DO 110 J = RANK + 1, N
  195. DO 100 I = RANK + 1, J
  196. AFAC( I, J ) = ZERO
  197. 100 CONTINUE
  198. 110 CONTINUE
  199. END IF
  200. *
  201. DO 120 K = N, 1, -1
  202. *
  203. * Compute the (K,K) element of the result.
  204. *
  205. T = SDOT( K, AFAC( 1, K ), 1, AFAC( 1, K ), 1 )
  206. AFAC( K, K ) = T
  207. *
  208. * Compute the rest of column K.
  209. *
  210. CALL STRMV( 'Upper', 'Transpose', 'Non-unit', K-1, AFAC,
  211. $ LDAFAC, AFAC( 1, K ), 1 )
  212. *
  213. 120 CONTINUE
  214. *
  215. * Compute the product L*L', overwriting L.
  216. *
  217. ELSE
  218. *
  219. IF( RANK.LT.N ) THEN
  220. DO 140 J = RANK + 1, N
  221. DO 130 I = J, N
  222. AFAC( I, J ) = ZERO
  223. 130 CONTINUE
  224. 140 CONTINUE
  225. END IF
  226. *
  227. DO 150 K = N, 1, -1
  228. * Add a multiple of column K of the factor L to each of
  229. * columns K+1 through N.
  230. *
  231. IF( K+1.LE.N )
  232. $ CALL SSYR( 'Lower', N-K, ONE, AFAC( K+1, K ), 1,
  233. $ AFAC( K+1, K+1 ), LDAFAC )
  234. *
  235. * Scale column K by the diagonal element.
  236. *
  237. T = AFAC( K, K )
  238. CALL SSCAL( N-K+1, T, AFAC( K, K ), 1 )
  239. 150 CONTINUE
  240. *
  241. END IF
  242. *
  243. * Form P*L*L'*P' or P*U'*U*P'
  244. *
  245. IF( LSAME( UPLO, 'U' ) ) THEN
  246. *
  247. DO 170 J = 1, N
  248. DO 160 I = 1, N
  249. IF( PIV( I ).LE.PIV( J ) ) THEN
  250. IF( I.LE.J ) THEN
  251. PERM( PIV( I ), PIV( J ) ) = AFAC( I, J )
  252. ELSE
  253. PERM( PIV( I ), PIV( J ) ) = AFAC( J, I )
  254. END IF
  255. END IF
  256. 160 CONTINUE
  257. 170 CONTINUE
  258. *
  259. *
  260. ELSE
  261. *
  262. DO 190 J = 1, N
  263. DO 180 I = 1, N
  264. IF( PIV( I ).GE.PIV( J ) ) THEN
  265. IF( I.GE.J ) THEN
  266. PERM( PIV( I ), PIV( J ) ) = AFAC( I, J )
  267. ELSE
  268. PERM( PIV( I ), PIV( J ) ) = AFAC( J, I )
  269. END IF
  270. END IF
  271. 180 CONTINUE
  272. 190 CONTINUE
  273. *
  274. END IF
  275. *
  276. * Compute the difference P*L*L'*P' - A (or P*U'*U*P' - A).
  277. *
  278. IF( LSAME( UPLO, 'U' ) ) THEN
  279. DO 210 J = 1, N
  280. DO 200 I = 1, J
  281. PERM( I, J ) = PERM( I, J ) - A( I, J )
  282. 200 CONTINUE
  283. 210 CONTINUE
  284. ELSE
  285. DO 230 J = 1, N
  286. DO 220 I = J, N
  287. PERM( I, J ) = PERM( I, J ) - A( I, J )
  288. 220 CONTINUE
  289. 230 CONTINUE
  290. END IF
  291. *
  292. * Compute norm( P*L*L'P - A ) / ( N * norm(A) * EPS ), or
  293. * ( P*U'*U*P' - A )/ ( N * norm(A) * EPS ).
  294. *
  295. RESID = SLANSY( '1', UPLO, N, PERM, LDAFAC, RWORK )
  296. *
  297. RESID = ( ( RESID / REAL( N ) ) / ANORM ) / EPS
  298. *
  299. RETURN
  300. *
  301. * End of SPST01
  302. *
  303. END