LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ 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*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
Here is the call graph for this function:
Here is the caller graph for this function: