LAPACK 3.3.1 Linear Algebra PACKage

# stfttr.f

Go to the documentation of this file.
```00001       SUBROUTINE STFTTR( TRANSR, UPLO, N, ARF, A, LDA, INFO )
00002 *
00003 *  -- LAPACK routine (version 3.3.1)                                    --
00004 *
00005 *  -- Contributed by Fred Gustavson of the IBM Watson Research Center --
00006 *  -- April 2011                                                      --
00007 *
00008 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00009 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00010 *
00011 *     .. Scalar Arguments ..
00012       CHARACTER          TRANSR, UPLO
00013       INTEGER            INFO, N, LDA
00014 *     ..
00015 *     .. Array Arguments ..
00016       REAL               A( 0: LDA-1, 0: * ), ARF( 0: * )
00017 *     ..
00018 *
00019 *  Purpose
00020 *  =======
00021 *
00022 *  STFTTR copies a triangular matrix A from rectangular full packed
00023 *  format (TF) to standard full format (TR).
00024 *
00025 *  Arguments
00026 *  =========
00027 *
00028 *  TRANSR   (input) CHARACTER*1
00029 *          = 'N':  ARF is in Normal format;
00030 *          = 'T':  ARF is in Transpose format.
00031 *
00032 *  UPLO    (input) CHARACTER*1
00033 *          = 'U':  A is upper triangular;
00034 *          = 'L':  A is lower triangular.
00035 *
00036 *  N       (input) INTEGER
00037 *          The order of the matrices ARF and A. N >= 0.
00038 *
00039 *  ARF     (input) REAL array, dimension (N*(N+1)/2).
00040 *          On entry, the upper (if UPLO = 'U') or lower (if UPLO = 'L')
00041 *          matrix A in RFP format. See the "Notes" below for more
00042 *          details.
00043 *
00044 *  A       (output) REAL array, dimension (LDA,N)
00045 *          On exit, the triangular matrix A.  If UPLO = 'U', the
00046 *          leading N-by-N upper triangular part of the array A contains
00047 *          the upper triangular matrix, and the strictly lower
00048 *          triangular part of A is not referenced.  If UPLO = 'L', the
00049 *          leading N-by-N lower triangular part of the array A contains
00050 *          the lower triangular matrix, and the strictly upper
00051 *          triangular part of A is not referenced.
00052 *
00053 *  LDA     (input) INTEGER
00054 *          The leading dimension of the array A.  LDA >= max(1,N).
00055 *
00056 *  INFO    (output) INTEGER
00057 *          = 0:  successful exit
00058 *          < 0:  if INFO = -i, the i-th argument had an illegal value
00059 *
00060 *  Further Details
00061 *  ===============
00062 *
00063 *  We first consider Rectangular Full Packed (RFP) Format when N is
00064 *  even. We give an example where N = 6.
00065 *
00066 *      AP is Upper             AP is Lower
00067 *
00068 *   00 01 02 03 04 05       00
00069 *      11 12 13 14 15       10 11
00070 *         22 23 24 25       20 21 22
00071 *            33 34 35       30 31 32 33
00072 *               44 45       40 41 42 43 44
00073 *                  55       50 51 52 53 54 55
00074 *
00075 *
00076 *  Let TRANSR = 'N'. RFP holds AP as follows:
00077 *  For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last
00078 *  three columns of AP upper. The lower triangle A(4:6,0:2) consists of
00079 *  the transpose of the first three columns of AP upper.
00080 *  For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first
00081 *  three columns of AP lower. The upper triangle A(0:2,0:2) consists of
00082 *  the transpose of the last three columns of AP lower.
00083 *  This covers the case N even and TRANSR = 'N'.
00084 *
00085 *         RFP A                   RFP A
00086 *
00087 *        03 04 05                33 43 53
00088 *        13 14 15                00 44 54
00089 *        23 24 25                10 11 55
00090 *        33 34 35                20 21 22
00091 *        00 44 45                30 31 32
00092 *        01 11 55                40 41 42
00093 *        02 12 22                50 51 52
00094 *
00095 *  Now let TRANSR = 'T'. RFP A in both UPLO cases is just the
00096 *  transpose of RFP A above. One therefore gets:
00097 *
00098 *
00099 *           RFP A                   RFP A
00100 *
00101 *     03 13 23 33 00 01 02    33 00 10 20 30 40 50
00102 *     04 14 24 34 44 11 12    43 44 11 21 31 41 51
00103 *     05 15 25 35 45 55 22    53 54 55 22 32 42 52
00104 *
00105 *
00106 *  We then consider Rectangular Full Packed (RFP) Format when N is
00107 *  odd. We give an example where N = 5.
00108 *
00109 *     AP is Upper                 AP is Lower
00110 *
00111 *   00 01 02 03 04              00
00112 *      11 12 13 14              10 11
00113 *         22 23 24              20 21 22
00114 *            33 34              30 31 32 33
00115 *               44              40 41 42 43 44
00116 *
00117 *
00118 *  Let TRANSR = 'N'. RFP holds AP as follows:
00119 *  For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last
00120 *  three columns of AP upper. The lower triangle A(3:4,0:1) consists of
00121 *  the transpose of the first two columns of AP upper.
00122 *  For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first
00123 *  three columns of AP lower. The upper triangle A(0:1,1:2) consists of
00124 *  the transpose of the last two columns of AP lower.
00125 *  This covers the case N odd and TRANSR = 'N'.
00126 *
00127 *         RFP A                   RFP A
00128 *
00129 *        02 03 04                00 33 43
00130 *        12 13 14                10 11 44
00131 *        22 23 24                20 21 22
00132 *        00 33 34                30 31 32
00133 *        01 11 44                40 41 42
00134 *
00135 *  Now let TRANSR = 'T'. RFP A in both UPLO cases is just the
00136 *  transpose of RFP A above. One therefore gets:
00137 *
00138 *           RFP A                   RFP A
00139 *
00140 *     02 12 22 00 01             00 10 20 30 40 50
00141 *     03 13 23 33 11             33 11 21 31 41 51
00142 *     04 14 24 34 44             43 44 22 32 42 52
00143 *
00144 *  Reference
00145 *  =========
00146 *
00147 *  =====================================================================
00148 *
00149 *     ..
00150 *     .. Local Scalars ..
00151       LOGICAL            LOWER, NISODD, NORMALTRANSR
00152       INTEGER            N1, N2, K, NT, NX2, NP1X2
00153       INTEGER            I, J, L, IJ
00154 *     ..
00155 *     .. External Functions ..
00156       LOGICAL            LSAME
00157       EXTERNAL           LSAME
00158 *     ..
00159 *     .. External Subroutines ..
00160       EXTERNAL           XERBLA
00161 *     ..
00162 *     .. Intrinsic Functions ..
00163       INTRINSIC          MAX, MOD
00164 *     ..
00165 *     .. Executable Statements ..
00166 *
00167 *     Test the input parameters.
00168 *
00169       INFO = 0
00170       NORMALTRANSR = LSAME( TRANSR, 'N' )
00171       LOWER = LSAME( UPLO, 'L' )
00172       IF( .NOT.NORMALTRANSR .AND. .NOT.LSAME( TRANSR, 'T' ) ) THEN
00173          INFO = -1
00174       ELSE IF( .NOT.LOWER .AND. .NOT.LSAME( UPLO, 'U' ) ) THEN
00175          INFO = -2
00176       ELSE IF( N.LT.0 ) THEN
00177          INFO = -3
00178       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
00179          INFO = -6
00180       END IF
00181       IF( INFO.NE.0 ) THEN
00182          CALL XERBLA( 'STFTTR', -INFO )
00183          RETURN
00184       END IF
00185 *
00186 *     Quick return if possible
00187 *
00188       IF( N.LE.1 ) THEN
00189          IF( N.EQ.1 ) THEN
00190             A( 0, 0 ) = ARF( 0 )
00191          END IF
00192          RETURN
00193       END IF
00194 *
00195 *     Size of array ARF(0:nt-1)
00196 *
00197       NT = N*( N+1 ) / 2
00198 *
00199 *     set N1 and N2 depending on LOWER: for N even N1=N2=K
00200 *
00201       IF( LOWER ) THEN
00202          N2 = N / 2
00203          N1 = N - N2
00204       ELSE
00205          N1 = N / 2
00206          N2 = N - N1
00207       END IF
00208 *
00209 *     If N is odd, set NISODD = .TRUE., LDA=N+1 and A is (N+1)--by--K2.
00210 *     If N is even, set K = N/2 and NISODD = .FALSE., LDA=N and A is
00211 *     N--by--(N+1)/2.
00212 *
00213       IF( MOD( N, 2 ).EQ.0 ) THEN
00214          K = N / 2
00215          NISODD = .FALSE.
00216          IF( .NOT.LOWER )
00217      \$      NP1X2 = N + N + 2
00218       ELSE
00219          NISODD = .TRUE.
00220          IF( .NOT.LOWER )
00221      \$      NX2 = N + N
00222       END IF
00223 *
00224       IF( NISODD ) THEN
00225 *
00226 *        N is odd
00227 *
00228          IF( NORMALTRANSR ) THEN
00229 *
00230 *           N is odd and TRANSR = 'N'
00231 *
00232             IF( LOWER ) THEN
00233 *
00234 *              N is odd, TRANSR = 'N', and UPLO = 'L'
00235 *
00236                IJ = 0
00237                DO J = 0, N2
00238                   DO I = N1, N2 + J
00239                      A( N2+J, I ) = ARF( IJ )
00240                      IJ = IJ + 1
00241                   END DO
00242                   DO I = J, N - 1
00243                      A( I, J ) = ARF( IJ )
00244                      IJ = IJ + 1
00245                   END DO
00246                END DO
00247 *
00248             ELSE
00249 *
00250 *              N is odd, TRANSR = 'N', and UPLO = 'U'
00251 *
00252                IJ = NT - N
00253                DO J = N - 1, N1, -1
00254                   DO I = 0, J
00255                      A( I, J ) = ARF( IJ )
00256                      IJ = IJ + 1
00257                   END DO
00258                   DO L = J - N1, N1 - 1
00259                      A( J-N1, L ) = ARF( IJ )
00260                      IJ = IJ + 1
00261                   END DO
00262                   IJ = IJ - NX2
00263                END DO
00264 *
00265             END IF
00266 *
00267          ELSE
00268 *
00269 *           N is odd and TRANSR = 'T'
00270 *
00271             IF( LOWER ) THEN
00272 *
00273 *              N is odd, TRANSR = 'T', and UPLO = 'L'
00274 *
00275                IJ = 0
00276                DO J = 0, N2 - 1
00277                   DO I = 0, J
00278                      A( J, I ) = ARF( IJ )
00279                      IJ = IJ + 1
00280                   END DO
00281                   DO I = N1 + J, N - 1
00282                      A( I, N1+J ) = ARF( IJ )
00283                      IJ = IJ + 1
00284                   END DO
00285                END DO
00286                DO J = N2, N - 1
00287                   DO I = 0, N1 - 1
00288                      A( J, I ) = ARF( IJ )
00289                      IJ = IJ + 1
00290                   END DO
00291                END DO
00292 *
00293             ELSE
00294 *
00295 *              N is odd, TRANSR = 'T', and UPLO = 'U'
00296 *
00297                IJ = 0
00298                DO J = 0, N1
00299                   DO I = N1, N - 1
00300                      A( J, I ) = ARF( IJ )
00301                      IJ = IJ + 1
00302                   END DO
00303                END DO
00304                DO J = 0, N1 - 1
00305                   DO I = 0, J
00306                      A( I, J ) = ARF( IJ )
00307                      IJ = IJ + 1
00308                   END DO
00309                   DO L = N2 + J, N - 1
00310                      A( N2+J, L ) = ARF( IJ )
00311                      IJ = IJ + 1
00312                   END DO
00313                END DO
00314 *
00315             END IF
00316 *
00317          END IF
00318 *
00319       ELSE
00320 *
00321 *        N is even
00322 *
00323          IF( NORMALTRANSR ) THEN
00324 *
00325 *           N is even and TRANSR = 'N'
00326 *
00327             IF( LOWER ) THEN
00328 *
00329 *              N is even, TRANSR = 'N', and UPLO = 'L'
00330 *
00331                IJ = 0
00332                DO J = 0, K - 1
00333                   DO I = K, K + J
00334                      A( K+J, I ) = ARF( IJ )
00335                      IJ = IJ + 1
00336                   END DO
00337                   DO I = J, N - 1
00338                      A( I, J ) = ARF( IJ )
00339                      IJ = IJ + 1
00340                   END DO
00341                END DO
00342 *
00343             ELSE
00344 *
00345 *              N is even, TRANSR = 'N', and UPLO = 'U'
00346 *
00347                IJ = NT - N - 1
00348                DO J = N - 1, K, -1
00349                   DO I = 0, J
00350                      A( I, J ) = ARF( IJ )
00351                      IJ = IJ + 1
00352                   END DO
00353                   DO L = J - K, K - 1
00354                      A( J-K, L ) = ARF( IJ )
00355                      IJ = IJ + 1
00356                   END DO
00357                   IJ = IJ - NP1X2
00358                END DO
00359 *
00360             END IF
00361 *
00362          ELSE
00363 *
00364 *           N is even and TRANSR = 'T'
00365 *
00366             IF( LOWER ) THEN
00367 *
00368 *              N is even, TRANSR = 'T', and UPLO = 'L'
00369 *
00370                IJ = 0
00371                J = K
00372                DO I = K, N - 1
00373                   A( I, J ) = ARF( IJ )
00374                   IJ = IJ + 1
00375                END DO
00376                DO J = 0, K - 2
00377                   DO I = 0, J
00378                      A( J, I ) = ARF( IJ )
00379                      IJ = IJ + 1
00380                   END DO
00381                   DO I = K + 1 + J, N - 1
00382                      A( I, K+1+J ) = ARF( IJ )
00383                      IJ = IJ + 1
00384                   END DO
00385                END DO
00386                DO J = K - 1, N - 1
00387                   DO I = 0, K - 1
00388                      A( J, I ) = ARF( IJ )
00389                      IJ = IJ + 1
00390                   END DO
00391                END DO
00392 *
00393             ELSE
00394 *
00395 *              N is even, TRANSR = 'T', and UPLO = 'U'
00396 *
00397                IJ = 0
00398                DO J = 0, K
00399                   DO I = K, N - 1
00400                      A( J, I ) = ARF( IJ )
00401                      IJ = IJ + 1
00402                   END DO
00403                END DO
00404                DO J = 0, K - 2
00405                   DO I = 0, J
00406                      A( I, J ) = ARF( IJ )
00407                      IJ = IJ + 1
00408                   END DO
00409                   DO L = K + 1 + J, N - 1
00410                      A( K+1+J, L ) = ARF( IJ )
00411                      IJ = IJ + 1
00412                   END DO
00413                END DO
00414 *              Note that here, on exit of the loop, J = K-1
00415                DO I = 0, J
00416                   A( I, J ) = ARF( IJ )
00417                   IJ = IJ + 1
00418                END DO
00419 *
00420             END IF
00421 *
00422          END IF
00423 *
00424       END IF
00425 *
00426       RETURN
00427 *
00428 *     End of STFTTR
00429 *
00430       END
```