1       SUBROUTINE DTPTTF( TRANSR, UPLO, N, AP, ARF, INFO )
  2 *
  3 *  -- LAPACK routine (version 3.3.1)                                    --
  4 *
  5 *  -- Contributed by Fred Gustavson of the IBM Watson Research Center --
  6 *  -- April 2011                                                      --
  7 *
  8 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  9 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 10 *
 11 *     ..
 12 *     .. Scalar Arguments ..
 13       CHARACTER          TRANSR, UPLO
 14       INTEGER            INFO, N
 15 *     ..
 16 *     .. Array Arguments ..
 17       DOUBLE PRECISION   AP( 0* ), ARF( 0* )
 18 *
 19 *  Purpose
 20 *  =======
 21 *
 22 *  DTPTTF copies a triangular matrix A from standard packed format (TP)
 23 *  to rectangular full packed format (TF).
 24 *
 25 *  Arguments
 26 *  =========
 27 *
 28 *  TRANSR  (input) CHARACTER*1
 29 *          = 'N':  ARF in Normal format is wanted;
 30 *          = 'T':  ARF in Conjugate-transpose format is wanted.
 31 *
 32 *  UPLO    (input) CHARACTER*1
 33 *          = 'U':  A is upper triangular;
 34 *          = 'L':  A is lower triangular.
 35 *
 36 *  N       (input) INTEGER
 37 *          The order of the matrix A.  N >= 0.
 38 *
 39 *  AP      (input) DOUBLE PRECISION array, dimension ( N*(N+1)/2 ),
 40 *          On entry, the upper or lower triangular matrix A, packed
 41 *          columnwise in a linear array. The j-th column of A is stored
 42 *          in the array AP as follows:
 43 *          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
 44 *          if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.
 45 *
 46 *  ARF     (output) DOUBLE PRECISION array, dimension ( N*(N+1)/2 ),
 47 *          On exit, the upper or lower triangular matrix A stored in
 48 *          RFP format. For a further discussion see Notes below.
 49 *
 50 *  INFO    (output) INTEGER
 51 *          = 0:  successful exit
 52 *          < 0:  if INFO = -i, the i-th argument had an illegal value
 53 *
 54 *  Further Details
 55 *  ===============
 56 *
 57 *  We first consider Rectangular Full Packed (RFP) Format when N is
 58 *  even. We give an example where N = 6.
 59 *
 60 *      AP is Upper             AP is Lower
 61 *
 62 *   00 01 02 03 04 05       00
 63 *      11 12 13 14 15       10 11
 64 *         22 23 24 25       20 21 22
 65 *            33 34 35       30 31 32 33
 66 *               44 45       40 41 42 43 44
 67 *                  55       50 51 52 53 54 55
 68 *
 69 *
 70 *  Let TRANSR = 'N'. RFP holds AP as follows:
 71 *  For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last
 72 *  three columns of AP upper. The lower triangle A(4:6,0:2) consists of
 73 *  the transpose of the first three columns of AP upper.
 74 *  For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first
 75 *  three columns of AP lower. The upper triangle A(0:2,0:2) consists of
 76 *  the transpose of the last three columns of AP lower.
 77 *  This covers the case N even and TRANSR = 'N'.
 78 *
 79 *         RFP A                   RFP A
 80 *
 81 *        03 04 05                33 43 53
 82 *        13 14 15                00 44 54
 83 *        23 24 25                10 11 55
 84 *        33 34 35                20 21 22
 85 *        00 44 45                30 31 32
 86 *        01 11 55                40 41 42
 87 *        02 12 22                50 51 52
 88 *
 89 *  Now let TRANSR = 'T'. RFP A in both UPLO cases is just the
 90 *  transpose of RFP A above. One therefore gets:
 91 *
 92 *
 93 *           RFP A                   RFP A
 94 *
 95 *     03 13 23 33 00 01 02    33 00 10 20 30 40 50
 96 *     04 14 24 34 44 11 12    43 44 11 21 31 41 51
 97 *     05 15 25 35 45 55 22    53 54 55 22 32 42 52
 98 *
 99 *
