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.

ztrsmf.f 16 kB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457
  1. SUBROUTINE ZTRSMF ( SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA,
  2. $ B, LDB )
  3. * .. Scalar Arguments ..
  4. IMPLICIT NONE
  5. CHARACTER*1 SIDE, UPLO, TRANSA, DIAG
  6. INTEGER M, N, LDA, LDB
  7. COMPLEX*16 ALPHA
  8. * .. Array Arguments ..
  9. COMPLEX*16 A( LDA, * ), B( LDB, * )
  10. * ..
  11. *
  12. * Purpose
  13. * =======
  14. *
  15. * ZTRSM solves one of the matrix equations
  16. *
  17. * op( A )*X = alpha*B, or X*op( A ) = alpha*B,
  18. *
  19. * where alpha is a scalar, X and B are m by n matrices, A is a unit, or
  20. * non-unit, upper or lower triangular matrix and op( A ) is one of
  21. *
  22. * op( A ) = A or op( A ) = A' or op( A ) = conjg( A' ).
  23. *
  24. * The matrix X is overwritten on B.
  25. *
  26. * Parameters
  27. * ==========
  28. *
  29. * SIDE - CHARACTER*1.
  30. * On entry, SIDE specifies whether op( A ) appears on the left
  31. * or right of X as follows:
  32. *
  33. * SIDE = 'L' or 'l' op( A )*X = alpha*B.
  34. *
  35. * SIDE = 'R' or 'r' X*op( A ) = alpha*B.
  36. *
  37. * Unchanged on exit.
  38. *
  39. * UPLO - CHARACTER*1.
  40. * On entry, UPLO specifies whether the matrix A is an upper or
  41. * lower triangular matrix as follows:
  42. *
  43. * UPLO = 'U' or 'u' A is an upper triangular matrix.
  44. *
  45. * UPLO = 'L' or 'l' A is a lower triangular matrix.
  46. *
  47. * Unchanged on exit.
  48. *
  49. * TRANSA - CHARACTER*1.
  50. * On entry, TRANSA specifies the form of op( A ) to be used in
  51. * the matrix multiplication as follows:
  52. *
  53. * TRANSA = 'N' or 'n' op( A ) = A.
  54. *
  55. * TRANSA = 'T' or 't' op( A ) = A'.
  56. *
  57. * TRANSA = 'C' or 'c' op( A ) = conjg( A' ).
  58. *
  59. * Unchanged on exit.
  60. *
  61. * DIAG - CHARACTER*1.
  62. * On entry, DIAG specifies whether or not A is unit triangular
  63. * as follows:
  64. *
  65. * DIAG = 'U' or 'u' A is assumed to be unit triangular.
  66. *
  67. * DIAG = 'N' or 'n' A is not assumed to be unit
  68. * triangular.
  69. *
  70. * Unchanged on exit.
  71. *
  72. * M - INTEGER.
  73. * On entry, M specifies the number of rows of B. M must be at
  74. * least zero.
  75. * Unchanged on exit.
  76. *
  77. * N - INTEGER.
  78. * On entry, N specifies the number of columns of B. N must be
  79. * at least zero.
  80. * Unchanged on exit.
  81. *
  82. * ALPHA - COMPLEX*16 .
  83. * On entry, ALPHA specifies the scalar alpha. When alpha is
  84. * zero then A is not referenced and B need not be set before
  85. * entry.
  86. * Unchanged on exit.
  87. *
  88. * A - COMPLEX*16 array of DIMENSION ( LDA, k ), where k is m
  89. * when SIDE = 'L' or 'l' and is n when SIDE = 'R' or 'r'.
  90. * Before entry with UPLO = 'U' or 'u', the leading k by k
  91. * upper triangular part of the array A must contain the upper
  92. * triangular matrix and the strictly lower triangular part of
  93. * A is not referenced.
  94. * Before entry with UPLO = 'L' or 'l', the leading k by k
  95. * lower triangular part of the array A must contain the lower
  96. * triangular matrix and the strictly upper triangular part of
  97. * A is not referenced.
  98. * Note that when DIAG = 'U' or 'u', the diagonal elements of
  99. * A are not referenced either, but are assumed to be unity.
  100. * Unchanged on exit.
  101. *
  102. * LDA - INTEGER.
  103. * On entry, LDA specifies the first dimension of A as declared
  104. * in the calling (sub) program. When SIDE = 'L' or 'l' then
  105. * LDA must be at least max( 1, m ), when SIDE = 'R' or 'r'
  106. * then LDA must be at least max( 1, n ).
  107. * Unchanged on exit.
  108. *
  109. * B - COMPLEX*16 array of DIMENSION ( LDB, n ).
  110. * Before entry, the leading m by n part of the array B must
  111. * contain the right-hand side matrix B, and on exit is
  112. * overwritten by the solution matrix X.
  113. *
  114. * LDB - INTEGER.
  115. * On entry, LDB specifies the first dimension of B as declared
  116. * in the calling (sub) program. LDB must be at least
  117. * max( 1, m ).
  118. * Unchanged on exit.
  119. *
  120. *
  121. * Level 3 Blas routine.
  122. *
  123. * -- Written on 8-February-1989.
  124. * Jack Dongarra, Argonne National Laboratory.
  125. * Iain Duff, AERE Harwell.
  126. * Jeremy Du Croz, Numerical Algorithms Group Ltd.
  127. * Sven Hammarling, Numerical Algorithms Group Ltd.
  128. *
  129. *
  130. * .. External Functions ..
  131. LOGICAL LSAME
  132. EXTERNAL LSAME
  133. * .. External Subroutines ..
  134. EXTERNAL XERBLA
  135. * .. Intrinsic Functions ..
  136. INTRINSIC DCONJG, MAX
  137. * .. Local Scalars ..
  138. LOGICAL LSIDE, NOCONJ, NOUNIT, UPPER
  139. INTEGER I, INFO, J, K, NROWA
  140. COMPLEX*16 TEMP
  141. * .. Parameters ..
  142. COMPLEX*16 ONE
  143. PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ) )
  144. COMPLEX*16 ZERO
  145. PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) )
  146. * ..
  147. * .. Executable Statements ..
  148. *
  149. * Test the input parameters.
  150. *
  151. LSIDE = LSAME( SIDE , 'L' )
  152. IF( LSIDE )THEN
  153. NROWA = M
  154. ELSE
  155. NROWA = N
  156. END IF
  157. NOCONJ = (LSAME( TRANSA, 'N' ) .OR. LSAME( TRANSA, 'T' ))
  158. NOUNIT = LSAME( DIAG , 'N' )
  159. UPPER = LSAME( UPLO , 'U' )
  160. *
  161. INFO = 0
  162. IF( ( .NOT.LSIDE ).AND.
  163. $ ( .NOT.LSAME( SIDE , 'R' ) ) )THEN
  164. INFO = 1
  165. ELSE IF( ( .NOT.UPPER ).AND.
  166. $ ( .NOT.LSAME( UPLO , 'L' ) ) )THEN
  167. INFO = 2
  168. ELSE IF( ( .NOT.LSAME( TRANSA, 'N' ) ).AND.
  169. $ ( .NOT.LSAME( TRANSA, 'T' ) ).AND.
  170. $ ( .NOT.LSAME( TRANSA, 'R' ) ).AND.
  171. $ ( .NOT.LSAME( TRANSA, 'C' ) ) )THEN
  172. INFO = 3
  173. ELSE IF( ( .NOT.LSAME( DIAG , 'U' ) ).AND.
  174. $ ( .NOT.LSAME( DIAG , 'N' ) ) )THEN
  175. INFO = 4
  176. ELSE IF( M .LT.0 )THEN
  177. INFO = 5
  178. ELSE IF( N .LT.0 )THEN
  179. INFO = 6
  180. ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN
  181. INFO = 9
  182. ELSE IF( LDB.LT.MAX( 1, M ) )THEN
  183. INFO = 11
  184. END IF
  185. IF( INFO.NE.0 )THEN
  186. CALL XERBLA( 'ZTRSM ', INFO )
  187. RETURN
  188. END IF
  189. *
  190. * Quick return if possible.
  191. *
  192. IF( N.EQ.0 )
  193. $ RETURN
  194. *
  195. * And when alpha.eq.zero.
  196. *
  197. IF( ALPHA.EQ.ZERO )THEN
  198. DO 20, J = 1, N
  199. DO 10, I = 1, M
  200. B( I, J ) = ZERO
  201. 10 CONTINUE
  202. 20 CONTINUE
  203. RETURN
  204. END IF
  205. *
  206. * Start the operations.
  207. *
  208. IF( LSIDE )THEN
  209. IF( LSAME( TRANSA, 'N' ) .OR. LSAME( TRANSA, 'R' ) )THEN
  210. *
  211. * Form B := alpha*inv( A )*B.
  212. *
  213. IF( UPPER )THEN
  214. DO 60, J = 1, N
  215. IF( ALPHA.NE.ONE )THEN
  216. DO 30, I = 1, M
  217. B( I, J ) = ALPHA*B( I, J )
  218. 30 CONTINUE
  219. END IF
  220. DO 50, K = M, 1, -1
  221. IF( B( K, J ).NE.ZERO )THEN
  222. IF( NOUNIT ) THEN
  223. IF (NOCONJ) THEN
  224. B( K, J ) = B( K, J )/A( K, K )
  225. ELSE
  226. B( K, J ) = B( K, J )/DCONJG(A( K, K ))
  227. ENDIF
  228. ENDIF
  229. IF (NOCONJ) THEN
  230. DO 40, I = 1, K - 1
  231. B( I, J ) = B( I, J ) - B( K, J )*A( I, K )
  232. 40 CONTINUE
  233. ELSE
  234. DO 45, I = 1, K - 1
  235. B( I, J ) = B( I, J ) - B( K, J )*DCONJG(A( I, K ))
  236. 45 CONTINUE
  237. ENDIF
  238. END IF
  239. 50 CONTINUE
  240. 60 CONTINUE
  241. ELSE
  242. DO 100, J = 1, N
  243. IF( ALPHA.NE.ONE )THEN
  244. DO 70, I = 1, M
  245. B( I, J ) = ALPHA*B( I, J )
  246. 70 CONTINUE
  247. END IF
  248. DO 90 K = 1, M
  249. IF (NOCONJ) THEN
  250. IF( B( K, J ).NE.ZERO )THEN
  251. IF( NOUNIT )
  252. $ B( K, J ) = B( K, J )/A( K, K )
  253. DO 80, I = K + 1, M
  254. B( I, J ) = B( I, J ) - B( K, J )*A( I, K )
  255. 80 CONTINUE
  256. END IF
  257. ELSE
  258. IF( B( K, J ).NE.ZERO )THEN
  259. IF( NOUNIT )
  260. $ B( K, J ) = B( K, J )/DCONJG(A( K, K ))
  261. DO 85, I = K + 1, M
  262. B( I, J ) = B( I, J ) - B( K, J )*DCONJG(A( I, K ))
  263. 85 CONTINUE
  264. END IF
  265. ENDIF
  266. 90 CONTINUE
  267. 100 CONTINUE
  268. END IF
  269. ELSE
  270. *
  271. * Form B := alpha*inv( A' )*B
  272. * or B := alpha*inv( conjg( A' ) )*B.
  273. *
  274. IF( UPPER )THEN
  275. DO 140, J = 1, N
  276. DO 130, I = 1, M
  277. TEMP = ALPHA*B( I, J )
  278. IF( NOCONJ )THEN
  279. DO 110, K = 1, I - 1
  280. TEMP = TEMP - A( K, I )*B( K, J )
  281. 110 CONTINUE
  282. IF( NOUNIT )
  283. $ TEMP = TEMP/A( I, I )
  284. ELSE
  285. DO 120, K = 1, I - 1
  286. TEMP = TEMP - DCONJG( A( K, I ) )*B( K, J )
  287. 120 CONTINUE
  288. IF( NOUNIT )
  289. $ TEMP = TEMP/DCONJG( A( I, I ) )
  290. END IF
  291. B( I, J ) = TEMP
  292. 130 CONTINUE
  293. 140 CONTINUE
  294. ELSE
  295. DO 180, J = 1, N
  296. DO 170, I = M, 1, -1
  297. TEMP = ALPHA*B( I, J )
  298. IF( NOCONJ )THEN
  299. DO 150, K = I + 1, M
  300. TEMP = TEMP - A( K, I )*B( K, J )
  301. 150 CONTINUE
  302. IF( NOUNIT )
  303. $ TEMP = TEMP/A( I, I )
  304. ELSE
  305. DO 160, K = I + 1, M
  306. TEMP = TEMP - DCONJG( A( K, I ) )*B( K, J )
  307. 160 CONTINUE
  308. IF( NOUNIT )
  309. $ TEMP = TEMP/DCONJG( A( I, I ) )
  310. END IF
  311. B( I, J ) = TEMP
  312. 170 CONTINUE
  313. 180 CONTINUE
  314. END IF
  315. END IF
  316. ELSE
  317. IF( LSAME( TRANSA, 'N' ) .OR. LSAME( TRANSA, 'R' ) )THEN
  318. *
  319. * Form B := alpha*B*inv( A ).
  320. *
  321. IF( UPPER )THEN
  322. DO 230, J = 1, N
  323. IF( ALPHA.NE.ONE )THEN
  324. DO 190, I = 1, M
  325. B( I, J ) = ALPHA*B( I, J )
  326. 190 CONTINUE
  327. END IF
  328. DO 210, K = 1, J - 1
  329. IF( A( K, J ).NE.ZERO )THEN
  330. IF (NOCONJ) THEN
  331. DO 200, I = 1, M
  332. B( I, J ) = B( I, J ) - A( K, J )*B( I, K )
  333. 200 CONTINUE
  334. ELSE
  335. DO 205, I = 1, M
  336. B( I, J ) = B( I, J ) - DCONJG(A( K, J ))*B( I, K )
  337. 205 CONTINUE
  338. ENDIF
  339. END IF
  340. 210 CONTINUE
  341. IF( NOUNIT )THEN
  342. IF (NOCONJ) THEN
  343. TEMP = ONE/A( J, J )
  344. ELSE
  345. TEMP = ONE/DCONJG(A( J, J ))
  346. ENDIF
  347. DO 220, I = 1, M
  348. B( I, J ) = TEMP*B( I, J )
  349. 220 CONTINUE
  350. END IF
  351. 230 CONTINUE
  352. ELSE
  353. DO 280, J = N, 1, -1
  354. IF( ALPHA.NE.ONE )THEN
  355. DO 240, I = 1, M
  356. B( I, J ) = ALPHA*B( I, J )
  357. 240 CONTINUE
  358. END IF
  359. DO 260, K = J + 1, N
  360. IF( A( K, J ).NE.ZERO )THEN
  361. IF (NOCONJ) THEN
  362. DO 250, I = 1, M
  363. B( I, J ) = B( I, J ) - A( K, J )*B( I, K )
  364. 250 CONTINUE
  365. ELSE
  366. DO 255, I = 1, M
  367. B( I, J ) = B( I, J ) - DCONJG(A( K, J ))*B( I, K )
  368. 255 CONTINUE
  369. ENDIF
  370. END IF
  371. 260 CONTINUE
  372. IF( NOUNIT )THEN
  373. IF (NOCONJ) THEN
  374. TEMP = ONE/A( J, J )
  375. ELSE
  376. TEMP = ONE/DCONJG(A( J, J ))
  377. ENDIF
  378. DO 270, I = 1, M
  379. B( I, J ) = TEMP*B( I, J )
  380. 270 CONTINUE
  381. END IF
  382. 280 CONTINUE
  383. END IF
  384. ELSE
  385. *
  386. * Form B := alpha*B*inv( A' )
  387. * or B := alpha*B*inv( conjg( A' ) ).
  388. *
  389. IF( UPPER )THEN
  390. DO 330, K = N, 1, -1
  391. IF( NOUNIT )THEN
  392. IF( NOCONJ )THEN
  393. TEMP = ONE/A( K, K )
  394. ELSE
  395. TEMP = ONE/DCONJG( A( K, K ) )
  396. END IF
  397. DO 290, I = 1, M
  398. B( I, K ) = TEMP*B( I, K )
  399. 290 CONTINUE
  400. END IF
  401. DO 310, J = 1, K - 1
  402. IF( A( J, K ).NE.ZERO )THEN
  403. IF( NOCONJ )THEN
  404. TEMP = A( J, K )
  405. ELSE
  406. TEMP = DCONJG( A( J, K ) )
  407. END IF
  408. DO 300, I = 1, M
  409. B( I, J ) = B( I, J ) - TEMP*B( I, K )
  410. 300 CONTINUE
  411. END IF
  412. 310 CONTINUE
  413. IF( ALPHA.NE.ONE )THEN
  414. DO 320, I = 1, M
  415. B( I, K ) = ALPHA*B( I, K )
  416. 320 CONTINUE
  417. END IF
  418. 330 CONTINUE
  419. ELSE
  420. DO 380, K = 1, N
  421. IF( NOUNIT )THEN
  422. IF( NOCONJ )THEN
  423. TEMP = ONE/A( K, K )
  424. ELSE
  425. TEMP = ONE/DCONJG( A( K, K ) )
  426. END IF
  427. DO 340, I = 1, M
  428. B( I, K ) = TEMP*B( I, K )
  429. 340 CONTINUE
  430. END IF
  431. DO 360, J = K + 1, N
  432. IF( A( J, K ).NE.ZERO )THEN
  433. IF( NOCONJ )THEN
  434. TEMP = A( J, K )
  435. ELSE
  436. TEMP = DCONJG( A( J, K ) )
  437. END IF
  438. DO 350, I = 1, M
  439. B( I, J ) = B( I, J ) - TEMP*B( I, K )
  440. 350 CONTINUE
  441. END IF
  442. 360 CONTINUE
  443. IF( ALPHA.NE.ONE )THEN
  444. DO 370, I = 1, M
  445. B( I, K ) = ALPHA*B( I, K )
  446. 370 CONTINUE
  447. END IF
  448. 380 CONTINUE
  449. END IF
  450. END IF
  451. END IF
  452. *
  453. RETURN
  454. *
  455. * End of ZTRSM .
  456. *
  457. END