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

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