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


swiftnav
Author(s):
autogenerated on Sat Jun 8 2019 18:56:15