ztpttf.c
Go to the documentation of this file.
00001 /* ztpttf.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 ztpttf_(char *transr, char *uplo, integer *n, 
00017         doublecomplex *ap, doublecomplex *arf, integer *info)
00018 {
00019     /* System generated locals */
00020     integer i__1, i__2, i__3, i__4;
00021     doublecomplex z__1;
00022 
00023     /* Builtin functions */
00024     void d_cnjg(doublecomplex *, doublecomplex *);
00025 
00026     /* Local variables */
00027     integer i__, j, k, n1, n2, ij, jp, js, nt, lda, ijp;
00028     logical normaltransr;
00029     extern logical lsame_(char *, char *);
00030     logical lower;
00031     extern /* Subroutine */ int xerbla_(char *, integer *);
00032     logical nisodd;
00033 
00034 
00035 /*  -- LAPACK routine (version 3.2)                                    -- */
00036 
00037 /*  -- Contributed by Fred Gustavson of the IBM Watson Research Center -- */
00038 /*  -- November 2008                                                   -- */
00039 
00040 /*  -- LAPACK is a software package provided by Univ. of Tennessee,    -- */
00041 /*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */
00042 
00043 /*     .. */
00044 /*     .. Scalar Arguments .. */
00045 /*     .. */
00046 /*     .. Array Arguments .. */
00047 
00048 /*  Purpose */
00049 /*  ======= */
00050 
00051 /*  ZTPTTF copies a triangular matrix A from standard packed format (TP) */
00052 /*  to rectangular full packed format (TF). */
00053 
00054 /*  Arguments */
00055 /*  ========= */
00056 
00057 /*  TRANSR   (input) CHARACTER */
00058 /*          = 'N':  ARF in Normal format is wanted; */
00059 /*          = 'C':  ARF in Conjugate-transpose format is wanted. */
00060 
00061 /*  UPLO    (input) CHARACTER */
00062 /*          = 'U':  A is upper triangular; */
00063 /*          = 'L':  A is lower triangular. */
00064 
00065 /*  N       (input) INTEGER */
00066 /*          The order of the matrix A.  N >= 0. */
00067 
00068 /*  AP      (input) COMPLEX*16 array, dimension ( N*(N+1)/2 ), */
00069 /*          On entry, the upper or lower triangular matrix A, packed */
00070 /*          columnwise in a linear array. The j-th column of A is stored */
00071 /*          in the array AP as follows: */
00072 /*          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j; */
00073 /*          if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n. */
00074 
00075 /*  ARF     (output) COMPLEX*16 array, dimension ( N*(N+1)/2 ), */
00076 /*          On exit, the upper or lower triangular matrix A stored in */
00077 /*          RFP format. For a further discussion see Notes below. */
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 Standard Packed Format when N is even. */
00087 /*  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 /*  conjugate-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 /*  conjugate-transpose of the last three columns of AP lower. */
00106 /*  To denote conjugate we place -- above the element. This covers the */
00107 /*  case N even and TRANSR = 'N'. */
00108 
00109 /*         RFP A                   RFP A */
00110 
00111 /*                                -- -- -- */
00112 /*        03 04 05                33 43 53 */
00113 /*                                   -- -- */
00114 /*        13 14 15                00 44 54 */
00115 /*                                      -- */
00116 /*        23 24 25                10 11 55 */
00117 
00118 /*        33 34 35                20 21 22 */
00119 /*        -- */
00120 /*        00 44 45                30 31 32 */
00121 /*        -- -- */
00122 /*        01 11 55                40 41 42 */
00123 /*        -- -- -- */
00124 /*        02 12 22                50 51 52 */
00125 
00126 /*  Now let TRANSR = 'C'. RFP A in both UPLO cases is just the conjugate- */
00127 /*  transpose of RFP A above. One therefore gets: */
00128 
00129 
00130 /*           RFP A                   RFP A */
00131 
00132 /*     -- -- -- --                -- -- -- -- -- -- */
00133 /*     03 13 23 33 00 01 02    33 00 10 20 30 40 50 */
00134 /*     -- -- -- -- --                -- -- -- -- -- */
00135 /*     04 14 24 34 44 11 12    43 44 11 21 31 41 51 */
00136 /*     -- -- -- -- -- --                -- -- -- -- */
00137 /*     05 15 25 35 45 55 22    53 54 55 22 32 42 52 */
00138 
00139 
00140 /*  We next  consider Standard Packed Format when N is odd. */
00141 /*  We give an example where N = 5. */
00142 
00143 /*     AP is Upper                 AP is Lower */
00144 
00145 /*   00 01 02 03 04              00 */
00146 /*      11 12 13 14              10 11 */
00147 /*         22 23 24              20 21 22 */
00148 /*            33 34              30 31 32 33 */
00149 /*               44              40 41 42 43 44 */
00150 
00151 
00152 /*  Let TRANSR = 'N'. RFP holds AP as follows: */
00153 /*  For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last */
00154 /*  three columns of AP upper. The lower triangle A(3:4,0:1) consists of */
00155 /*  conjugate-transpose of the first two   columns of AP upper. */
00156 /*  For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first */
00157 /*  three columns of AP lower. The upper triangle A(0:1,1:2) consists of */
00158 /*  conjugate-transpose of the last two   columns of AP lower. */
00159 /*  To denote conjugate we place -- above the element. This covers the */
00160 /*  case N odd  and TRANSR = 'N'. */
00161 
00162 /*         RFP A                   RFP A */
00163 
00164 /*                                   -- -- */
00165 /*        02 03 04                00 33 43 */
00166 /*                                      -- */
00167 /*        12 13 14                10 11 44 */
00168 
00169 /*        22 23 24                20 21 22 */
00170 /*        -- */
00171 /*        00 33 34                30 31 32 */
00172 /*        -- -- */
00173 /*        01 11 44                40 41 42 */
00174 
00175 /*  Now let TRANSR = 'C'. RFP A in both UPLO cases is just the conjugate- */
00176 /*  transpose of RFP A above. One therefore gets: */
00177 
00178 
00179 /*           RFP A                   RFP A */
00180 
00181 /*     -- -- --                   -- -- -- -- -- -- */
00182 /*     02 12 22 00 01             00 10 20 30 40 50 */
00183 /*     -- -- -- --                   -- -- -- -- -- */
00184 /*     03 13 23 33 11             33 11 21 31 41 51 */
00185 /*     -- -- -- -- --                   -- -- -- -- */
00186 /*     04 14 24 34 44             43 44 22 32 42 52 */
00187 
00188 /*  ===================================================================== */
00189 
00190 /*     .. Parameters .. */
00191 /*     .. */
00192 /*     .. Local Scalars .. */
00193 /*     .. */
00194 /*     .. External Functions .. */
00195 /*     .. */
00196 /*     .. External Subroutines .. */
00197 /*     .. */
00198 /*     .. Intrinsic Functions .. */
00199 /*     .. */
00200 /*     .. Executable Statements .. */
00201 
00202 /*     Test the input parameters. */
00203 
00204     *info = 0;
00205     normaltransr = lsame_(transr, "N");
00206     lower = lsame_(uplo, "L");
00207     if (! normaltransr && ! lsame_(transr, "C")) {
00208         *info = -1;
00209     } else if (! lower && ! lsame_(uplo, "U")) {
00210         *info = -2;
00211     } else if (*n < 0) {
00212         *info = -3;
00213     }
00214     if (*info != 0) {
00215         i__1 = -(*info);
00216         xerbla_("ZTPTTF", &i__1);
00217         return 0;
00218     }
00219 
00220 /*     Quick return if possible */
00221 
00222     if (*n == 0) {
00223         return 0;
00224     }
00225 
00226     if (*n == 1) {
00227         if (normaltransr) {
00228             arf[0].r = ap[0].r, arf[0].i = ap[0].i;
00229         } else {
00230             d_cnjg(&z__1, ap);
00231             arf[0].r = z__1.r, arf[0].i = z__1.i;
00232         }
00233         return 0;
00234     }
00235 
00236 /*     Size of array ARF(0:NT-1) */
00237 
00238     nt = *n * (*n + 1) / 2;
00239 
00240 /*     Set N1 and N2 depending on LOWER */
00241 
00242     if (lower) {
00243         n2 = *n / 2;
00244         n1 = *n - n2;
00245     } else {
00246         n1 = *n / 2;
00247         n2 = *n - n1;
00248     }
00249 
00250 /*     If N is odd, set NISODD = .TRUE. */
00251 /*     If N is even, set K = N/2 and NISODD = .FALSE. */
00252 
00253 /*     set lda of ARF^C; ARF^C is (0:(N+1)/2-1,0:N-noe) */
00254 /*     where noe = 0 if n is even, noe = 1 if n is odd */
00255 
00256     if (*n % 2 == 0) {
00257         k = *n / 2;
00258         nisodd = FALSE_;
00259         lda = *n + 1;
00260     } else {
00261         nisodd = TRUE_;
00262         lda = *n;
00263     }
00264 
00265 /*     ARF^C has lda rows and n+1-noe cols */
00266 
00267     if (! normaltransr) {
00268         lda = (*n + 1) / 2;
00269     }
00270 
00271 /*     start execution: there are eight cases */
00272 
00273     if (nisodd) {
00274 
00275 /*        N is odd */
00276 
00277         if (normaltransr) {
00278 
00279 /*           N is odd and TRANSR = 'N' */
00280 
00281             if (lower) {
00282 
00283 /*             SRPA for LOWER, NORMAL and N is odd ( a(0:n-1,0:n1-1) ) */
00284 /*             T1 -> a(0,0), T2 -> a(0,1), S -> a(n1,0) */
00285 /*             T1 -> a(0), T2 -> a(n), S -> a(n1); lda = n */
00286 
00287                 ijp = 0;
00288                 jp = 0;
00289                 i__1 = n2;
00290                 for (j = 0; j <= i__1; ++j) {
00291                     i__2 = *n - 1;
00292                     for (i__ = j; i__ <= i__2; ++i__) {
00293                         ij = i__ + jp;
00294                         i__3 = ij;
00295                         i__4 = ijp;
00296                         arf[i__3].r = ap[i__4].r, arf[i__3].i = ap[i__4].i;
00297                         ++ijp;
00298                     }
00299                     jp += lda;
00300                 }
00301                 i__1 = n2 - 1;
00302                 for (i__ = 0; i__ <= i__1; ++i__) {
00303                     i__2 = n2;
00304                     for (j = i__ + 1; j <= i__2; ++j) {
00305                         ij = i__ + j * lda;
00306                         i__3 = ij;
00307                         d_cnjg(&z__1, &ap[ijp]);
00308                         arf[i__3].r = z__1.r, arf[i__3].i = z__1.i;
00309                         ++ijp;
00310                     }
00311                 }
00312 
00313             } else {
00314 
00315 /*             SRPA for UPPER, NORMAL and N is odd ( a(0:n-1,0:n2-1) */
00316 /*             T1 -> a(n1+1,0), T2 -> a(n1,0), S -> a(0,0) */
00317 /*             T1 -> a(n2), T2 -> a(n1), S -> a(0) */
00318 
00319                 ijp = 0;
00320                 i__1 = n1 - 1;
00321                 for (j = 0; j <= i__1; ++j) {
00322                     ij = n2 + j;
00323                     i__2 = j;
00324                     for (i__ = 0; i__ <= i__2; ++i__) {
00325                         i__3 = ij;
00326                         d_cnjg(&z__1, &ap[ijp]);
00327                         arf[i__3].r = z__1.r, arf[i__3].i = z__1.i;
00328                         ++ijp;
00329                         ij += lda;
00330                     }
00331                 }
00332                 js = 0;
00333                 i__1 = *n - 1;
00334                 for (j = n1; j <= i__1; ++j) {
00335                     ij = js;
00336                     i__2 = js + j;
00337                     for (ij = js; ij <= i__2; ++ij) {
00338                         i__3 = ij;
00339                         i__4 = ijp;
00340                         arf[i__3].r = ap[i__4].r, arf[i__3].i = ap[i__4].i;
00341                         ++ijp;
00342                     }
00343                     js += lda;
00344                 }
00345 
00346             }
00347 
00348         } else {
00349 
00350 /*           N is odd and TRANSR = 'C' */
00351 
00352             if (lower) {
00353 
00354 /*              SRPA for LOWER, TRANSPOSE and N is odd */
00355 /*              T1 -> A(0,0) , T2 -> A(1,0) , S -> A(0,n1) */
00356 /*              T1 -> a(0+0) , T2 -> a(1+0) , S -> a(0+n1*n1); lda=n1 */
00357 
00358                 ijp = 0;
00359                 i__1 = n2;
00360                 for (i__ = 0; i__ <= i__1; ++i__) {
00361                     i__2 = *n * lda - 1;
00362                     i__3 = lda;
00363                     for (ij = i__ * (lda + 1); i__3 < 0 ? ij >= i__2 : ij <= 
00364                             i__2; ij += i__3) {
00365                         i__4 = ij;
00366                         d_cnjg(&z__1, &ap[ijp]);
00367                         arf[i__4].r = z__1.r, arf[i__4].i = z__1.i;
00368                         ++ijp;
00369                     }
00370                 }
00371                 js = 1;
00372                 i__1 = n2 - 1;
00373                 for (j = 0; j <= i__1; ++j) {
00374                     i__3 = js + n2 - j - 1;
00375                     for (ij = js; ij <= i__3; ++ij) {
00376                         i__2 = ij;
00377                         i__4 = ijp;
00378                         arf[i__2].r = ap[i__4].r, arf[i__2].i = ap[i__4].i;
00379                         ++ijp;
00380                     }
00381                     js = js + lda + 1;
00382                 }
00383 
00384             } else {
00385 
00386 /*              SRPA for UPPER, TRANSPOSE and N is odd */
00387 /*              T1 -> A(0,n1+1), T2 -> A(0,n1), S -> A(0,0) */
00388 /*              T1 -> a(n2*n2), T2 -> a(n1*n2), S -> a(0); lda = n2 */
00389 
00390                 ijp = 0;
00391                 js = n2 * lda;
00392                 i__1 = n1 - 1;
00393                 for (j = 0; j <= i__1; ++j) {
00394                     i__3 = js + j;
00395                     for (ij = js; ij <= i__3; ++ij) {
00396                         i__2 = ij;
00397                         i__4 = ijp;
00398                         arf[i__2].r = ap[i__4].r, arf[i__2].i = ap[i__4].i;
00399                         ++ijp;
00400                     }
00401                     js += lda;
00402                 }
00403                 i__1 = n1;
00404                 for (i__ = 0; i__ <= i__1; ++i__) {
00405                     i__3 = i__ + (n1 + i__) * lda;
00406                     i__2 = lda;
00407                     for (ij = i__; i__2 < 0 ? ij >= i__3 : ij <= i__3; ij += 
00408                             i__2) {
00409                         i__4 = ij;
00410                         d_cnjg(&z__1, &ap[ijp]);
00411                         arf[i__4].r = z__1.r, arf[i__4].i = z__1.i;
00412                         ++ijp;
00413                     }
00414                 }
00415 
00416             }
00417 
00418         }
00419 
00420     } else {
00421 
00422 /*        N is even */
00423 
00424         if (normaltransr) {
00425 
00426 /*           N is even and TRANSR = 'N' */
00427 
00428             if (lower) {
00429 
00430 /*              SRPA for LOWER, NORMAL, and N is even ( a(0:n,0:k-1) ) */
00431 /*              T1 -> a(1,0), T2 -> a(0,0), S -> a(k+1,0) */
00432 /*              T1 -> a(1), T2 -> a(0), S -> a(k+1) */
00433 
00434                 ijp = 0;
00435                 jp = 0;
00436                 i__1 = k - 1;
00437                 for (j = 0; j <= i__1; ++j) {
00438                     i__2 = *n - 1;
00439                     for (i__ = j; i__ <= i__2; ++i__) {
00440                         ij = i__ + 1 + jp;
00441                         i__3 = ij;
00442                         i__4 = ijp;
00443                         arf[i__3].r = ap[i__4].r, arf[i__3].i = ap[i__4].i;
00444                         ++ijp;
00445                     }
00446                     jp += lda;
00447                 }
00448                 i__1 = k - 1;
00449                 for (i__ = 0; i__ <= i__1; ++i__) {
00450                     i__2 = k - 1;
00451                     for (j = i__; j <= i__2; ++j) {
00452                         ij = i__ + j * lda;
00453                         i__3 = ij;
00454                         d_cnjg(&z__1, &ap[ijp]);
00455                         arf[i__3].r = z__1.r, arf[i__3].i = z__1.i;
00456                         ++ijp;
00457                     }
00458                 }
00459 
00460             } else {
00461 
00462 /*              SRPA for UPPER, NORMAL, and N is even ( a(0:n,0:k-1) ) */
00463 /*              T1 -> a(k+1,0) ,  T2 -> a(k,0),   S -> a(0,0) */
00464 /*              T1 -> a(k+1), T2 -> a(k), S -> a(0) */
00465 
00466                 ijp = 0;
00467                 i__1 = k - 1;
00468                 for (j = 0; j <= i__1; ++j) {
00469                     ij = k + 1 + j;
00470                     i__2 = j;
00471                     for (i__ = 0; i__ <= i__2; ++i__) {
00472                         i__3 = ij;
00473                         d_cnjg(&z__1, &ap[ijp]);
00474                         arf[i__3].r = z__1.r, arf[i__3].i = z__1.i;
00475                         ++ijp;
00476                         ij += lda;
00477                     }
00478                 }
00479                 js = 0;
00480                 i__1 = *n - 1;
00481                 for (j = k; j <= i__1; ++j) {
00482                     ij = js;
00483                     i__2 = js + j;
00484                     for (ij = js; ij <= i__2; ++ij) {
00485                         i__3 = ij;
00486                         i__4 = ijp;
00487                         arf[i__3].r = ap[i__4].r, arf[i__3].i = ap[i__4].i;
00488                         ++ijp;
00489                     }
00490                     js += lda;
00491                 }
00492 
00493             }
00494 
00495         } else {
00496 
00497 /*           N is even and TRANSR = 'C' */
00498 
00499             if (lower) {
00500 
00501 /*              SRPA for LOWER, TRANSPOSE and N is even (see paper) */
00502 /*              T1 -> B(0,1), T2 -> B(0,0), S -> B(0,k+1) */
00503 /*              T1 -> a(0+k), T2 -> a(0+0), S -> a(0+k*(k+1)); lda=k */
00504 
00505                 ijp = 0;
00506                 i__1 = k - 1;
00507                 for (i__ = 0; i__ <= i__1; ++i__) {
00508                     i__2 = (*n + 1) * lda - 1;
00509                     i__3 = lda;
00510                     for (ij = i__ + (i__ + 1) * lda; i__3 < 0 ? ij >= i__2 : 
00511                             ij <= i__2; ij += i__3) {
00512                         i__4 = ij;
00513                         d_cnjg(&z__1, &ap[ijp]);
00514                         arf[i__4].r = z__1.r, arf[i__4].i = z__1.i;
00515                         ++ijp;
00516                     }
00517                 }
00518                 js = 0;
00519                 i__1 = k - 1;
00520                 for (j = 0; j <= i__1; ++j) {
00521                     i__3 = js + k - j - 1;
00522                     for (ij = js; ij <= i__3; ++ij) {
00523                         i__2 = ij;
00524                         i__4 = ijp;
00525                         arf[i__2].r = ap[i__4].r, arf[i__2].i = ap[i__4].i;
00526                         ++ijp;
00527                     }
00528                     js = js + lda + 1;
00529                 }
00530 
00531             } else {
00532 
00533 /*              SRPA for UPPER, TRANSPOSE and N is even (see paper) */
00534 /*              T1 -> B(0,k+1),     T2 -> B(0,k),   S -> B(0,0) */
00535 /*              T1 -> a(0+k*(k+1)), T2 -> a(0+k*k), S -> a(0+0)); lda=k */
00536 
00537                 ijp = 0;
00538                 js = (k + 1) * lda;
00539                 i__1 = k - 1;
00540                 for (j = 0; j <= i__1; ++j) {
00541                     i__3 = js + j;
00542                     for (ij = js; ij <= i__3; ++ij) {
00543                         i__2 = ij;
00544                         i__4 = ijp;
00545                         arf[i__2].r = ap[i__4].r, arf[i__2].i = ap[i__4].i;
00546                         ++ijp;
00547                     }
00548                     js += lda;
00549                 }
00550                 i__1 = k - 1;
00551                 for (i__ = 0; i__ <= i__1; ++i__) {
00552                     i__3 = i__ + (k + i__) * lda;
00553                     i__2 = lda;
00554                     for (ij = i__; i__2 < 0 ? ij >= i__3 : ij <= i__3; ij += 
00555                             i__2) {
00556                         i__4 = ij;
00557                         d_cnjg(&z__1, &ap[ijp]);
00558                         arf[i__4].r = z__1.r, arf[i__4].i = z__1.i;
00559                         ++ijp;
00560                     }
00561                 }
00562 
00563             }
00564 
00565         }
00566 
00567     }
00568 
00569     return 0;
00570 
00571 /*     End of ZTPTTF */
00572 
00573 } /* ztpttf_ */

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