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.

dtprfb.f 26 kB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808
  1. *> \brief \b DTPRFB applies a real "triangular-pentagonal" block reflector to a real matrix, which is composed of two blocks.
  2. *
  3. * =========== DOCUMENTATION ===========
  4. *
  5. * Online html documentation available at
  6. * http://www.netlib.org/lapack/explore-html/
  7. *
  8. *> \htmlonly
  9. *> Download DTPRFB + dependencies
  10. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtprfb.f">
  11. *> [TGZ]</a>
  12. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtprfb.f">
  13. *> [ZIP]</a>
  14. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtprfb.f">
  15. *> [TXT]</a>
  16. *> \endhtmlonly
  17. *
  18. * Definition:
  19. * ===========
  20. *
  21. * SUBROUTINE DTPRFB( SIDE, TRANS, DIRECT, STOREV, M, N, K, L,
  22. * V, LDV, T, LDT, A, LDA, B, LDB, WORK, LDWORK )
  23. *
  24. * .. Scalar Arguments ..
  25. * CHARACTER DIRECT, SIDE, STOREV, TRANS
  26. * INTEGER K, L, LDA, LDB, LDT, LDV, LDWORK, M, N
  27. * ..
  28. * .. Array Arguments ..
  29. * DOUBLE PRECISION A( LDA, * ), B( LDB, * ), T( LDT, * ),
  30. * $ V( LDV, * ), WORK( LDWORK, * )
  31. * ..
  32. *
  33. *
  34. *> \par Purpose:
  35. * =============
  36. *>
  37. *> \verbatim
  38. *>
  39. *> DTPRFB applies a real "triangular-pentagonal" block reflector H or its
  40. *> transpose H**T to a real matrix C, which is composed of two
  41. *> blocks A and B, either from the left or right.
  42. *>
  43. *> \endverbatim
  44. *
  45. * Arguments:
  46. * ==========
  47. *
  48. *> \param[in] SIDE
  49. *> \verbatim
  50. *> SIDE is CHARACTER*1
  51. *> = 'L': apply H or H**T from the Left
  52. *> = 'R': apply H or H**T from the Right
  53. *> \endverbatim
  54. *>
  55. *> \param[in] TRANS
  56. *> \verbatim
  57. *> TRANS is CHARACTER*1
  58. *> = 'N': apply H (No transpose)
  59. *> = 'T': apply H**T (Transpose)
  60. *> \endverbatim
  61. *>
  62. *> \param[in] DIRECT
  63. *> \verbatim
  64. *> DIRECT is CHARACTER*1
  65. *> Indicates how H is formed from a product of elementary
  66. *> reflectors
  67. *> = 'F': H = H(1) H(2) . . . H(k) (Forward)
  68. *> = 'B': H = H(k) . . . H(2) H(1) (Backward)
  69. *> \endverbatim
  70. *>
  71. *> \param[in] STOREV
  72. *> \verbatim
  73. *> STOREV is CHARACTER*1
  74. *> Indicates how the vectors which define the elementary
  75. *> reflectors are stored:
  76. *> = 'C': Columns
  77. *> = 'R': Rows
  78. *> \endverbatim
  79. *>
  80. *> \param[in] M
  81. *> \verbatim
  82. *> M is INTEGER
  83. *> The number of rows of the matrix B.
  84. *> M >= 0.
  85. *> \endverbatim
  86. *>
  87. *> \param[in] N
  88. *> \verbatim
  89. *> N is INTEGER
  90. *> The number of columns of the matrix B.
  91. *> N >= 0.
  92. *> \endverbatim
  93. *>
  94. *> \param[in] K
  95. *> \verbatim
  96. *> K is INTEGER
  97. *> The order of the matrix T, i.e. the number of elementary
  98. *> reflectors whose product defines the block reflector.
  99. *> K >= 0.
  100. *> \endverbatim
  101. *>
  102. *> \param[in] L
  103. *> \verbatim
  104. *> L is INTEGER
  105. *> The order of the trapezoidal part of V.
  106. *> K >= L >= 0. See Further Details.
  107. *> \endverbatim
  108. *>
  109. *> \param[in] V
  110. *> \verbatim
  111. *> V is DOUBLE PRECISION array, dimension
  112. *> (LDV,K) if STOREV = 'C'
  113. *> (LDV,M) if STOREV = 'R' and SIDE = 'L'
  114. *> (LDV,N) if STOREV = 'R' and SIDE = 'R'
  115. *> The pentagonal matrix V, which contains the elementary reflectors
  116. *> H(1), H(2), ..., H(K). See Further Details.
  117. *> \endverbatim
  118. *>
  119. *> \param[in] LDV
  120. *> \verbatim
  121. *> LDV is INTEGER
  122. *> The leading dimension of the array V.
  123. *> If STOREV = 'C' and SIDE = 'L', LDV >= max(1,M);
  124. *> if STOREV = 'C' and SIDE = 'R', LDV >= max(1,N);
  125. *> if STOREV = 'R', LDV >= K.
  126. *> \endverbatim
  127. *>
  128. *> \param[in] T
  129. *> \verbatim
  130. *> T is DOUBLE PRECISION array, dimension (LDT,K)
  131. *> The triangular K-by-K matrix T in the representation of the
  132. *> block reflector.
  133. *> \endverbatim
  134. *>
  135. *> \param[in] LDT
  136. *> \verbatim
  137. *> LDT is INTEGER
  138. *> The leading dimension of the array T.
  139. *> LDT >= K.
  140. *> \endverbatim
  141. *>
  142. *> \param[in,out] A
  143. *> \verbatim
  144. *> A is DOUBLE PRECISION array, dimension
  145. *> (LDA,N) if SIDE = 'L' or (LDA,K) if SIDE = 'R'
  146. *> On entry, the K-by-N or M-by-K matrix A.
  147. *> On exit, A is overwritten by the corresponding block of
  148. *> H*C or H**T*C or C*H or C*H**T. See Further Details.
  149. *> \endverbatim
  150. *>
  151. *> \param[in] LDA
  152. *> \verbatim
  153. *> LDA is INTEGER
  154. *> The leading dimension of the array A.
  155. *> If SIDE = 'L', LDA >= max(1,K);
  156. *> If SIDE = 'R', LDA >= max(1,M).
  157. *> \endverbatim
  158. *>
  159. *> \param[in,out] B
  160. *> \verbatim
  161. *> B is DOUBLE PRECISION array, dimension (LDB,N)
  162. *> On entry, the M-by-N matrix B.
  163. *> On exit, B is overwritten by the corresponding block of
  164. *> H*C or H**T*C or C*H or C*H**T. See Further Details.
  165. *> \endverbatim
  166. *>
  167. *> \param[in] LDB
  168. *> \verbatim
  169. *> LDB is INTEGER
  170. *> The leading dimension of the array B.
  171. *> LDB >= max(1,M).
  172. *> \endverbatim
  173. *>
  174. *> \param[out] WORK
  175. *> \verbatim
  176. *> WORK is DOUBLE PRECISION array, dimension
  177. *> (LDWORK,N) if SIDE = 'L',
  178. *> (LDWORK,K) if SIDE = 'R'.
  179. *> \endverbatim
  180. *>
  181. *> \param[in] LDWORK
  182. *> \verbatim
  183. *> LDWORK is INTEGER
  184. *> The leading dimension of the array WORK.
  185. *> If SIDE = 'L', LDWORK >= K;
  186. *> if SIDE = 'R', LDWORK >= M.
  187. *> \endverbatim
  188. *
  189. * Authors:
  190. * ========
  191. *
  192. *> \author Univ. of Tennessee
  193. *> \author Univ. of California Berkeley
  194. *> \author Univ. of Colorado Denver
  195. *> \author NAG Ltd.
  196. *
  197. *> \ingroup doubleOTHERauxiliary
  198. *
  199. *> \par Further Details:
  200. * =====================
  201. *>
  202. *> \verbatim
  203. *>
  204. *> The matrix C is a composite matrix formed from blocks A and B.
  205. *> The block B is of size M-by-N; if SIDE = 'R', A is of size M-by-K,
  206. *> and if SIDE = 'L', A is of size K-by-N.
  207. *>
  208. *> If SIDE = 'R' and DIRECT = 'F', C = [A B].
  209. *>
  210. *> If SIDE = 'L' and DIRECT = 'F', C = [A]
  211. *> [B].
  212. *>
  213. *> If SIDE = 'R' and DIRECT = 'B', C = [B A].
  214. *>
  215. *> If SIDE = 'L' and DIRECT = 'B', C = [B]
  216. *> [A].
  217. *>
  218. *> The pentagonal matrix V is composed of a rectangular block V1 and a
  219. *> trapezoidal block V2. The size of the trapezoidal block is determined by
  220. *> the parameter L, where 0<=L<=K. If L=K, the V2 block of V is triangular;
  221. *> if L=0, there is no trapezoidal block, thus V = V1 is rectangular.
  222. *>
  223. *> If DIRECT = 'F' and STOREV = 'C': V = [V1]
  224. *> [V2]
  225. *> - V2 is upper trapezoidal (first L rows of K-by-K upper triangular)
  226. *>
  227. *> If DIRECT = 'F' and STOREV = 'R': V = [V1 V2]
  228. *>
  229. *> - V2 is lower trapezoidal (first L columns of K-by-K lower triangular)
  230. *>
  231. *> If DIRECT = 'B' and STOREV = 'C': V = [V2]
  232. *> [V1]
  233. *> - V2 is lower trapezoidal (last L rows of K-by-K lower triangular)
  234. *>
  235. *> If DIRECT = 'B' and STOREV = 'R': V = [V2 V1]
  236. *>
  237. *> - V2 is upper trapezoidal (last L columns of K-by-K upper triangular)
  238. *>
  239. *> If STOREV = 'C' and SIDE = 'L', V is M-by-K with V2 L-by-K.
  240. *>
  241. *> If STOREV = 'C' and SIDE = 'R', V is N-by-K with V2 L-by-K.
  242. *>
  243. *> If STOREV = 'R' and SIDE = 'L', V is K-by-M with V2 K-by-L.
  244. *>
  245. *> If STOREV = 'R' and SIDE = 'R', V is K-by-N with V2 K-by-L.
  246. *> \endverbatim
  247. *>
  248. * =====================================================================
  249. SUBROUTINE DTPRFB( SIDE, TRANS, DIRECT, STOREV, M, N, K, L,
  250. $ V, LDV, T, LDT, A, LDA, B, LDB, WORK, LDWORK )
  251. *
  252. * -- LAPACK auxiliary routine --
  253. * -- LAPACK is a software package provided by Univ. of Tennessee, --
  254. * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  255. *
  256. * .. Scalar Arguments ..
  257. CHARACTER DIRECT, SIDE, STOREV, TRANS
  258. INTEGER K, L, LDA, LDB, LDT, LDV, LDWORK, M, N
  259. * ..
  260. * .. Array Arguments ..
  261. DOUBLE PRECISION A( LDA, * ), B( LDB, * ), T( LDT, * ),
  262. $ V( LDV, * ), WORK( LDWORK, * )
  263. * ..
  264. *
  265. * ==========================================================================
  266. *
  267. * .. Parameters ..
  268. DOUBLE PRECISION ONE, ZERO
  269. PARAMETER ( ONE = 1.0, ZERO = 0.0 )
  270. * ..
  271. * .. Local Scalars ..
  272. INTEGER I, J, MP, NP, KP
  273. LOGICAL LEFT, FORWARD, COLUMN, RIGHT, BACKWARD, ROW
  274. * ..
  275. * .. External Functions ..
  276. LOGICAL LSAME
  277. EXTERNAL LSAME
  278. * ..
  279. * .. External Subroutines ..
  280. EXTERNAL DGEMM, DTRMM
  281. * ..
  282. * .. Executable Statements ..
  283. *
  284. * Quick return if possible
  285. *
  286. IF( M.LE.0 .OR. N.LE.0 .OR. K.LE.0 .OR. L.LT.0 ) RETURN
  287. *
  288. IF( LSAME( STOREV, 'C' ) ) THEN
  289. COLUMN = .TRUE.
  290. ROW = .FALSE.
  291. ELSE IF ( LSAME( STOREV, 'R' ) ) THEN
  292. COLUMN = .FALSE.
  293. ROW = .TRUE.
  294. ELSE
  295. COLUMN = .FALSE.
  296. ROW = .FALSE.
  297. END IF
  298. *
  299. IF( LSAME( SIDE, 'L' ) ) THEN
  300. LEFT = .TRUE.
  301. RIGHT = .FALSE.
  302. ELSE IF( LSAME( SIDE, 'R' ) ) THEN
  303. LEFT = .FALSE.
  304. RIGHT = .TRUE.
  305. ELSE
  306. LEFT = .FALSE.
  307. RIGHT = .FALSE.
  308. END IF
  309. *
  310. IF( LSAME( DIRECT, 'F' ) ) THEN
  311. FORWARD = .TRUE.
  312. BACKWARD = .FALSE.
  313. ELSE IF( LSAME( DIRECT, 'B' ) ) THEN
  314. FORWARD = .FALSE.
  315. BACKWARD = .TRUE.
  316. ELSE
  317. FORWARD = .FALSE.
  318. BACKWARD = .FALSE.
  319. END IF
  320. *
  321. * ---------------------------------------------------------------------------
  322. *
  323. IF( COLUMN .AND. FORWARD .AND. LEFT ) THEN
  324. *
  325. * ---------------------------------------------------------------------------
  326. *
  327. * Let W = [ I ] (K-by-K)
  328. * [ V ] (M-by-K)
  329. *
  330. * Form H C or H**T C where C = [ A ] (K-by-N)
  331. * [ B ] (M-by-N)
  332. *
  333. * H = I - W T W**T or H**T = I - W T**T W**T
  334. *
  335. * A = A - T (A + V**T B) or A = A - T**T (A + V**T B)
  336. * B = B - V T (A + V**T B) or B = B - V T**T (A + V**T B)
  337. *
  338. * ---------------------------------------------------------------------------
  339. *
  340. MP = MIN( M-L+1, M )
  341. KP = MIN( L+1, K )
  342. *
  343. DO J = 1, N
  344. DO I = 1, L
  345. WORK( I, J ) = B( M-L+I, J )
  346. END DO
  347. END DO
  348. CALL DTRMM( 'L', 'U', 'T', 'N', L, N, ONE, V( MP, 1 ), LDV,
  349. $ WORK, LDWORK )
  350. CALL DGEMM( 'T', 'N', L, N, M-L, ONE, V, LDV, B, LDB,
  351. $ ONE, WORK, LDWORK )
  352. CALL DGEMM( 'T', 'N', K-L, N, M, ONE, V( 1, KP ), LDV,
  353. $ B, LDB, ZERO, WORK( KP, 1 ), LDWORK )
  354. *
  355. DO J = 1, N
  356. DO I = 1, K
  357. WORK( I, J ) = WORK( I, J ) + A( I, J )
  358. END DO
  359. END DO
  360. *
  361. CALL DTRMM( 'L', 'U', TRANS, 'N', K, N, ONE, T, LDT,
  362. $ WORK, LDWORK )
  363. *
  364. DO J = 1, N
  365. DO I = 1, K
  366. A( I, J ) = A( I, J ) - WORK( I, J )
  367. END DO
  368. END DO
  369. *
  370. CALL DGEMM( 'N', 'N', M-L, N, K, -ONE, V, LDV, WORK, LDWORK,
  371. $ ONE, B, LDB )
  372. CALL DGEMM( 'N', 'N', L, N, K-L, -ONE, V( MP, KP ), LDV,
  373. $ WORK( KP, 1 ), LDWORK, ONE, B( MP, 1 ), LDB )
  374. CALL DTRMM( 'L', 'U', 'N', 'N', L, N, ONE, V( MP, 1 ), LDV,
  375. $ WORK, LDWORK )
  376. DO J = 1, N
  377. DO I = 1, L
  378. B( M-L+I, J ) = B( M-L+I, J ) - WORK( I, J )
  379. END DO
  380. END DO
  381. *
  382. * ---------------------------------------------------------------------------
  383. *
  384. ELSE IF( COLUMN .AND. FORWARD .AND. RIGHT ) THEN
  385. *
  386. * ---------------------------------------------------------------------------
  387. *
  388. * Let W = [ I ] (K-by-K)
  389. * [ V ] (N-by-K)
  390. *
  391. * Form C H or C H**T where C = [ A B ] (A is M-by-K, B is M-by-N)
  392. *
  393. * H = I - W T W**T or H**T = I - W T**T W**T
  394. *
  395. * A = A - (A + B V) T or A = A - (A + B V) T**T
  396. * B = B - (A + B V) T V**T or B = B - (A + B V) T**T V**T
  397. *
  398. * ---------------------------------------------------------------------------
  399. *
  400. NP = MIN( N-L+1, N )
  401. KP = MIN( L+1, K )
  402. *
  403. DO J = 1, L
  404. DO I = 1, M
  405. WORK( I, J ) = B( I, N-L+J )
  406. END DO
  407. END DO
  408. CALL DTRMM( 'R', 'U', 'N', 'N', M, L, ONE, V( NP, 1 ), LDV,
  409. $ WORK, LDWORK )
  410. CALL DGEMM( 'N', 'N', M, L, N-L, ONE, B, LDB,
  411. $ V, LDV, ONE, WORK, LDWORK )
  412. CALL DGEMM( 'N', 'N', M, K-L, N, ONE, B, LDB,
  413. $ V( 1, KP ), LDV, ZERO, WORK( 1, KP ), LDWORK )
  414. *
  415. DO J = 1, K
  416. DO I = 1, M
  417. WORK( I, J ) = WORK( I, J ) + A( I, J )
  418. END DO
  419. END DO
  420. *
  421. CALL DTRMM( 'R', 'U', TRANS, 'N', M, K, ONE, T, LDT,
  422. $ WORK, LDWORK )
  423. *
  424. DO J = 1, K
  425. DO I = 1, M
  426. A( I, J ) = A( I, J ) - WORK( I, J )
  427. END DO
  428. END DO
  429. *
  430. CALL DGEMM( 'N', 'T', M, N-L, K, -ONE, WORK, LDWORK,
  431. $ V, LDV, ONE, B, LDB )
  432. CALL DGEMM( 'N', 'T', M, L, K-L, -ONE, WORK( 1, KP ), LDWORK,
  433. $ V( NP, KP ), LDV, ONE, B( 1, NP ), LDB )
  434. CALL DTRMM( 'R', 'U', 'T', 'N', M, L, ONE, V( NP, 1 ), LDV,
  435. $ WORK, LDWORK )
  436. DO J = 1, L
  437. DO I = 1, M
  438. B( I, N-L+J ) = B( I, N-L+J ) - WORK( I, J )
  439. END DO
  440. END DO
  441. *
  442. * ---------------------------------------------------------------------------
  443. *
  444. ELSE IF( COLUMN .AND. BACKWARD .AND. LEFT ) THEN
  445. *
  446. * ---------------------------------------------------------------------------
  447. *
  448. * Let W = [ V ] (M-by-K)
  449. * [ I ] (K-by-K)
  450. *
  451. * Form H C or H**T C where C = [ B ] (M-by-N)
  452. * [ A ] (K-by-N)
  453. *
  454. * H = I - W T W**T or H**T = I - W T**T W**T
  455. *
  456. * A = A - T (A + V**T B) or A = A - T**T (A + V**T B)
  457. * B = B - V T (A + V**T B) or B = B - V T**T (A + V**T B)
  458. *
  459. * ---------------------------------------------------------------------------
  460. *
  461. MP = MIN( L+1, M )
  462. KP = MIN( K-L+1, K )
  463. *
  464. DO J = 1, N
  465. DO I = 1, L
  466. WORK( K-L+I, J ) = B( I, J )
  467. END DO
  468. END DO
  469. *
  470. CALL DTRMM( 'L', 'L', 'T', 'N', L, N, ONE, V( 1, KP ), LDV,
  471. $ WORK( KP, 1 ), LDWORK )
  472. CALL DGEMM( 'T', 'N', L, N, M-L, ONE, V( MP, KP ), LDV,
  473. $ B( MP, 1 ), LDB, ONE, WORK( KP, 1 ), LDWORK )
  474. CALL DGEMM( 'T', 'N', K-L, N, M, ONE, V, LDV,
  475. $ B, LDB, ZERO, WORK, LDWORK )
  476. *
  477. DO J = 1, N
  478. DO I = 1, K
  479. WORK( I, J ) = WORK( I, J ) + A( I, J )
  480. END DO
  481. END DO
  482. *
  483. CALL DTRMM( 'L', 'L', TRANS, 'N', K, N, ONE, T, LDT,
  484. $ WORK, LDWORK )
  485. *
  486. DO J = 1, N
  487. DO I = 1, K
  488. A( I, J ) = A( I, J ) - WORK( I, J )
  489. END DO
  490. END DO
  491. *
  492. CALL DGEMM( 'N', 'N', M-L, N, K, -ONE, V( MP, 1 ), LDV,
  493. $ WORK, LDWORK, ONE, B( MP, 1 ), LDB )
  494. CALL DGEMM( 'N', 'N', L, N, K-L, -ONE, V, LDV,
  495. $ WORK, LDWORK, ONE, B, LDB )
  496. CALL DTRMM( 'L', 'L', 'N', 'N', L, N, ONE, V( 1, KP ), LDV,
  497. $ WORK( KP, 1 ), LDWORK )
  498. DO J = 1, N
  499. DO I = 1, L
  500. B( I, J ) = B( I, J ) - WORK( K-L+I, J )
  501. END DO
  502. END DO
  503. *
  504. * ---------------------------------------------------------------------------
  505. *
  506. ELSE IF( COLUMN .AND. BACKWARD .AND. RIGHT ) THEN
  507. *
  508. * ---------------------------------------------------------------------------
  509. *
  510. * Let W = [ V ] (N-by-K)
  511. * [ I ] (K-by-K)
  512. *
  513. * Form C H or C H**T where C = [ B A ] (B is M-by-N, A is M-by-K)
  514. *
  515. * H = I - W T W**T or H**T = I - W T**T W**T
  516. *
  517. * A = A - (A + B V) T or A = A - (A + B V) T**T
  518. * B = B - (A + B V) T V**T or B = B - (A + B V) T**T V**T
  519. *
  520. * ---------------------------------------------------------------------------
  521. *
  522. NP = MIN( L+1, N )
  523. KP = MIN( K-L+1, K )
  524. *
  525. DO J = 1, L
  526. DO I = 1, M
  527. WORK( I, K-L+J ) = B( I, J )
  528. END DO
  529. END DO
  530. CALL DTRMM( 'R', 'L', 'N', 'N', M, L, ONE, V( 1, KP ), LDV,
  531. $ WORK( 1, KP ), LDWORK )
  532. CALL DGEMM( 'N', 'N', M, L, N-L, ONE, B( 1, NP ), LDB,
  533. $ V( NP, KP ), LDV, ONE, WORK( 1, KP ), LDWORK )
  534. CALL DGEMM( 'N', 'N', M, K-L, N, ONE, B, LDB,
  535. $ V, LDV, ZERO, WORK, LDWORK )
  536. *
  537. DO J = 1, K
  538. DO I = 1, M
  539. WORK( I, J ) = WORK( I, J ) + A( I, J )
  540. END DO
  541. END DO
  542. *
  543. CALL DTRMM( 'R', 'L', TRANS, 'N', M, K, ONE, T, LDT,
  544. $ WORK, LDWORK )
  545. *
  546. DO J = 1, K
  547. DO I = 1, M
  548. A( I, J ) = A( I, J ) - WORK( I, J )
  549. END DO
  550. END DO
  551. *
  552. CALL DGEMM( 'N', 'T', M, N-L, K, -ONE, WORK, LDWORK,
  553. $ V( NP, 1 ), LDV, ONE, B( 1, NP ), LDB )
  554. CALL DGEMM( 'N', 'T', M, L, K-L, -ONE, WORK, LDWORK,
  555. $ V, LDV, ONE, B, LDB )
  556. CALL DTRMM( 'R', 'L', 'T', 'N', M, L, ONE, V( 1, KP ), LDV,
  557. $ WORK( 1, KP ), LDWORK )
  558. DO J = 1, L
  559. DO I = 1, M
  560. B( I, J ) = B( I, J ) - WORK( I, K-L+J )
  561. END DO
  562. END DO
  563. *
  564. * ---------------------------------------------------------------------------
  565. *
  566. ELSE IF( ROW .AND. FORWARD .AND. LEFT ) THEN
  567. *
  568. * ---------------------------------------------------------------------------
  569. *
  570. * Let W = [ I V ] ( I is K-by-K, V is K-by-M )
  571. *
  572. * Form H C or H**T C where C = [ A ] (K-by-N)
  573. * [ B ] (M-by-N)
  574. *
  575. * H = I - W**T T W or H**T = I - W**T T**T W
  576. *
  577. * A = A - T (A + V B) or A = A - T**T (A + V B)
  578. * B = B - V**T T (A + V B) or B = B - V**T T**T (A + V B)
  579. *
  580. * ---------------------------------------------------------------------------
  581. *
  582. MP = MIN( M-L+1, M )
  583. KP = MIN( L+1, K )
  584. *
  585. DO J = 1, N
  586. DO I = 1, L
  587. WORK( I, J ) = B( M-L+I, J )
  588. END DO
  589. END DO
  590. CALL DTRMM( 'L', 'L', 'N', 'N', L, N, ONE, V( 1, MP ), LDV,
  591. $ WORK, LDB )
  592. CALL DGEMM( 'N', 'N', L, N, M-L, ONE, V, LDV,B, LDB,
  593. $ ONE, WORK, LDWORK )
  594. CALL DGEMM( 'N', 'N', K-L, N, M, ONE, V( KP, 1 ), LDV,
  595. $ B, LDB, ZERO, WORK( KP, 1 ), LDWORK )
  596. *
  597. DO J = 1, N
  598. DO I = 1, K
  599. WORK( I, J ) = WORK( I, J ) + A( I, J )
  600. END DO
  601. END DO
  602. *
  603. CALL DTRMM( 'L', 'U', TRANS, 'N', K, N, ONE, T, LDT,
  604. $ WORK, LDWORK )
  605. *
  606. DO J = 1, N
  607. DO I = 1, K
  608. A( I, J ) = A( I, J ) - WORK( I, J )
  609. END DO
  610. END DO
  611. *
  612. CALL DGEMM( 'T', 'N', M-L, N, K, -ONE, V, LDV, WORK, LDWORK,
  613. $ ONE, B, LDB )
  614. CALL DGEMM( 'T', 'N', L, N, K-L, -ONE, V( KP, MP ), LDV,
  615. $ WORK( KP, 1 ), LDWORK, ONE, B( MP, 1 ), LDB )
  616. CALL DTRMM( 'L', 'L', 'T', 'N', L, N, ONE, V( 1, MP ), LDV,
  617. $ WORK, LDWORK )
  618. DO J = 1, N
  619. DO I = 1, L
  620. B( M-L+I, J ) = B( M-L+I, J ) - WORK( I, J )
  621. END DO
  622. END DO
  623. *
  624. * ---------------------------------------------------------------------------
  625. *
  626. ELSE IF( ROW .AND. FORWARD .AND. RIGHT ) THEN
  627. *
  628. * ---------------------------------------------------------------------------
  629. *
  630. * Let W = [ I V ] ( I is K-by-K, V is K-by-N )
  631. *
  632. * Form C H or C H**T where C = [ A B ] (A is M-by-K, B is M-by-N)
  633. *
  634. * H = I - W**T T W or H**T = I - W**T T**T W
  635. *
  636. * A = A - (A + B V**T) T or A = A - (A + B V**T) T**T
  637. * B = B - (A + B V**T) T V or B = B - (A + B V**T) T**T V
  638. *
  639. * ---------------------------------------------------------------------------
  640. *
  641. NP = MIN( N-L+1, N )
  642. KP = MIN( L+1, K )
  643. *
  644. DO J = 1, L
  645. DO I = 1, M
  646. WORK( I, J ) = B( I, N-L+J )
  647. END DO
  648. END DO
  649. CALL DTRMM( 'R', 'L', 'T', 'N', M, L, ONE, V( 1, NP ), LDV,
  650. $ WORK, LDWORK )
  651. CALL DGEMM( 'N', 'T', M, L, N-L, ONE, B, LDB, V, LDV,
  652. $ ONE, WORK, LDWORK )
  653. CALL DGEMM( 'N', 'T', M, K-L, N, ONE, B, LDB,
  654. $ V( KP, 1 ), LDV, ZERO, WORK( 1, KP ), LDWORK )
  655. *
  656. DO J = 1, K
  657. DO I = 1, M
  658. WORK( I, J ) = WORK( I, J ) + A( I, J )
  659. END DO
  660. END DO
  661. *
  662. CALL DTRMM( 'R', 'U', TRANS, 'N', M, K, ONE, T, LDT,
  663. $ WORK, LDWORK )
  664. *
  665. DO J = 1, K
  666. DO I = 1, M
  667. A( I, J ) = A( I, J ) - WORK( I, J )
  668. END DO
  669. END DO
  670. *
  671. CALL DGEMM( 'N', 'N', M, N-L, K, -ONE, WORK, LDWORK,
  672. $ V, LDV, ONE, B, LDB )
  673. CALL DGEMM( 'N', 'N', M, L, K-L, -ONE, WORK( 1, KP ), LDWORK,
  674. $ V( KP, NP ), LDV, ONE, B( 1, NP ), LDB )
  675. CALL DTRMM( 'R', 'L', 'N', 'N', M, L, ONE, V( 1, NP ), LDV,
  676. $ WORK, LDWORK )
  677. DO J = 1, L
  678. DO I = 1, M
  679. B( I, N-L+J ) = B( I, N-L+J ) - WORK( I, J )
  680. END DO
  681. END DO
  682. *
  683. * ---------------------------------------------------------------------------
  684. *
  685. ELSE IF( ROW .AND. BACKWARD .AND. LEFT ) THEN
  686. *
  687. * ---------------------------------------------------------------------------
  688. *
  689. * Let W = [ V I ] ( I is K-by-K, V is K-by-M )
  690. *
  691. * Form H C or H**T C where C = [ B ] (M-by-N)
  692. * [ A ] (K-by-N)
  693. *
  694. * H = I - W**T T W or H**T = I - W**T T**T W
  695. *
  696. * A = A - T (A + V B) or A = A - T**T (A + V B)
  697. * B = B - V**T T (A + V B) or B = B - V**T T**T (A + V B)
  698. *
  699. * ---------------------------------------------------------------------------
  700. *
  701. MP = MIN( L+1, M )
  702. KP = MIN( K-L+1, K )
  703. *
  704. DO J = 1, N
  705. DO I = 1, L
  706. WORK( K-L+I, J ) = B( I, J )
  707. END DO
  708. END DO
  709. CALL DTRMM( 'L', 'U', 'N', 'N', L, N, ONE, V( KP, 1 ), LDV,
  710. $ WORK( KP, 1 ), LDWORK )
  711. CALL DGEMM( 'N', 'N', L, N, M-L, ONE, V( KP, MP ), LDV,
  712. $ B( MP, 1 ), LDB, ONE, WORK( KP, 1 ), LDWORK )
  713. CALL DGEMM( 'N', 'N', K-L, N, M, ONE, V, LDV, B, LDB,
  714. $ ZERO, WORK, LDWORK )
  715. *
  716. DO J = 1, N
  717. DO I = 1, K
  718. WORK( I, J ) = WORK( I, J ) + A( I, J )
  719. END DO
  720. END DO
  721. *
  722. CALL DTRMM( 'L', 'L ', TRANS, 'N', K, N, ONE, T, LDT,
  723. $ WORK, LDWORK )
  724. *
  725. DO J = 1, N
  726. DO I = 1, K
  727. A( I, J ) = A( I, J ) - WORK( I, J )
  728. END DO
  729. END DO
  730. *
  731. CALL DGEMM( 'T', 'N', M-L, N, K, -ONE, V( 1, MP ), LDV,
  732. $ WORK, LDWORK, ONE, B( MP, 1 ), LDB )
  733. CALL DGEMM( 'T', 'N', L, N, K-L, -ONE, V, LDV,
  734. $ WORK, LDWORK, ONE, B, LDB )
  735. CALL DTRMM( 'L', 'U', 'T', 'N', L, N, ONE, V( KP, 1 ), LDV,
  736. $ WORK( KP, 1 ), LDWORK )
  737. DO J = 1, N
  738. DO I = 1, L
  739. B( I, J ) = B( I, J ) - WORK( K-L+I, J )
  740. END DO
  741. END DO
  742. *
  743. * ---------------------------------------------------------------------------
  744. *
  745. ELSE IF( ROW .AND. BACKWARD .AND. RIGHT ) THEN
  746. *
  747. * ---------------------------------------------------------------------------
  748. *
  749. * Let W = [ V I ] ( I is K-by-K, V is K-by-N )
  750. *
  751. * Form C H or C H**T where C = [ B A ] (A is M-by-K, B is M-by-N)
  752. *
  753. * H = I - W**T T W or H**T = I - W**T T**T W
  754. *
  755. * A = A - (A + B V**T) T or A = A - (A + B V**T) T**T
  756. * B = B - (A + B V**T) T V or B = B - (A + B V**T) T**T V
  757. *
  758. * ---------------------------------------------------------------------------
  759. *
  760. NP = MIN( L+1, N )
  761. KP = MIN( K-L+1, K )
  762. *
  763. DO J = 1, L
  764. DO I = 1, M
  765. WORK( I, K-L+J ) = B( I, J )
  766. END DO
  767. END DO
  768. CALL DTRMM( 'R', 'U', 'T', 'N', M, L, ONE, V( KP, 1 ), LDV,
  769. $ WORK( 1, KP ), LDWORK )
  770. CALL DGEMM( 'N', 'T', M, L, N-L, ONE, B( 1, NP ), LDB,
  771. $ V( KP, NP ), LDV, ONE, WORK( 1, KP ), LDWORK )
  772. CALL DGEMM( 'N', 'T', M, K-L, N, ONE, B, LDB, V, LDV,
  773. $ ZERO, WORK, LDWORK )
  774. *
  775. DO J = 1, K
  776. DO I = 1, M
  777. WORK( I, J ) = WORK( I, J ) + A( I, J )
  778. END DO
  779. END DO
  780. *
  781. CALL DTRMM( 'R', 'L', TRANS, 'N', M, K, ONE, T, LDT,
  782. $ WORK, LDWORK )
  783. *
  784. DO J = 1, K
  785. DO I = 1, M
  786. A( I, J ) = A( I, J ) - WORK( I, J )
  787. END DO
  788. END DO
  789. *
  790. CALL DGEMM( 'N', 'N', M, N-L, K, -ONE, WORK, LDWORK,
  791. $ V( 1, NP ), LDV, ONE, B( 1, NP ), LDB )
  792. CALL DGEMM( 'N', 'N', M, L, K-L , -ONE, WORK, LDWORK,
  793. $ V, LDV, ONE, B, LDB )
  794. CALL DTRMM( 'R', 'U', 'N', 'N', M, L, ONE, V( KP, 1 ), LDV,
  795. $ WORK( 1, KP ), LDWORK )
  796. DO J = 1, L
  797. DO I = 1, M
  798. B( I, J ) = B( I, J ) - WORK( I, K-L+J )
  799. END DO
  800. END DO
  801. *
  802. END IF
  803. *
  804. RETURN
  805. *
  806. * End of DTPRFB
  807. *
  808. END