Télécharger dtrsv.eso

Retour à la liste

Numérotation des lignes :

  1. C DTRSV SOURCE CHAT 06/03/29 21:18:59 5360
  2. SUBROUTINE DTRSV ( UPLO, TRANS, DIAG, N, A, LDA, X, INCX )
  3. IMPLICIT INTEGER(I-N)
  4. IMPLICIT REAL*8 (A-H,O-Z)
  5. -INC CCOPTIO
  6. * .. Scalar Arguments ..
  7. INTEGER INCX, LDA, N
  8. CHARACTER*1 DIAG, TRANS, UPLO
  9. * .. Array Arguments ..
  10. REAL*8 A( LDA, * ), X( * )
  11. * ..
  12. *
  13. * Purpose
  14. * =======
  15. *
  16. * DTRSV solves one of the systems of equations
  17. *
  18. * A*x = b, or A'*x = b,
  19. *
  20. * where b and x are n element vectors and A is an n by n unit, or
  21. * non-unit, upper or lower triangular matrix.
  22. *
  23. * No test for singularity or near-singularity is included in this
  24. * routine. Such tests must be performed before calling this routine.
  25. *
  26. * Parameters
  27. * ==========
  28. *
  29. * UPLO - CHARACTER*1.
  30. * On entry, UPLO specifies whether the matrix is an upper or
  31. * lower triangular matrix as follows:
  32. *
  33. * UPLO = 'U' or 'u' A is an upper triangular matrix.
  34. *
  35. * UPLO = 'L' or 'l' A is a lower triangular matrix.
  36. *
  37. * Unchanged on exit.
  38. *
  39. * TRANS - CHARACTER*1.
  40. * On entry, TRANS specifies the equations to be solved as
  41. * follows:
  42. *
  43. * TRANS = 'N' or 'n' A*x = b.
  44. *
  45. * TRANS = 'T' or 't' A'*x = b.
  46. *
  47. * TRANS = 'C' or 'c' A'*x = b.
  48. *
  49. * Unchanged on exit.
  50. *
  51. * DIAG - CHARACTER*1.
  52. * On entry, DIAG specifies whether or not A is unit
  53. * triangular as follows:
  54. *
  55. * DIAG = 'U' or 'u' A is assumed to be unit triangular.
  56. *
  57. * DIAG = 'N' or 'n' A is not assumed to be unit
  58. * triangular.
  59. *
  60. * Unchanged on exit.
  61. *
  62. * N - INTEGER.
  63. * On entry, N specifies the order of the matrix A.
  64. * N must be at least zero.
  65. * Unchanged on exit.
  66. *
  67. * A - DOUBLE PRECISION array of DIMENSION ( LDA, n ).
  68. * Before entry with UPLO = 'U' or 'u', the leading n by n
  69. * upper triangular part of the array A must contain the upper
  70. * triangular matrix and the strictly lower triangular part of
  71. * A is not referenced.
  72. * Before entry with UPLO = 'L' or 'l', the leading n by n
  73. * lower triangular part of the array A must contain the lower
  74. * triangular matrix and the strictly upper triangular part of
  75. * A is not referenced.
  76. * Note that when DIAG = 'U' or 'u', the diagonal elements of
  77. * A are not referenced either, but are assumed to be unity.
  78. * Unchanged on exit.
  79. *
  80. * LDA - INTEGER.
  81. * On entry, LDA specifies the first dimension of A as declared
  82. * in the calling (sub) program. LDA must be at least
  83. * max( 1, n ).
  84. * Unchanged on exit.
  85. *
  86. * X - DOUBLE PRECISION array of dimension at least
  87. * ( 1 + ( n - 1 )*abs( INCX ) ).
  88. * Before entry, the incremented array X must contain the n
  89. * element right-hand side vector b. On exit, X is overwritten
  90. * with the solution vector x.
  91. *
  92. * INCX - INTEGER.
  93. * On entry, INCX specifies the increment for the elements of
  94. * X. INCX must not be zero.
  95. * Unchanged on exit.
  96. *
  97. *
  98. * Level 2 Blas routine.
  99. *
  100. * -- Written on 22-October-1986.
  101. * Jack Dongarra, Argonne National Lab.
  102. * Jeremy Du Croz, Nag Central Office.
  103. * Sven Hammarling, Nag Central Office.
  104. * Richard Hanson, Sandia National Labs.
  105. C modified 5/11/98 : double precision -> real*8
  106. * commented INTRINSIC (not esope compatible)
  107. C added error handling
  108. *
  109. *
  110. * .. Parameters ..
  111. REAL*8 ZERO
  112. PARAMETER ( ZERO = 0.0D+0 )
  113. * .. Local Scalars ..
  114. REAL*8 TEMP
  115. INTEGER I, INFO, IX, J, JX, KX
  116. LOGICAL NOUNIT
  117. * .. External Functions ..
  118. LOGICAL LSAME
  119. EXTERNAL LSAME
  120. * .. Intrinsic Functions ..
  121. * INTRINSIC MAX
  122. * ..
  123. * .. Executable Statements ..
  124. *
  125. * Test the input parameters.
  126. *
  127. INFO = 0
  128. IF ( .NOT.LSAME( UPLO , 'U' ).AND.
  129. $ .NOT.LSAME( UPLO , 'L' ) )THEN
  130. INFO = 1
  131. ELSE IF( .NOT.LSAME( TRANS, 'N' ).AND.
  132. $ .NOT.LSAME( TRANS, 'T' ).AND.
  133. $ .NOT.LSAME( TRANS, 'C' ) )THEN
  134. INFO = 2
  135. ELSE IF( .NOT.LSAME( DIAG , 'U' ).AND.
  136. $ .NOT.LSAME( DIAG , 'N' ) )THEN
  137. INFO = 3
  138. ELSE IF( N.LT.0 )THEN
  139. INFO = 4
  140. ELSE IF( LDA.LT.MAX( 1, N ) )THEN
  141. INFO = 6
  142. ELSE IF( INCX.EQ.0 )THEN
  143. INFO = 8
  144. END IF
  145. IF( INFO.NE.0 )THEN
  146. GOTO 9999
  147. END IF
  148. *
  149. * Quick return if possible.
  150. *
  151. IF( N.EQ.0 )
  152. $ RETURN
  153. *
  154. NOUNIT = LSAME( DIAG, 'N' )
  155. *
  156. * Set up the start point in X if the increment is not unity. This
  157. * will be ( N - 1 )*INCX too small for descending loops.
  158. *
  159. IF( INCX.LE.0 )THEN
  160. KX = 1 - ( N - 1 )*INCX
  161. ELSE IF( INCX.NE.1 )THEN
  162. KX = 1
  163. END IF
  164. *
  165. * Start the operations. In this version the elements of A are
  166. * accessed sequentially with one pass through A.
  167. *
  168. IF( LSAME( TRANS, 'N' ) )THEN
  169. *
  170. * Form x := inv( A )*x.
  171. *
  172. IF( LSAME( UPLO, 'U' ) )THEN
  173. IF( INCX.EQ.1 )THEN
  174. DO 20, J = N, 1, -1
  175. IF( X( J ).NE.ZERO )THEN
  176. IF( NOUNIT )
  177. $ X( J ) = X( J )/A( J, J )
  178. TEMP = X( J )
  179. DO 10, I = J - 1, 1, -1
  180. X( I ) = X( I ) - TEMP*A( I, J )
  181. 10 CONTINUE
  182. END IF
  183. 20 CONTINUE
  184. ELSE
  185. JX = KX + ( N - 1 )*INCX
  186. DO 40, J = N, 1, -1
  187. IF( X( JX ).NE.ZERO )THEN
  188. IF( NOUNIT )
  189. $ X( JX ) = X( JX )/A( J, J )
  190. TEMP = X( JX )
  191. IX = JX
  192. DO 30, I = J - 1, 1, -1
  193. IX = IX - INCX
  194. X( IX ) = X( IX ) - TEMP*A( I, J )
  195. 30 CONTINUE
  196. END IF
  197. JX = JX - INCX
  198. 40 CONTINUE
  199. END IF
  200. ELSE
  201. IF( INCX.EQ.1 )THEN
  202. DO 60, J = 1, N
  203. IF( X( J ).NE.ZERO )THEN
  204. IF( NOUNIT )
  205. $ X( J ) = X( J )/A( J, J )
  206. TEMP = X( J )
  207. DO 50, I = J + 1, N
  208. X( I ) = X( I ) - TEMP*A( I, J )
  209. 50 CONTINUE
  210. END IF
  211. 60 CONTINUE
  212. ELSE
  213. JX = KX
  214. DO 80, J = 1, N
  215. IF( X( JX ).NE.ZERO )THEN
  216. IF( NOUNIT )
  217. $ X( JX ) = X( JX )/A( J, J )
  218. TEMP = X( JX )
  219. IX = JX
  220. DO 70, I = J + 1, N
  221. IX = IX + INCX
  222. X( IX ) = X( IX ) - TEMP*A( I, J )
  223. 70 CONTINUE
  224. END IF
  225. JX = JX + INCX
  226. 80 CONTINUE
  227. END IF
  228. END IF
  229. ELSE
  230. *
  231. * Form x := inv( A' )*x.
  232. *
  233. IF( LSAME( UPLO, 'U' ) )THEN
  234. IF( INCX.EQ.1 )THEN
  235. DO 100, J = 1, N
  236. TEMP = X( J )
  237. DO 90, I = 1, J - 1
  238. TEMP = TEMP - A( I, J )*X( I )
  239. 90 CONTINUE
  240. IF( NOUNIT )
  241. $ TEMP = TEMP/A( J, J )
  242. X( J ) = TEMP
  243. 100 CONTINUE
  244. ELSE
  245. JX = KX
  246. DO 120, J = 1, N
  247. TEMP = X( JX )
  248. IX = KX
  249. DO 110, I = 1, J - 1
  250. TEMP = TEMP - A( I, J )*X( IX )
  251. IX = IX + INCX
  252. 110 CONTINUE
  253. IF( NOUNIT )
  254. $ TEMP = TEMP/A( J, J )
  255. X( JX ) = TEMP
  256. JX = JX + INCX
  257. 120 CONTINUE
  258. END IF
  259. ELSE
  260. IF( INCX.EQ.1 )THEN
  261. DO 140, J = N, 1, -1
  262. TEMP = X( J )
  263. DO 130, I = N, J + 1, -1
  264. TEMP = TEMP - A( I, J )*X( I )
  265. 130 CONTINUE
  266. IF( NOUNIT )
  267. $ TEMP = TEMP/A( J, J )
  268. X( J ) = TEMP
  269. 140 CONTINUE
  270. ELSE
  271. KX = KX + ( N - 1 )*INCX
  272. JX = KX
  273. DO 160, J = N, 1, -1
  274. TEMP = X( JX )
  275. IX = KX
  276. DO 150, I = N, J + 1, -1
  277. TEMP = TEMP - A( I, J )*X( IX )
  278. IX = IX - INCX
  279. 150 CONTINUE
  280. IF( NOUNIT )
  281. $ TEMP = TEMP/A( J, J )
  282. X( JX ) = TEMP
  283. JX = JX - INCX
  284. 160 CONTINUE
  285. END IF
  286. END IF
  287. END IF
  288. *
  289. C
  290. C Normal termination
  291. C
  292. RETURN
  293. C
  294. C Error handling
  295. C
  296. 9999 CONTINUE
  297. WRITE(IOIMP,*) 'Error number :',INFO
  298. WRITE(IOIMP,*) 'in subroutine dtrsv.'
  299. call erreur(21)
  300. RETURN
  301. *
  302. * End of DTRSV .
  303. *
  304. END
  305.  
  306.  
  307.  
  308.  
  309.  
  310.  
  311.  
  312.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales