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

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