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.

ctfttp.f 16 kB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540
  1. *> \brief \b CTFTTP copies a triangular matrix from the rectangular full packed format (TF) to the standard packed format (TP).
  2. *
  3. * =========== DOCUMENTATION ===========
  4. *
  5. * Online html documentation available at
  6. * http://www.netlib.org/lapack/explore-html/
  7. *
  8. *> \htmlonly
  9. *> Download CTFTTP + dependencies
  10. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ctfttp.f">
  11. *> [TGZ]</a>
  12. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ctfttp.f">
  13. *> [ZIP]</a>
  14. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ctfttp.f">
  15. *> [TXT]</a>
  16. *> \endhtmlonly
  17. *
  18. * Definition:
  19. * ===========
  20. *
  21. * SUBROUTINE CTFTTP( TRANSR, UPLO, N, ARF, AP, INFO )
  22. *
  23. * .. Scalar Arguments ..
  24. * CHARACTER TRANSR, UPLO
  25. * INTEGER INFO, N
  26. * ..
  27. * .. Array Arguments ..
  28. * COMPLEX AP( 0: * ), ARF( 0: * )
  29. * ..
  30. *
  31. *
  32. *> \par Purpose:
  33. * =============
  34. *>
  35. *> \verbatim
  36. *>
  37. *> CTFTTP copies a triangular matrix A from rectangular full packed
  38. *> format (TF) to standard packed format (TP).
  39. *> \endverbatim
  40. *
  41. * Arguments:
  42. * ==========
  43. *
  44. *> \param[in] TRANSR
  45. *> \verbatim
  46. *> TRANSR is CHARACTER*1
  47. *> = 'N': ARF is in Normal format;
  48. *> = 'C': ARF is in Conjugate-transpose format;
  49. *> \endverbatim
  50. *>
  51. *> \param[in] UPLO
  52. *> \verbatim
  53. *> UPLO is CHARACTER*1
  54. *> = 'U': A is upper triangular;
  55. *> = 'L': A is lower triangular.
  56. *> \endverbatim
  57. *>
  58. *> \param[in] N
  59. *> \verbatim
  60. *> N is INTEGER
  61. *> The order of the matrix A. N >= 0.
  62. *> \endverbatim
  63. *>
  64. *> \param[in] ARF
  65. *> \verbatim
  66. *> ARF is COMPLEX array, dimension ( N*(N+1)/2 ),
  67. *> On entry, the upper or lower triangular matrix A stored in
  68. *> RFP format. For a further discussion see Notes below.
  69. *> \endverbatim
  70. *>
  71. *> \param[out] AP
  72. *> \verbatim
  73. *> AP is COMPLEX array, dimension ( N*(N+1)/2 ),
  74. *> On exit, the upper or lower triangular matrix A, packed
  75. *> columnwise in a linear array. The j-th column of A is stored
  76. *> in the array AP as follows:
  77. *> if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
  78. *> if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.
  79. *> \endverbatim
  80. *>
  81. *> \param[out] INFO
  82. *> \verbatim
  83. *> INFO is INTEGER
  84. *> = 0: successful exit
  85. *> < 0: if INFO = -i, the i-th argument had an illegal value
  86. *> \endverbatim
  87. *
  88. * Authors:
  89. * ========
  90. *
  91. *> \author Univ. of Tennessee
  92. *> \author Univ. of California Berkeley
  93. *> \author Univ. of Colorado Denver
  94. *> \author NAG Ltd.
  95. *
  96. *> \ingroup complexOTHERcomputational
  97. *
  98. *> \par Further Details:
  99. * =====================
  100. *>
  101. *> \verbatim
  102. *>
  103. *> We first consider Standard Packed Format when N is even.
  104. *> We give an example where N = 6.
  105. *>
  106. *> AP is Upper AP is Lower
  107. *>
  108. *> 00 01 02 03 04 05 00
  109. *> 11 12 13 14 15 10 11
  110. *> 22 23 24 25 20 21 22
  111. *> 33 34 35 30 31 32 33
  112. *> 44 45 40 41 42 43 44
  113. *> 55 50 51 52 53 54 55
  114. *>
  115. *>
  116. *> Let TRANSR = 'N'. RFP holds AP as follows:
  117. *> For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last
  118. *> three columns of AP upper. The lower triangle A(4:6,0:2) consists of
  119. *> conjugate-transpose of the first three columns of AP upper.
  120. *> For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first
  121. *> three columns of AP lower. The upper triangle A(0:2,0:2) consists of
  122. *> conjugate-transpose of the last three columns of AP lower.
  123. *> To denote conjugate we place -- above the element. This covers the
  124. *> case N even and TRANSR = 'N'.
  125. *>
  126. *> RFP A RFP A
  127. *>
  128. *> -- -- --
  129. *> 03 04 05 33 43 53
  130. *> -- --
  131. *> 13 14 15 00 44 54
  132. *> --
  133. *> 23 24 25 10 11 55
  134. *>
  135. *> 33 34 35 20 21 22
  136. *> --
  137. *> 00 44 45 30 31 32
  138. *> -- --
  139. *> 01 11 55 40 41 42
  140. *> -- -- --
  141. *> 02 12 22 50 51 52
  142. *>
  143. *> Now let TRANSR = 'C'. RFP A in both UPLO cases is just the conjugate-
  144. *> transpose of RFP A above. One therefore gets:
  145. *>
  146. *>
  147. *> RFP A RFP A
  148. *>
  149. *> -- -- -- -- -- -- -- -- -- --
  150. *> 03 13 23 33 00 01 02 33 00 10 20 30 40 50
  151. *> -- -- -- -- -- -- -- -- -- --
  152. *> 04 14 24 34 44 11 12 43 44 11 21 31 41 51
  153. *> -- -- -- -- -- -- -- -- -- --
  154. *> 05 15 25 35 45 55 22 53 54 55 22 32 42 52
  155. *>
  156. *>
  157. *> We next consider Standard Packed Format when N is odd.
  158. *> We give an example where N = 5.
  159. *>
  160. *> AP is Upper AP is Lower
  161. *>
  162. *> 00 01 02 03 04 00
  163. *> 11 12 13 14 10 11
  164. *> 22 23 24 20 21 22
  165. *> 33 34 30 31 32 33
  166. *> 44 40 41 42 43 44
  167. *>
  168. *>
  169. *> Let TRANSR = 'N'. RFP holds AP as follows:
  170. *> For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last
  171. *> three columns of AP upper. The lower triangle A(3:4,0:1) consists of
  172. *> conjugate-transpose of the first two columns of AP upper.
  173. *> For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first
  174. *> three columns of AP lower. The upper triangle A(0:1,1:2) consists of
  175. *> conjugate-transpose of the last two columns of AP lower.
  176. *> To denote conjugate we place -- above the element. This covers the
  177. *> case N odd and TRANSR = 'N'.
  178. *>
  179. *> RFP A RFP A
  180. *>
  181. *> -- --
  182. *> 02 03 04 00 33 43
  183. *> --
  184. *> 12 13 14 10 11 44
  185. *>
  186. *> 22 23 24 20 21 22
  187. *> --
  188. *> 00 33 34 30 31 32
  189. *> -- --
  190. *> 01 11 44 40 41 42
  191. *>
  192. *> Now let TRANSR = 'C'. RFP A in both UPLO cases is just the conjugate-
  193. *> transpose of RFP A above. One therefore gets:
  194. *>
  195. *>
  196. *> RFP A RFP A
  197. *>
  198. *> -- -- -- -- -- -- -- -- --
  199. *> 02 12 22 00 01 00 10 20 30 40 50
  200. *> -- -- -- -- -- -- -- -- --
  201. *> 03 13 23 33 11 33 11 21 31 41 51
  202. *> -- -- -- -- -- -- -- -- --
  203. *> 04 14 24 34 44 43 44 22 32 42 52
  204. *> \endverbatim
  205. *>
  206. * =====================================================================
  207. SUBROUTINE CTFTTP( TRANSR, UPLO, N, ARF, AP, INFO )
  208. *
  209. * -- LAPACK computational routine --
  210. * -- LAPACK is a software package provided by Univ. of Tennessee, --
  211. * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  212. *
  213. * .. Scalar Arguments ..
  214. CHARACTER TRANSR, UPLO
  215. INTEGER INFO, N
  216. * ..
  217. * .. Array Arguments ..
  218. COMPLEX AP( 0: * ), ARF( 0: * )
  219. * ..
  220. *
  221. * =====================================================================
  222. *
  223. * .. Parameters ..
  224. * ..
  225. * .. Local Scalars ..
  226. LOGICAL LOWER, NISODD, NORMALTRANSR
  227. INTEGER N1, N2, K, NT
  228. INTEGER I, J, IJ
  229. INTEGER IJP, JP, LDA, JS
  230. * ..
  231. * .. External Functions ..
  232. LOGICAL LSAME
  233. EXTERNAL LSAME
  234. * ..
  235. * .. External Subroutines ..
  236. EXTERNAL XERBLA
  237. * ..
  238. * .. Intrinsic Functions ..
  239. INTRINSIC CONJG
  240. * ..
  241. * .. Intrinsic Functions ..
  242. * ..
  243. * .. Executable Statements ..
  244. *
  245. * Test the input parameters.
  246. *
  247. INFO = 0
  248. NORMALTRANSR = LSAME( TRANSR, 'N' )
  249. LOWER = LSAME( UPLO, 'L' )
  250. IF( .NOT.NORMALTRANSR .AND. .NOT.LSAME( TRANSR, 'C' ) ) THEN
  251. INFO = -1
  252. ELSE IF( .NOT.LOWER .AND. .NOT.LSAME( UPLO, 'U' ) ) THEN
  253. INFO = -2
  254. ELSE IF( N.LT.0 ) THEN
  255. INFO = -3
  256. END IF
  257. IF( INFO.NE.0 ) THEN
  258. CALL XERBLA( 'CTFTTP', -INFO )
  259. RETURN
  260. END IF
  261. *
  262. * Quick return if possible
  263. *
  264. IF( N.EQ.0 )
  265. $ RETURN
  266. *
  267. IF( N.EQ.1 ) THEN
  268. IF( NORMALTRANSR ) THEN
  269. AP( 0 ) = ARF( 0 )
  270. ELSE
  271. AP( 0 ) = CONJG( ARF( 0 ) )
  272. END IF
  273. RETURN
  274. END IF
  275. *
  276. * Size of array ARF(0:NT-1)
  277. *
  278. NT = N*( N+1 ) / 2
  279. *
  280. * Set N1 and N2 depending on LOWER
  281. *
  282. IF( LOWER ) THEN
  283. N2 = N / 2
  284. N1 = N - N2
  285. ELSE
  286. N1 = N / 2
  287. N2 = N - N1
  288. END IF
  289. *
  290. * If N is odd, set NISODD = .TRUE.
  291. * If N is even, set K = N/2 and NISODD = .FALSE.
  292. *
  293. * set lda of ARF^C; ARF^C is (0:(N+1)/2-1,0:N-noe)
  294. * where noe = 0 if n is even, noe = 1 if n is odd
  295. *
  296. IF( MOD( N, 2 ).EQ.0 ) THEN
  297. K = N / 2
  298. NISODD = .FALSE.
  299. LDA = N + 1
  300. ELSE
  301. NISODD = .TRUE.
  302. LDA = N
  303. END IF
  304. *
  305. * ARF^C has lda rows and n+1-noe cols
  306. *
  307. IF( .NOT.NORMALTRANSR )
  308. $ LDA = ( N+1 ) / 2
  309. *
  310. * start execution: there are eight cases
  311. *
  312. IF( NISODD ) THEN
  313. *
  314. * N is odd
  315. *
  316. IF( NORMALTRANSR ) THEN
  317. *
  318. * N is odd and TRANSR = 'N'
  319. *
  320. IF( LOWER ) THEN
  321. *
  322. * SRPA for LOWER, NORMAL and N is odd ( a(0:n-1,0:n1-1) )
  323. * T1 -> a(0,0), T2 -> a(0,1), S -> a(n1,0)
  324. * T1 -> a(0), T2 -> a(n), S -> a(n1); lda = n
  325. *
  326. IJP = 0
  327. JP = 0
  328. DO J = 0, N2
  329. DO I = J, N - 1
  330. IJ = I + JP
  331. AP( IJP ) = ARF( IJ )
  332. IJP = IJP + 1
  333. END DO
  334. JP = JP + LDA
  335. END DO
  336. DO I = 0, N2 - 1
  337. DO J = 1 + I, N2
  338. IJ = I + J*LDA
  339. AP( IJP ) = CONJG( ARF( IJ ) )
  340. IJP = IJP + 1
  341. END DO
  342. END DO
  343. *
  344. ELSE
  345. *
  346. * SRPA for UPPER, NORMAL and N is odd ( a(0:n-1,0:n2-1)
  347. * T1 -> a(n1+1,0), T2 -> a(n1,0), S -> a(0,0)
  348. * T1 -> a(n2), T2 -> a(n1), S -> a(0)
  349. *
  350. IJP = 0
  351. DO J = 0, N1 - 1
  352. IJ = N2 + J
  353. DO I = 0, J
  354. AP( IJP ) = CONJG( ARF( IJ ) )
  355. IJP = IJP + 1
  356. IJ = IJ + LDA
  357. END DO
  358. END DO
  359. JS = 0
  360. DO J = N1, N - 1
  361. IJ = JS
  362. DO IJ = JS, JS + J
  363. AP( IJP ) = ARF( IJ )
  364. IJP = IJP + 1
  365. END DO
  366. JS = JS + LDA
  367. END DO
  368. *
  369. END IF
  370. *
  371. ELSE
  372. *
  373. * N is odd and TRANSR = 'C'
  374. *
  375. IF( LOWER ) THEN
  376. *
  377. * SRPA for LOWER, TRANSPOSE and N is odd
  378. * T1 -> A(0,0) , T2 -> A(1,0) , S -> A(0,n1)
  379. * T1 -> a(0+0) , T2 -> a(1+0) , S -> a(0+n1*n1); lda=n1
  380. *
  381. IJP = 0
  382. DO I = 0, N2
  383. DO IJ = I*( LDA+1 ), N*LDA - 1, LDA
  384. AP( IJP ) = CONJG( ARF( IJ ) )
  385. IJP = IJP + 1
  386. END DO
  387. END DO
  388. JS = 1
  389. DO J = 0, N2 - 1
  390. DO IJ = JS, JS + N2 - J - 1
  391. AP( IJP ) = ARF( IJ )
  392. IJP = IJP + 1
  393. END DO
  394. JS = JS + LDA + 1
  395. END DO
  396. *
  397. ELSE
  398. *
  399. * SRPA for UPPER, TRANSPOSE and N is odd
  400. * T1 -> A(0,n1+1), T2 -> A(0,n1), S -> A(0,0)
  401. * T1 -> a(n2*n2), T2 -> a(n1*n2), S -> a(0); lda = n2
  402. *
  403. IJP = 0
  404. JS = N2*LDA
  405. DO J = 0, N1 - 1
  406. DO IJ = JS, JS + J
  407. AP( IJP ) = ARF( IJ )
  408. IJP = IJP + 1
  409. END DO
  410. JS = JS + LDA
  411. END DO
  412. DO I = 0, N1
  413. DO IJ = I, I + ( N1+I )*LDA, LDA
  414. AP( IJP ) = CONJG( ARF( IJ ) )
  415. IJP = IJP + 1
  416. END DO
  417. END DO
  418. *
  419. END IF
  420. *
  421. END IF
  422. *
  423. ELSE
  424. *
  425. * N is even
  426. *
  427. IF( NORMALTRANSR ) THEN
  428. *
  429. * N is even and TRANSR = 'N'
  430. *
  431. IF( LOWER ) THEN
  432. *
  433. * SRPA for LOWER, NORMAL, and N is even ( a(0:n,0:k-1) )
  434. * T1 -> a(1,0), T2 -> a(0,0), S -> a(k+1,0)
  435. * T1 -> a(1), T2 -> a(0), S -> a(k+1)
  436. *
  437. IJP = 0
  438. JP = 0
  439. DO J = 0, K - 1
  440. DO I = J, N - 1
  441. IJ = 1 + I + JP
  442. AP( IJP ) = ARF( IJ )
  443. IJP = IJP + 1
  444. END DO
  445. JP = JP + LDA
  446. END DO
  447. DO I = 0, K - 1
  448. DO J = I, K - 1
  449. IJ = I + J*LDA
  450. AP( IJP ) = CONJG( ARF( IJ ) )
  451. IJP = IJP + 1
  452. END DO
  453. END DO
  454. *
  455. ELSE
  456. *
  457. * SRPA for UPPER, NORMAL, and N is even ( a(0:n,0:k-1) )
  458. * T1 -> a(k+1,0) , T2 -> a(k,0), S -> a(0,0)
  459. * T1 -> a(k+1), T2 -> a(k), S -> a(0)
  460. *
  461. IJP = 0
  462. DO J = 0, K - 1
  463. IJ = K + 1 + J
  464. DO I = 0, J
  465. AP( IJP ) = CONJG( ARF( IJ ) )
  466. IJP = IJP + 1
  467. IJ = IJ + LDA
  468. END DO
  469. END DO
  470. JS = 0
  471. DO J = K, N - 1
  472. IJ = JS
  473. DO IJ = JS, JS + J
  474. AP( IJP ) = ARF( IJ )
  475. IJP = IJP + 1
  476. END DO
  477. JS = JS + LDA
  478. END DO
  479. *
  480. END IF
  481. *
  482. ELSE
  483. *
  484. * N is even and TRANSR = 'C'
  485. *
  486. IF( LOWER ) THEN
  487. *
  488. * SRPA for LOWER, TRANSPOSE and N is even (see paper)
  489. * T1 -> B(0,1), T2 -> B(0,0), S -> B(0,k+1)
  490. * T1 -> a(0+k), T2 -> a(0+0), S -> a(0+k*(k+1)); lda=k
  491. *
  492. IJP = 0
  493. DO I = 0, K - 1
  494. DO IJ = I + ( I+1 )*LDA, ( N+1 )*LDA - 1, LDA
  495. AP( IJP ) = CONJG( ARF( IJ ) )
  496. IJP = IJP + 1
  497. END DO
  498. END DO
  499. JS = 0
  500. DO J = 0, K - 1
  501. DO IJ = JS, JS + K - J - 1
  502. AP( IJP ) = ARF( IJ )
  503. IJP = IJP + 1
  504. END DO
  505. JS = JS + LDA + 1
  506. END DO
  507. *
  508. ELSE
  509. *
  510. * SRPA for UPPER, TRANSPOSE and N is even (see paper)
  511. * T1 -> B(0,k+1), T2 -> B(0,k), S -> B(0,0)
  512. * T1 -> a(0+k*(k+1)), T2 -> a(0+k*k), S -> a(0+0)); lda=k
  513. *
  514. IJP = 0
  515. JS = ( K+1 )*LDA
  516. DO J = 0, K - 1
  517. DO IJ = JS, JS + J
  518. AP( IJP ) = ARF( IJ )
  519. IJP = IJP + 1
  520. END DO
  521. JS = JS + LDA
  522. END DO
  523. DO I = 0, K - 1
  524. DO IJ = I, I + ( K+I )*LDA, LDA
  525. AP( IJP ) = CONJG( ARF( IJ ) )
  526. IJP = IJP + 1
  527. END DO
  528. END DO
  529. *
  530. END IF
  531. *
  532. END IF
  533. *
  534. END IF
  535. *
  536. RETURN
  537. *
  538. * End of CTFTTP
  539. *
  540. END