LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine stfttp ( character  TRANSR,
character  UPLO,
integer  N,
real, dimension( 0: * )  ARF,
real, dimension( 0: * )  AP,
integer  INFO 
)

STFTTP copies a triangular matrix from the rectangular full packed format (TF) to the standard packed format (TP).

Download STFTTP + dependencies [TGZ] [ZIP] [TXT]

Purpose: 
 STFTTP copies a triangular matrix A from rectangular full packed
 format (TF) to standard packed format (TP).
Parameters
[in]TRANSR
          TRANSR is CHARACTER*1
          = 'N':  ARF is in Normal format;
          = 'T':  ARF is in Transpose format;
[in]UPLO
          UPLO is CHARACTER*1
          = 'U':  A is upper triangular;
          = 'L':  A is lower triangular.
[in]N
          N is INTEGER
          The order of the matrix A. N >= 0.
[in]ARF
          ARF is REAL array, dimension ( N*(N+1)/2 ),
          On entry, the upper or lower triangular matrix A stored in
          RFP format. For a further discussion see Notes below.
[out]AP
          AP is REAL array, dimension ( N*(N+1)/2 ),
          On exit, the upper or lower triangular matrix A, packed
          columnwise in a linear array. The j-th column of A is stored
          in the array AP as follows:
          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
          if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.
[out]INFO
          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
September 2012
Further Details: 
  We first consider Rectangular Full Packed (RFP) Format when N is
  even. We give an example where N = 6.

      AP is Upper             AP is Lower

   00 01 02 03 04 05       00
      11 12 13 14 15       10 11
         22 23 24 25       20 21 22
            33 34 35       30 31 32 33
               44 45       40 41 42 43 44
                  55       50 51 52 53 54 55


  Let TRANSR = 'N'. RFP holds AP as follows:
  For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last
  three columns of AP upper. The lower triangle A(4:6,0:2) consists of
  the transpose of the first three columns of AP upper.
  For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first
  three columns of AP lower. The upper triangle A(0:2,0:2) consists of
  the transpose of the last three columns of AP lower.
  This covers the case N even and TRANSR = 'N'.

         RFP A                   RFP A

        03 04 05                33 43 53
        13 14 15                00 44 54
        23 24 25                10 11 55
        33 34 35                20 21 22
        00 44 45                30 31 32
        01 11 55                40 41 42
        02 12 22                50 51 52

  Now let TRANSR = 'T'. RFP A in both UPLO cases is just the
  transpose of RFP A above. One therefore gets:


           RFP A                   RFP A

     03 13 23 33 00 01 02    33 00 10 20 30 40 50
     04 14 24 34 44 11 12    43 44 11 21 31 41 51
     05 15 25 35 45 55 22    53 54 55 22 32 42 52


  We then consider Rectangular Full Packed (RFP) Format when N is
  odd. We give an example where N = 5.

     AP is Upper                 AP is Lower

   00 01 02 03 04              00
      11 12 13 14              10 11
         22 23 24              20 21 22
            33 34              30 31 32 33
               44              40 41 42 43 44


  Let TRANSR = 'N'. RFP holds AP as follows:
  For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last
  three columns of AP upper. The lower triangle A(3:4,0:1) consists of
  the transpose of the first two columns of AP upper.
  For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first
  three columns of AP lower. The upper triangle A(0:1,1:2) consists of
  the transpose of the last two columns of AP lower.
  This covers the case N odd and TRANSR = 'N'.

         RFP A                   RFP A

        02 03 04                00 33 43
        12 13 14                10 11 44
        22 23 24                20 21 22
        00 33 34                30 31 32
        01 11 44                40 41 42

  Now let TRANSR = 'T'. RFP A in both UPLO cases is just the
  transpose of RFP A above. One therefore gets:

           RFP A                   RFP A

     02 12 22 00 01             00 10 20 30 40 50
     03 13 23 33 11             33 11 21 31 41 51
     04 14 24 34 44             43 44 22 32 42 52

Definition at line 189 of file stfttp.f.

