swiftnav: ctrttf.c Source File

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