dtfttr.c
Go to the documentation of this file.
00001 /* dtfttr.f -- translated by f2c (version 20061008).
00002    You must link the resulting object file with libf2c:
00003         on Microsoft Windows system, link with libf2c.lib;
00004         on Linux or Unix systems, link with .../path/to/libf2c.a -lm
00005         or, if you install libf2c.a in a standard place, with -lf2c -lm
00006         -- in that order, at the end of the command line, as in
00007                 cc *.o -lf2c -lm
00008         Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
00009 
00010                 http://www.netlib.org/f2c/libf2c.zip
00011 */
00012 
00013 #include "f2c.h"
00014 #include "blaswrap.h"
00015 
00016 /* Subroutine */ int dtfttr_(char *transr, char *uplo, integer *n, doublereal 
00017         *arf, doublereal *a, integer *lda, integer *info)
00018 {
00019     /* System generated locals */
00020     integer a_dim1, a_offset, i__1, i__2;
00021 
00022     /* Local variables */
00023     integer i__, j, k, l, n1, n2, ij, nt, nx2, np1x2;
00024     logical normaltransr;
00025     extern logical lsame_(char *, char *);
00026     logical lower;
00027     extern /* Subroutine */ int xerbla_(char *, integer *);
00028     logical nisodd;
00029 
00030 
00031 /*  -- LAPACK routine (version 3.2)                                    -- */
00032 
00033 /*  -- Contributed by Fred Gustavson of the IBM Watson Research Center -- */
00034 /*  -- November 2008                                                   -- */
00035 
00036 /*  -- LAPACK is a software package provided by Univ. of Tennessee,    -- */
00037 /*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */
00038 
00039 /*     .. Scalar Arguments .. */
00040 /*     .. */
00041 /*     .. Array Arguments .. */
00042 /*     .. */
00043 
00044 /*  Purpose */
00045 /*  ======= */
00046 
00047 /*  DTFTTR copies a triangular matrix A from rectangular full packed */
00048 /*  format (TF) to standard full format (TR). */
00049 
00050 /*  Arguments */
00051 /*  ========= */
00052 
00053 /*  TRANSR   (input) CHARACTER */
00054 /*          = 'N':  ARF is in Normal format; */
00055 /*          = 'T':  ARF is in Transpose format. */
00056 
00057 /*  UPLO    (input) CHARACTER */
00058 /*          = 'U':  A is upper triangular; */
00059 /*          = 'L':  A is lower triangular. */
00060 
00061 /*  N       (input) INTEGER */
00062 /*          The order of the matrices ARF and A. N >= 0. */
00063 
00064 /*  ARF     (input) DOUBLE PRECISION array, dimension (N*(N+1)/2). */
00065 /*          On entry, the upper (if UPLO = 'U') or lower (if UPLO = 'L') */
00066 /*          matrix A in RFP format. See the "Notes" below for more */
00067 /*          details. */
00068 
00069 /*  A       (output) DOUBLE PRECISION array, dimension (LDA,N) */
00070 /*          On exit, the triangular matrix A.  If UPLO = 'U', the */
00071 /*          leading N-by-N upper triangular part of the array A contains */
00072 /*          the upper triangular matrix, and the strictly lower */
00073 /*          triangular part of A is not referenced.  If UPLO = 'L', the */
00074 /*          leading N-by-N lower triangular part of the array A contains */
00075 /*          the lower triangular matrix, and the strictly upper */
00076 /*          triangular part of A is not referenced. */
00077 
00078 /*  LDA     (input) INTEGER */
00079 /*          The leading dimension of the array A.  LDA >= max(1,N). */
00080 
00081 /*  INFO    (output) INTEGER */
00082 /*          = 0:  successful exit */
00083 /*          < 0:  if INFO = -i, the i-th argument had an illegal value */
00084 
00085 /*  Notes */
00086 /*  ===== */
00087 
00088 /*  We first consider Rectangular Full Packed (RFP) Format when N is */
00089 /*  even. We give an example where N = 6. */
00090 
00091 /*      AP is Upper             AP is Lower */
00092 
00093 /*   00 01 02 03 04 05       00 */
00094 /*      11 12 13 14 15       10 11 */
00095 /*         22 23 24 25       20 21 22 */
00096 /*            33 34 35       30 31 32 33 */
00097 /*               44 45       40 41 42 43 44 */
00098 /*                  55       50 51 52 53 54 55 */
00099 
00100 
00101 /*  Let TRANSR = 'N'. RFP holds AP as follows: */
00102 /*  For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last */
00103 /*  three columns of AP upper. The lower triangle A(4:6,0:2) consists of */
00104 /*  the transpose of the first three columns of AP upper. */
00105 /*  For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first */
00106 /*  three columns of AP lower. The upper triangle A(0:2,0:2) consists of */
00107 /*  the transpose of the last three columns of AP lower. */
00108 /*  This covers the case N even and TRANSR = 'N'. */
00109 
00110 /*         RFP A                   RFP A */
00111 
00112 /*        03 04 05                33 43 53 */
00113 /*        13 14 15                00 44 54 */
00114 /*        23 24 25                10 11 55 */
00115 /*        33 34 35                20 21 22 */
00116 /*        00 44 45                30 31 32 */
00117 /*        01 11 55                40 41 42 */
00118 /*        02 12 22                50 51 52 */
00119 
00120 /*  Now let TRANSR = 'T'. RFP A in both UPLO cases is just the */
00121 /*  transpose of RFP A above. One therefore gets: */
00122 
00123 
00124 /*           RFP A                   RFP A */
00125 
00126 /*     03 13 23 33 00 01 02    33 00 10 20 30 40 50 */
00127 /*     04 14 24 34 44 11 12    43 44 11 21 31 41 51 */
00128 /*     05 15 25 35 45 55 22    53 54 55 22 32 42 52 */
00129 
00130 
00131 /*  We first consider Rectangular Full Packed (RFP) Format when N is */
00132 /*  odd. We give an example where N = 5. */
00133 
00134 /*     AP is Upper                 AP is Lower */
00135 
00136 /*   00 01 02 03 04              00 */
00137 /*      11 12 13 14              10 11 */
00138 /*         22 23 24              20 21 22 */
00139 /*            33 34              30 31 32 33 */
00140 /*               44              40 41 42 43 44 */
00141 
00142 
00143 /*  Let TRANSR = 'N'. RFP holds AP as follows: */
00144 /*  For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last */
00145 /*  three columns of AP upper. The lower triangle A(3:4,0:1) consists of */
00146 /*  the transpose of the first two columns of AP upper. */
00147 /*  For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first */
00148 /*  three columns of AP lower. The upper triangle A(0:1,1:2) consists of */
00149 /*  the transpose of the last two columns of AP lower. */
00150 /*  This covers the case N odd and TRANSR = 'N'. */
00151 
00152 /*         RFP A                   RFP A */
00153 
00154 /*        02 03 04                00 33 43 */
00155 /*        12 13 14                10 11 44 */
00156 /*        22 23 24                20 21 22 */
00157 /*        00 33 34                30 31 32 */
00158 /*        01 11 44                40 41 42 */
00159 
00160 /*  Now let TRANSR = 'T'. RFP A in both UPLO cases is just the */
00161 /*  transpose of RFP A above. One therefore gets: */
00162 
00163 /*           RFP A                   RFP A */
00164 
00165 /*     02 12 22 00 01             00 10 20 30 40 50 */
00166 /*     03 13 23 33 11             33 11 21 31 41 51 */
00167 /*     04 14 24 34 44             43 44 22 32 42 52 */
00168 
00169 /*  Reference */
00170 /*  ========= */
00171 
00172 /*  ===================================================================== */
00173 
00174 /*     .. */
00175 /*     .. Local Scalars .. */
00176 /*     .. */
00177 /*     .. External Functions .. */
00178 /*     .. */
00179 /*     .. External Subroutines .. */
00180 /*     .. */
00181 /*     .. Intrinsic Functions .. */
00182 /*     .. */
00183 /*     .. Executable Statements .. */
00184 
00185 /*     Test the input parameters. */
00186 
00187     /* Parameter adjustments */
00188     a_dim1 = *lda - 1 - 0 + 1;
00189     a_offset = 0 + a_dim1 * 0;
00190     a -= a_offset;
00191 
00192     /* Function Body */
00193     *info = 0;
00194     normaltransr = lsame_(transr, "N");
00195     lower = lsame_(uplo, "L");
00196     if (! normaltransr && ! lsame_(transr, "T")) {
00197         *info = -1;
00198     } else if (! lower && ! lsame_(uplo, "U")) {
00199         *info = -2;
00200     } else if (*n < 0) {
00201         *info = -3;
00202     } else if (*lda < max(1,*n)) {
00203         *info = -6;
00204     }
00205     if (*info != 0) {
00206         i__1 = -(*info);
00207         xerbla_("DTFTTR", &i__1);
00208         return 0;
00209     }
00210 
00211 /*     Quick return if possible */
00212 
00213     if (*n <= 1) {
00214         if (*n == 1) {
00215             a[0] = arf[0];
00216         }
00217         return 0;
00218     }
00219 
00220 /*     Size of array ARF(0:nt-1) */
00221 
00222     nt = *n * (*n + 1) / 2;
00223 
00224 /*     set N1 and N2 depending on LOWER: for N even N1=N2=K */
00225 
00226     if (lower) {
00227         n2 = *n / 2;
00228         n1 = *n - n2;
00229     } else {
00230         n1 = *n / 2;
00231         n2 = *n - n1;
00232     }
00233 
00234 /*     If N is odd, set NISODD = .TRUE., LDA=N+1 and A is (N+1)--by--K2. */
00235 /*     If N is even, set K = N/2 and NISODD = .FALSE., LDA=N and A is */
00236 /*     N--by--(N+1)/2. */
00237 
00238     if (*n % 2 == 0) {
00239         k = *n / 2;
00240         nisodd = FALSE_;
00241         if (! lower) {
00242             np1x2 = *n + *n + 2;
00243         }
00244     } else {
00245         nisodd = TRUE_;
00246         if (! lower) {
00247             nx2 = *n + *n;
00248         }
00249     }
00250 
00251     if (nisodd) {
00252 
00253 /*        N is odd */
00254 
00255         if (normaltransr) {
00256 
00257 /*           N is odd and TRANSR = 'N' */
00258 
00259             if (lower) {
00260 
00261 /*              N is odd, TRANSR = 'N', and UPLO = 'L' */
00262 
00263                 ij = 0;
00264                 i__1 = n2;
00265                 for (j = 0; j <= i__1; ++j) {
00266                     i__2 = n2 + j;
00267                     for (i__ = n1; i__ <= i__2; ++i__) {
00268                         a[n2 + j + i__ * a_dim1] = arf[ij];
00269                         ++ij;
00270                     }
00271                     i__2 = *n - 1;
00272                     for (i__ = j; i__ <= i__2; ++i__) {
00273                         a[i__ + j * a_dim1] = arf[ij];
00274                         ++ij;
00275                     }
00276                 }
00277 
00278             } else {
00279 
00280 /*              N is odd, TRANSR = 'N', and UPLO = 'U' */
00281 
00282                 ij = nt - *n;
00283                 i__1 = n1;
00284                 for (j = *n - 1; j >= i__1; --j) {
00285                     i__2 = j;
00286                     for (i__ = 0; i__ <= i__2; ++i__) {
00287                         a[i__ + j * a_dim1] = arf[ij];
00288                         ++ij;
00289                     }
00290                     i__2 = n1 - 1;
00291                     for (l = j - n1; l <= i__2; ++l) {
00292                         a[j - n1 + l * a_dim1] = arf[ij];
00293                         ++ij;
00294                     }
00295                     ij -= nx2;
00296                 }
00297 
00298             }
00299 
00300         } else {
00301 
00302 /*           N is odd and TRANSR = 'T' */
00303 
00304             if (lower) {
00305 
00306 /*              N is odd, TRANSR = 'T', and UPLO = 'L' */
00307 
00308                 ij = 0;
00309                 i__1 = n2 - 1;
00310                 for (j = 0; j <= i__1; ++j) {
00311                     i__2 = j;
00312                     for (i__ = 0; i__ <= i__2; ++i__) {
00313                         a[j + i__ * a_dim1] = arf[ij];
00314                         ++ij;
00315                     }
00316                     i__2 = *n - 1;
00317                     for (i__ = n1 + j; i__ <= i__2; ++i__) {
00318                         a[i__ + (n1 + j) * a_dim1] = arf[ij];
00319                         ++ij;
00320                     }
00321                 }
00322                 i__1 = *n - 1;
00323                 for (j = n2; j <= i__1; ++j) {
00324                     i__2 = n1 - 1;
00325                     for (i__ = 0; i__ <= i__2; ++i__) {
00326                         a[j + i__ * a_dim1] = arf[ij];
00327                         ++ij;
00328                     }
00329                 }
00330 
00331             } else {
00332 
00333 /*              N is odd, TRANSR = 'T', and UPLO = 'U' */
00334 
00335                 ij = 0;
00336                 i__1 = n1;
00337                 for (j = 0; j <= i__1; ++j) {
00338                     i__2 = *n - 1;
00339                     for (i__ = n1; i__ <= i__2; ++i__) {
00340                         a[j + i__ * a_dim1] = arf[ij];
00341                         ++ij;
00342                     }
00343                 }
00344                 i__1 = n1 - 1;
00345                 for (j = 0; j <= i__1; ++j) {
00346                     i__2 = j;
00347                     for (i__ = 0; i__ <= i__2; ++i__) {
00348                         a[i__ + j * a_dim1] = arf[ij];
00349                         ++ij;
00350                     }
00351                     i__2 = *n - 1;
00352                     for (l = n2 + j; l <= i__2; ++l) {
00353                         a[n2 + j + l * a_dim1] = arf[ij];
00354                         ++ij;
00355                     }
00356                 }
00357 
00358             }
00359 
00360         }
00361 
00362     } else {
00363 
00364 /*        N is even */
00365 
00366         if (normaltransr) {
00367 
00368 /*           N is even and TRANSR = 'N' */
00369 
00370             if (lower) {
00371 
00372 /*              N is even, TRANSR = 'N', and UPLO = 'L' */
00373 
00374                 ij = 0;
00375                 i__1 = k - 1;
00376                 for (j = 0; j <= i__1; ++j) {
00377                     i__2 = k + j;
00378                     for (i__ = k; i__ <= i__2; ++i__) {
00379                         a[k + j + i__ * a_dim1] = arf[ij];
00380                         ++ij;
00381                     }
00382                     i__2 = *n - 1;
00383                     for (i__ = j; i__ <= i__2; ++i__) {
00384                         a[i__ + j * a_dim1] = arf[ij];
00385                         ++ij;
00386                     }
00387                 }
00388 
00389             } else {
00390 
00391 /*              N is even, TRANSR = 'N', and UPLO = 'U' */
00392 
00393                 ij = nt - *n - 1;
00394                 i__1 = k;
00395                 for (j = *n - 1; j >= i__1; --j) {
00396                     i__2 = j;
00397                     for (i__ = 0; i__ <= i__2; ++i__) {
00398                         a[i__ + j * a_dim1] = arf[ij];
00399                         ++ij;
00400                     }
00401                     i__2 = k - 1;
00402                     for (l = j - k; l <= i__2; ++l) {
00403                         a[j - k + l * a_dim1] = arf[ij];
00404                         ++ij;
00405                     }
00406                     ij -= np1x2;
00407                 }
00408 
00409             }
00410 
00411         } else {
00412 
00413 /*           N is even and TRANSR = 'T' */
00414 
00415             if (lower) {
00416 
00417 /*              N is even, TRANSR = 'T', and UPLO = 'L' */
00418 
00419                 ij = 0;
00420                 j = k;
00421                 i__1 = *n - 1;
00422                 for (i__ = k; i__ <= i__1; ++i__) {
00423                     a[i__ + j * a_dim1] = arf[ij];
00424                     ++ij;
00425                 }
00426                 i__1 = k - 2;
00427                 for (j = 0; j <= i__1; ++j) {
00428                     i__2 = j;
00429                     for (i__ = 0; i__ <= i__2; ++i__) {
00430                         a[j + i__ * a_dim1] = arf[ij];
00431                         ++ij;
00432                     }
00433                     i__2 = *n - 1;
00434                     for (i__ = k + 1 + j; i__ <= i__2; ++i__) {
00435                         a[i__ + (k + 1 + j) * a_dim1] = arf[ij];
00436                         ++ij;
00437                     }
00438                 }
00439                 i__1 = *n - 1;
00440                 for (j = k - 1; j <= i__1; ++j) {
00441                     i__2 = k - 1;
00442                     for (i__ = 0; i__ <= i__2; ++i__) {
00443                         a[j + i__ * a_dim1] = arf[ij];
00444                         ++ij;
00445                     }
00446                 }
00447 
00448             } else {
00449 
00450 /*              N is even, TRANSR = 'T', and UPLO = 'U' */
00451 
00452                 ij = 0;
00453                 i__1 = k;
00454                 for (j = 0; j <= i__1; ++j) {
00455                     i__2 = *n - 1;
00456                     for (i__ = k; i__ <= i__2; ++i__) {
00457                         a[j + i__ * a_dim1] = arf[ij];
00458                         ++ij;
00459                     }
00460                 }
00461                 i__1 = k - 2;
00462                 for (j = 0; j <= i__1; ++j) {
00463                     i__2 = j;
00464                     for (i__ = 0; i__ <= i__2; ++i__) {
00465                         a[i__ + j * a_dim1] = arf[ij];
00466                         ++ij;
00467                     }
00468                     i__2 = *n - 1;
00469                     for (l = k + 1 + j; l <= i__2; ++l) {
00470                         a[k + 1 + j + l * a_dim1] = arf[ij];
00471                         ++ij;
00472                     }
00473                 }
00474 /*              Note that here, on exit of the loop, J = K-1 */
00475                 i__1 = j;
00476                 for (i__ = 0; i__ <= i__1; ++i__) {
00477                     a[i__ + j * a_dim1] = arf[ij];
00478                     ++ij;
00479                 }
00480 
00481             }
00482 
00483         }
00484 
00485     }
00486 
00487     return 0;
00488 
00489 /*     End of DTFTTR */
00490 
00491 } /* dtfttr_ */


swiftnav
Author(s):
autogenerated on Sat Jun 8 2019 18:55:49