LAPACK 3.3.1 Linear Algebra PACKage

strttf.f

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