100 *  We then consider Rectangular Full Packed (RFP) Format when N is
101 *  odd. We give an example where N = 5.
102 *
103 *     AP is Upper                 AP is Lower
104 *
105 *   00 01 02 03 04              00
106 *      11 12 13 14              10 11
107 *         22 23 24              20 21 22
108 *            33 34              30 31 32 33
109 *               44              40 41 42 43 44
110 *
111 *
112 *  Let TRANSR = 'N'. RFP holds AP as follows:
113 *  For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last
114 *  three columns of AP upper. The lower triangle A(3:4,0:1) consists of
115 *  the transpose of the first two columns of AP upper.
116 *  For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first
117 *  three columns of AP lower. The upper triangle A(0:1,1:2) consists of
118 *  the transpose of the last two columns of AP lower.
119 *  This covers the case N odd and TRANSR = 'N'.
120 *
121 *         RFP A                   RFP A
122 *
123 *        02 03 04                00 33 43
124 *        12 13 14                10 11 44
125 *        22 23 24                20 21 22
126 *        00 33 34                30 31 32
127 *        01 11 44                40 41 42
128 *
129 *  Now let TRANSR = 'T'. RFP A in both UPLO cases is just the
130 *  transpose of RFP A above. One therefore gets:
131 *
132 *           RFP A                   RFP A
133 *
134 *     02 12 22 00 01             00 10 20 30 40 50
135 *     03 13 23 33 11             33 11 21 31 41 51
136 *     04 14 24 34 44             43 44 22 32 42 52
137 *
138 *  =====================================================================
139 *
140 *     .. Parameters ..
141 *     ..
142 *     .. Local Scalars ..
143       LOGICAL            LOWER, NISODD, NORMALTRANSR
144       INTEGER            N1, N2, K, NT
145       INTEGER            I, J, IJ
146       INTEGER            IJP, JP, LDA, JS
147 *     ..
148 *     .. External Functions ..
149       LOGICAL            LSAME
150       EXTERNAL           LSAME
151 *     ..
152 *     .. External Subroutines ..
153       EXTERNAL           XERBLA
154 *     ..
155 *     .. Intrinsic Functions ..
156       INTRINSIC          MOD
157 *     ..
158 *     .. Executable Statements ..
159 *
160 *     Test the input parameters.
161 *
162       INFO = 0
163       NORMALTRANSR = LSAME( TRANSR, 'N' )
164       LOWER = LSAME( UPLO, 'L' )
165       IF.NOT.NORMALTRANSR .AND. .NOT.LSAME( TRANSR, 'T' ) ) THEN
166          INFO = -1
167       ELSE IF.NOT.LOWER .AND. .NOT.LSAME( UPLO, 'U' ) ) THEN
168          INFO = -2
169       ELSE IF( N.LT.0 ) THEN
170          INFO = -3
171       END IF
172       IF( INFO.NE.0 ) THEN
173          CALL XERBLA( 'DTPTTF'-INFO )
174          RETURN
175       END IF
176 *
177 *     Quick return if possible
178 *
179       IF( N.EQ.0 )
180      $   RETURN
181 *
182       IF( N.EQ.1 ) THEN
183          IF( NORMALTRANSR ) THEN
184             ARF( 0 ) = AP( 0 )
185          ELSE
186             ARF( 0 ) = AP( 0 )
187          END IF
188          RETURN
189       END IF
190 *
191 *     Size of array ARF(0:NT-1)
192 *
193       NT = N*( N+1 ) / 2
194 *
195 *     Set N1 and N2 depending on LOWER
196 *
197       IF( LOWER ) THEN
198          N2 = N / 2
199          N1 = N - N2
200       ELSE
201          N1 = N / 2
202          N2 = N - N1
203       END IF
204 *
205 *     If N is odd, set NISODD = .TRUE.
206 *     If N is even, set K = N/2 and NISODD = .FALSE.
207 *
208 *     set lda of ARF^C; ARF^C is (0:(N+1)/2-1,0:N-noe)
209 *     where noe = 0 if n is even, noe = 1 if n is odd
210 *
211       IFMOD( N, 2 ).EQ.0 ) THEN
212          K = N / 2
213          NISODD = .FALSE.
214          LDA = N + 1
215       ELSE
216          NISODD = .TRUE.
217          LDA = N
218       END IF
219 *
220 *     ARF^C has lda rows and n+1-noe cols
221 *
222       IF.NOT.NORMALTRANSR )
223      $   LDA = ( N+1 ) / 2
224 *
225 *     start execution: there are eight cases
226 *
227       IF( NISODD ) THEN
228 *
229 *        N is odd
230 *
231          IF( NORMALTRANSR ) THEN
232 *
233 *           N is odd and TRANSR = 'N'
234 *
235             IF( LOWER ) THEN
236 *
237 *              N is odd, TRANSR = 'N', and UPLO = 'L'
238 *
239                IJP = 0
240                JP = 0
241                DO J = 0, N2
242                   DO I = J, N - 1
243                      IJ = I + JP
244                      ARF( IJ ) = AP( IJP )
245                      IJP = IJP + 1
246                   END DO
247                   JP = JP + LDA
248                END DO
249                DO I = 0, N2 - 1
250                   DO J = 1 + I, N2
251                      IJ = I + J*LDA
252                      ARF( IJ ) = AP( IJP )
253                      IJP = IJP + 1
254                   END DO
255                END DO
256 *
257             ELSE
258 *
259 *              N is odd, TRANSR = 'N', and UPLO = 'U'
260 *
261                IJP = 0
262                DO J = 0, N1 - 1
263                   IJ = N2 + J
264                   DO I = 0, J
265                      ARF( IJ ) = AP( IJP )
266                      IJP = IJP + 1
267                      IJ = IJ + LDA
268                   END DO
269                END DO
270                JS = 0
271                DO J = N1, N - 1
272                   IJ = JS
273                   DO IJ = JS, JS + J
274                      ARF( IJ ) = AP( IJP )
275                      IJP = IJP + 1
276                   END DO
277                   JS = JS + LDA
278                END DO
279 *
280             END IF
281 *
282          ELSE
283 *
284 *           N is odd and TRANSR = 'T'
285 *
286             IF( LOWER ) THEN
287 *
288 *              N is odd, TRANSR = 'T', and UPLO = 'L'
289 *
290                IJP = 0
291                DO I = 0, N2
292                   DO IJ = I*( LDA+1 ), N*LDA - 1, LDA
293                      ARF( IJ ) = AP( IJP )
294                      IJP = IJP + 1
295                   END DO
296                END DO
297                JS = 1
298                DO J = 0, N2 - 1
299                   DO IJ = JS, JS + N2 - J - 1
300                      ARF( IJ ) = AP( IJP )
301                      IJP = IJP + 1
302                   END DO
303                   JS = JS + LDA + 1
304                END DO
305 *
306             ELSE
307 *
308 *              N is odd, TRANSR = 'T', and UPLO = 'U'
309 *
310                IJP = 0
311                JS = N2*LDA
312                DO J = 0, N1 - 1
313                   DO IJ = JS, JS + J
314                      ARF( IJ ) = AP( IJP )
315                      IJP = IJP + 1
316                   END DO
317                   JS = JS + LDA
318                END DO
319                DO I = 0, N1
320                   DO IJ = I, I + ( N1+I )*LDA, LDA
321                      ARF( IJ ) = AP( IJP )
322                      IJP = IJP + 1
323                   END DO
324                END DO
325 *
326             END IF
327 *
328          END IF
329 *
330       ELSE
331 *
332 *        N is even
333 *
334          IF( NORMALTRANSR ) THEN
335 *
336 *           N is even and TRANSR = 'N'
337 *
338             IF( LOWER ) THEN
339 *
340 *              N is even, TRANSR = 'N', and UPLO = 'L'
341 *
342                IJP = 0
343                JP = 0
344                DO J = 0, K - 1
345                   DO I = J, N - 1
346                      IJ = 1 + I + JP
347                      ARF( IJ ) = AP( IJP )
348                      IJP = IJP + 1
349                   END DO
350                   JP = JP + LDA
351                END DO
352                DO I = 0, K - 1
353                   DO J = I, K - 1
354                      IJ = I + J*LDA
355                      ARF( IJ ) = AP( IJP )
356                      IJP = IJP + 1
357                   END DO
358                END DO
359 *
360             ELSE
361 *
362 *              N is even, TRANSR = 'N', and UPLO = 'U'
363 *
364                IJP = 0
365                DO J = 0, K - 1
366                   IJ = K + 1 + J
367                   DO I = 0, J
368                      ARF( IJ ) = AP( IJP )
369                      IJP = IJP + 1
370                      IJ = IJ + LDA
371                   END DO
372                END DO
373                JS = 0
374                DO J = K, N - 1
375                   IJ = JS
376                   DO IJ = JS, JS + J
377                      ARF( IJ ) = AP( IJP )
378                      IJP = IJP + 1
379                   END DO
380                   JS = JS + LDA
381                END DO
382 *
383             END IF
384 *
385          ELSE
386 *
387 *           N is even and TRANSR = 'T'
388 *
389             IF( LOWER ) THEN
390 *
391 *              N is even, TRANSR = 'T', and UPLO = 'L'
392 *
393                IJP = 0
394                DO I = 0, K - 1
395                   DO IJ = I + ( I+1 )*LDA, ( N+1 )*LDA - 1, LDA
396                      ARF( IJ ) = AP( IJP )
397                      IJP = IJP + 1
398                   END DO
399                END DO
400                JS = 0
401                DO J = 0, K - 1
402                   DO IJ = JS, JS + K - J - 1
403                      ARF( IJ ) = AP( IJP )
404                      IJP = IJP + 1
405                   END DO
406                   JS = JS + LDA + 1
407                END DO
408 *
409             ELSE
410 *
411 *              N is even, TRANSR = 'T', and UPLO = 'U'
412 *
413                IJP = 0
414                JS = ( K+1 )*LDA
415                DO J = 0, K - 1
416                   DO IJ = JS, JS + J
417                      ARF( IJ ) = AP( IJP )
418                      IJP = IJP + 1
419                   END DO
420                   JS = JS + LDA
421                END DO
422                DO I = 0, K - 1
423                   DO IJ = I, I + ( K+I )*LDA, LDA
424                      ARF( IJ ) = AP( IJP )
425                      IJP = IJP + 1
426                   END DO
427                END DO
428 *
429             END IF
430 *
431          END IF
432 *
433       END IF
434 *
435       RETURN
436 *
437 *     End of DTPTTF
438 *
439       END