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


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