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


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