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.

ztpttf.f 16 kB

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