LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ stfttp()

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.
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 186 of file stfttp.f.

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