189 *
190 * -- LAPACK computational routine (version 3.4.2) --
191 * -- LAPACK is a software package provided by Univ. of Tennessee, --
192 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
193 * September 2012
194 *
195 * .. Scalar Arguments ..
196  CHARACTER transr, uplo
197  INTEGER info, n
198 * ..
199 * .. Array Arguments ..
200  REAL ap( 0: * ), arf( 0: * )
201 * ..
202 *
203 * =====================================================================
204 *
205 * .. Parameters ..
206 * ..
207 * .. Local Scalars ..
208  LOGICAL lower, nisodd, normaltransr
209  INTEGER n1, n2, k, nt
210  INTEGER i, j, ij
211  INTEGER ijp, jp, lda, js
212 * ..
213 * .. External Functions ..
214  LOGICAL lsame
215  EXTERNAL lsame
216 * ..
217 * .. External Subroutines ..
218  EXTERNAL xerbla
219 * ..
220 * .. Executable Statements ..
221 *
222 * Test the input parameters.
223 *
224  info = 0
225  normaltransr = lsame( transr, 'N' )
226  lower = lsame( uplo, 'L' )
227  IF( .NOT.normaltransr .AND. .NOT.lsame( transr, 'T' ) ) THEN
228  info = -1
229  ELSE IF( .NOT.lower .AND. .NOT.lsame( uplo, 'U' ) ) THEN
230  info = -2
231  ELSE IF( n.LT.0 ) THEN
232  info = -3
233  END IF
234  IF( info.NE.0 ) THEN
235  CALL xerbla( 'STFTTP', -info )
236  RETURN
237  END IF
238 *
239 * Quick return if possible
240 *
241  IF( n.EQ.0 )
242  $ RETURN
243 *
244  IF( n.EQ.1 ) THEN
245  IF( normaltransr ) THEN
246  ap( 0 ) = arf( 0 )
247  ELSE
248  ap( 0 ) = arf( 0 )
249  END IF
250  RETURN
251  END IF
252 *
253 * Size of array ARF(0:NT-1)
254 *
255  nt = n*( n+1 ) / 2
256 *
257 * Set N1 and N2 depending on LOWER
258 *
259  IF( lower ) THEN
260  n2 = n / 2
261  n1 = n - n2
262  ELSE
263  n1 = n / 2
264  n2 = n - n1
265  END IF
266 *
267 * If N is odd, set NISODD = .TRUE.
268 * If N is even, set K = N/2 and NISODD = .FALSE.
269 *
270 * set lda of ARF^C; ARF^C is (0:(N+1)/2-1,0:N-noe)
271 * where noe = 0 if n is even, noe = 1 if n is odd
272 *
273  IF( mod( n, 2 ).EQ.0 ) THEN
274  k = n / 2
275  nisodd = .false.
276  lda = n + 1
277  ELSE
278  nisodd = .true.
279  lda = n
280  END IF
281 *
282 * ARF^C has lda rows and n+1-noe cols
283 *
284  IF( .NOT.normaltransr )
285  $ lda = ( n+1 ) / 2
286 *
287 * start execution: there are eight cases
288 *
289  IF( nisodd ) THEN
290 *
291 * N is odd
292 *
293  IF( normaltransr ) THEN
294 *
295 * N is odd and TRANSR = 'N'
296 *
297  IF( lower ) THEN
298 *
299 * SRPA for LOWER, NORMAL and N is odd ( a(0:n-1,0:n1-1) )
300 * T1 -> a(0,0), T2 -> a(0,1), S -> a(n1,0)
301 * T1 -> a(0), T2 -> a(n), S -> a(n1); lda = n
302 *
303  ijp = 0
304  jp = 0
305  DO j = 0, n2
306  DO i = j, n - 1
307  ij = i + jp
308  ap( ijp ) = arf( ij )
309  ijp = ijp + 1
310  END DO
311  jp = jp + lda
312  END DO
313  DO i = 0, n2 - 1
314  DO j = 1 + i, n2
315  ij = i + j*lda
316  ap( ijp ) = arf( ij )
317  ijp = ijp + 1
318  END DO
319  END DO
320 *
321  ELSE
322 *
323 * SRPA for UPPER, NORMAL and N is odd ( a(0:n-1,0:n2-1)
324 * T1 -> a(n1+1,0), T2 -> a(n1,0), S -> a(0,0)
325 * T1 -> a(n2), T2 -> a(n1), S -> a(0)
326 *
327  ijp = 0
328  DO j = 0, n1 - 1
329  ij = n2 + j
330  DO i = 0, j
331  ap( ijp ) = arf( ij )
332  ijp = ijp + 1
333  ij = ij + lda
334  END DO
335  END DO
336  js = 0
337  DO j = n1, n - 1
338  ij = js
339  DO ij = js, js + j
340  ap( ijp ) = arf( ij )
341  ijp = ijp + 1
342  END DO
343  js = js + lda
344  END DO
345 *
346  END IF
347 *
348  ELSE
349 *
350 * N is odd and TRANSR = 'T'
351 *
352  IF( lower ) THEN
353 *
354 * SRPA for LOWER, TRANSPOSE and N is odd
355 * T1 -> A(0,0) , T2 -> A(1,0) , S -> A(0,n1)
356 * T1 -> a(0+0) , T2 -> a(1+0) , S -> a(0+n1*n1); lda=n1
357 *
358  ijp = 0
359  DO i = 0, n2
360  DO ij = i*( lda+1 ), n*lda - 1, lda
361  ap( ijp ) = arf( ij )
362  ijp = ijp + 1
363  END DO
364  END DO
365  js = 1
366  DO j = 0, n2 - 1
367  DO ij = js, js + n2 - j - 1
368  ap( ijp ) = arf( ij )
369  ijp = ijp + 1
370  END DO
371  js = js + lda + 1
372  END DO
373 *
374  ELSE
375 *
376 * SRPA for UPPER, TRANSPOSE and N is odd
377 * T1 -> A(0,n1+1), T2 -> A(0,n1), S -> A(0,0)
378 * T1 -> a(n2*n2), T2 -> a(n1*n2), S -> a(0); lda = n2
379 *
380  ijp = 0
381  js = n2*lda
382  DO j = 0, n1 - 1
383  DO ij = js, js + j
384  ap( ijp ) = arf( ij )
385  ijp = ijp + 1
386  END DO
387  js = js + lda
388  END DO
389  DO i = 0, n1
390  DO ij = i, i + ( n1+i )*lda, lda
391  ap( ijp ) = arf( ij )
392  ijp = ijp + 1
393  END DO
394  END DO
395 *
396  END IF
397 *
398  END IF
399 *
400  ELSE
401 *
402 * N is even
403 *
404  IF( normaltransr ) THEN
405 *
406 * N is even and TRANSR = 'N'
407 *
408  IF( lower ) THEN
409 *
410 * SRPA for LOWER, NORMAL, and N is even ( a(0:n,0:k-1) )
411 * T1 -> a(1,0), T2 -> a(0,0), S -> a(k+1,0)
412 * T1 -> a(1), T2 -> a(0), S -> a(k+1)
413 *
414  ijp = 0
415  jp = 0
416  DO j = 0, k - 1
417  DO i = j, n - 1
418  ij = 1 + i + jp
419  ap( ijp ) = arf( ij )
420  ijp = ijp + 1
421  END DO
422  jp = jp + lda
423  END DO
424  DO i = 0, k - 1
425  DO j = i, k - 1
426  ij = i + j*lda
427  ap( ijp ) = arf( ij )
428  ijp = ijp + 1
429  END DO
430  END DO
431 *
432  ELSE
433 *
434 * SRPA for UPPER, NORMAL, and N is even ( a(0:n,0:k-1) )
435 * T1 -> a(k+1,0) , T2 -> a(k,0), S -> a(0,0)
436 * T1 -> a(k+1), T2 -> a(k), S -> a(0)
437 *
438  ijp = 0
439  DO j = 0, k - 1
440  ij = k + 1 + j
441  DO i = 0, j
442  ap( ijp ) = arf( ij )
443  ijp = ijp + 1
444  ij = ij + lda
445  END DO
446  END DO
447  js = 0
448  DO j = k, n - 1
449  ij = js
450  DO ij = js, js + j
451  ap( ijp ) = arf( ij )
452  ijp = ijp + 1
453  END DO
454  js = js + lda
455  END DO
456 *
457  END IF
458 *
459  ELSE
460 *
461 * N is even and TRANSR = 'T'
462 *
463  IF( lower ) THEN
464 *
465 * SRPA for LOWER, TRANSPOSE and N is even (see paper)
466 * T1 -> B(0,1), T2 -> B(0,0), S -> B(0,k+1)
467 * T1 -> a(0+k), T2 -> a(0+0), S -> a(0+k*(k+1)); lda=k
468 *
469  ijp = 0
470  DO i = 0, k - 1
471  DO ij = i + ( i+1 )*lda, ( n+1 )*lda - 1, lda
472  ap( ijp ) = arf( ij )
473  ijp = ijp + 1
474  END DO
475  END DO
476  js = 0
477  DO j = 0, k - 1
478  DO ij = js, js + k - j - 1
479  ap( ijp ) = arf( ij )
480  ijp = ijp + 1
481  END DO
482  js = js + lda + 1
483  END DO
484 *
485  ELSE
486 *
487 * SRPA for UPPER, TRANSPOSE and N is even (see paper)
488 * T1 -> B(0,k+1), T2 -> B(0,k), S -> B(0,0)
489 * T1 -> a(0+k*(k+1)), T2 -> a(0+k*k), S -> a(0+0)); lda=k
490 *
491  ijp = 0
492  js = ( k+1 )*lda
493  DO j = 0, k - 1
494  DO ij = js, js + j
495  ap( ijp ) = arf( ij )
496  ijp = ijp + 1
497  END DO
498  js = js + lda
499  END DO
500  DO i = 0, k - 1
501  DO ij = i, i + ( k+i )*lda, lda
502  ap( ijp ) = arf( ij )
503  ijp = ijp + 1
504  END DO
505  END DO
506 *
507  END IF
508 *
509  END IF
510 *
511  END IF
512 *
513  RETURN
514 *
515 * End of STFTTP
516 *
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
lsame
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55

Here is the call graph for this function:

Here is the caller graph for this function: