zlatsp.c
Go to the documentation of this file.
00001 /* zlatsp.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 /* Table of constant values */
00017 
00018 static integer c__5 = 5;
00019 static integer c__2 = 2;
00020 
00021 /* Subroutine */ int zlatsp_(char *uplo, integer *n, doublecomplex *x, 
00022         integer *iseed)
00023 {
00024     /* System generated locals */
00025     integer i__1, i__2, i__3;
00026     doublecomplex z__1, z__2, z__3;
00027 
00028     /* Builtin functions */
00029     double sqrt(doublereal), z_abs(doublecomplex *);
00030 
00031     /* Local variables */
00032     doublecomplex a, b, c__;
00033     integer j;
00034     doublecomplex r__;
00035     integer n5, jj;
00036     doublereal beta, alpha, alpha3;
00037     extern /* Double Complex */ VOID zlarnd_(doublecomplex *, integer *, 
00038             integer *);
00039 
00040 
00041 /*  -- LAPACK auxiliary test routine (version 3.1) -- */
00042 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00043 /*     November 2006 */
00044 
00045 /*     .. Scalar Arguments .. */
00046 /*     .. */
00047 /*     .. Array Arguments .. */
00048 /*     .. */
00049 
00050 /*  Purpose */
00051 /*  ======= */
00052 
00053 /*  ZLATSP generates a special test matrix for the complex symmetric */
00054 /*  (indefinite) factorization for packed matrices.  The pivot blocks of */
00055 /*  the generated matrix will be in the following order: */
00056 /*     2x2 pivot block, non diagonalizable */
00057 /*     1x1 pivot block */
00058 /*     2x2 pivot block, diagonalizable */
00059 /*     (cycle repeats) */
00060 /*  A row interchange is required for each non-diagonalizable 2x2 block. */
00061 
00062 /*  Arguments */
00063 /*  ========= */
00064 
00065 /*  UPLO    (input) CHARACTER */
00066 /*          Specifies whether the generated matrix is to be upper or */
00067 /*          lower triangular. */
00068 /*          = 'U':  Upper triangular */
00069 /*          = 'L':  Lower triangular */
00070 
00071 /*  N       (input) INTEGER */
00072 /*          The dimension of the matrix to be generated. */
00073 
00074 /*  X       (output) COMPLEX*16 array, dimension (N*(N+1)/2) */
00075 /*          The generated matrix in packed storage format.  The matrix */
00076 /*          consists of 3x3 and 2x2 diagonal blocks which result in the */
00077 /*          pivot sequence given above.  The matrix outside these */
00078 /*          diagonal blocks is zero. */
00079 
00080 /*  ISEED   (input/output) INTEGER array, dimension (4) */
00081 /*          On entry, the seed for the random number generator.  The last */
00082 /*          of the four integers must be odd.  (modified on exit) */
00083 
00084 /*  ===================================================================== */
00085 
00086 /*     .. Parameters .. */
00087 /*     .. */
00088 /*     .. Local Scalars .. */
00089 /*     .. */
00090 /*     .. External Functions .. */
00091 /*     .. */
00092 /*     .. Intrinsic Functions .. */
00093 /*     .. */
00094 /*     .. Executable Statements .. */
00095 
00096 /*     Initialize constants */
00097 
00098     /* Parameter adjustments */
00099     --iseed;
00100     --x;
00101 
00102     /* Function Body */
00103     alpha = (sqrt(17.) + 1.) / 8.;
00104     beta = alpha - .001;
00105     alpha3 = alpha * alpha * alpha;
00106 
00107 /*     Fill the matrix with zeros. */
00108 
00109     i__1 = *n * (*n + 1) / 2;
00110     for (j = 1; j <= i__1; ++j) {
00111         i__2 = j;
00112         x[i__2].r = 0., x[i__2].i = 0.;
00113 /* L10: */
00114     }
00115 
00116 /*     UPLO = 'U':  Upper triangular storage */
00117 
00118     if (*(unsigned char *)uplo == 'U') {
00119         n5 = *n / 5;
00120         n5 = *n - n5 * 5 + 1;
00121 
00122         jj = *n * (*n + 1) / 2;
00123         i__1 = n5;
00124         for (j = *n; j >= i__1; j += -5) {
00125             zlarnd_(&z__2, &c__5, &iseed[1]);
00126             z__1.r = alpha3 * z__2.r, z__1.i = alpha3 * z__2.i;
00127             a.r = z__1.r, a.i = z__1.i;
00128             zlarnd_(&z__2, &c__5, &iseed[1]);
00129             z__1.r = z__2.r / alpha, z__1.i = z__2.i / alpha;
00130             b.r = z__1.r, b.i = z__1.i;
00131             z__3.r = b.r * 2., z__3.i = b.i * 2.;
00132             z__2.r = z__3.r * 0. - z__3.i * 1., z__2.i = z__3.r * 1. + z__3.i 
00133                     * 0.;
00134             z__1.r = a.r - z__2.r, z__1.i = a.i - z__2.i;
00135             c__.r = z__1.r, c__.i = z__1.i;
00136             z__1.r = c__.r / beta, z__1.i = c__.i / beta;
00137             r__.r = z__1.r, r__.i = z__1.i;
00138             i__2 = jj;
00139             x[i__2].r = a.r, x[i__2].i = a.i;
00140             i__2 = jj - 2;
00141             x[i__2].r = b.r, x[i__2].i = b.i;
00142             jj -= j;
00143             i__2 = jj;
00144             zlarnd_(&z__1, &c__2, &iseed[1]);
00145             x[i__2].r = z__1.r, x[i__2].i = z__1.i;
00146             i__2 = jj - 1;
00147             x[i__2].r = r__.r, x[i__2].i = r__.i;
00148             jj -= j - 1;
00149             i__2 = jj;
00150             x[i__2].r = c__.r, x[i__2].i = c__.i;
00151             jj -= j - 2;
00152             i__2 = jj;
00153             zlarnd_(&z__1, &c__2, &iseed[1]);
00154             x[i__2].r = z__1.r, x[i__2].i = z__1.i;
00155             jj -= j - 3;
00156             i__2 = jj;
00157             zlarnd_(&z__1, &c__2, &iseed[1]);
00158             x[i__2].r = z__1.r, x[i__2].i = z__1.i;
00159             if (z_abs(&x[jj + (j - 3)]) > z_abs(&x[jj])) {
00160                 i__2 = jj + (j - 4);
00161                 i__3 = jj + (j - 3);
00162                 z__1.r = x[i__3].r * 2., z__1.i = x[i__3].i * 2.;
00163                 x[i__2].r = z__1.r, x[i__2].i = z__1.i;
00164             } else {
00165                 i__2 = jj + (j - 4);
00166                 i__3 = jj;
00167                 z__1.r = x[i__3].r * 2., z__1.i = x[i__3].i * 2.;
00168                 x[i__2].r = z__1.r, x[i__2].i = z__1.i;
00169             }
00170             jj -= j - 4;
00171 /* L20: */
00172         }
00173 
00174 /*        Clean-up for N not a multiple of 5. */
00175 
00176         j = n5 - 1;
00177         if (j > 2) {
00178             zlarnd_(&z__2, &c__5, &iseed[1]);
00179             z__1.r = alpha3 * z__2.r, z__1.i = alpha3 * z__2.i;
00180             a.r = z__1.r, a.i = z__1.i;
00181             zlarnd_(&z__2, &c__5, &iseed[1]);
00182             z__1.r = z__2.r / alpha, z__1.i = z__2.i / alpha;
00183             b.r = z__1.r, b.i = z__1.i;
00184             z__3.r = b.r * 2., z__3.i = b.i * 2.;
00185             z__2.r = z__3.r * 0. - z__3.i * 1., z__2.i = z__3.r * 1. + z__3.i 
00186                     * 0.;
00187             z__1.r = a.r - z__2.r, z__1.i = a.i - z__2.i;
00188             c__.r = z__1.r, c__.i = z__1.i;
00189             z__1.r = c__.r / beta, z__1.i = c__.i / beta;
00190             r__.r = z__1.r, r__.i = z__1.i;
00191             i__1 = jj;
00192             x[i__1].r = a.r, x[i__1].i = a.i;
00193             i__1 = jj - 2;
00194             x[i__1].r = b.r, x[i__1].i = b.i;
00195             jj -= j;
00196             i__1 = jj;
00197             zlarnd_(&z__1, &c__2, &iseed[1]);
00198             x[i__1].r = z__1.r, x[i__1].i = z__1.i;
00199             i__1 = jj - 1;
00200             x[i__1].r = r__.r, x[i__1].i = r__.i;
00201             jj -= j - 1;
00202             i__1 = jj;
00203             x[i__1].r = c__.r, x[i__1].i = c__.i;
00204             jj -= j - 2;
00205             j += -3;
00206         }
00207         if (j > 1) {
00208             i__1 = jj;
00209             zlarnd_(&z__1, &c__2, &iseed[1]);
00210             x[i__1].r = z__1.r, x[i__1].i = z__1.i;
00211             i__1 = jj - j;
00212             zlarnd_(&z__1, &c__2, &iseed[1]);
00213             x[i__1].r = z__1.r, x[i__1].i = z__1.i;
00214             if (z_abs(&x[jj]) > z_abs(&x[jj - j])) {
00215                 i__1 = jj - 1;
00216                 i__2 = jj;
00217                 z__1.r = x[i__2].r * 2., z__1.i = x[i__2].i * 2.;
00218                 x[i__1].r = z__1.r, x[i__1].i = z__1.i;
00219             } else {
00220                 i__1 = jj - 1;
00221                 i__2 = jj - j;
00222                 z__1.r = x[i__2].r * 2., z__1.i = x[i__2].i * 2.;
00223                 x[i__1].r = z__1.r, x[i__1].i = z__1.i;
00224             }
00225             jj = jj - j - (j - 1);
00226             j += -2;
00227         } else if (j == 1) {
00228             i__1 = jj;
00229             zlarnd_(&z__1, &c__2, &iseed[1]);
00230             x[i__1].r = z__1.r, x[i__1].i = z__1.i;
00231             --j;
00232         }
00233 
00234 /*     UPLO = 'L':  Lower triangular storage */
00235 
00236     } else {
00237         n5 = *n / 5;
00238         n5 *= 5;
00239 
00240         jj = 1;
00241         i__1 = n5;
00242         for (j = 1; j <= i__1; j += 5) {
00243             zlarnd_(&z__2, &c__5, &iseed[1]);
00244             z__1.r = alpha3 * z__2.r, z__1.i = alpha3 * z__2.i;
00245             a.r = z__1.r, a.i = z__1.i;
00246             zlarnd_(&z__2, &c__5, &iseed[1]);
00247             z__1.r = z__2.r / alpha, z__1.i = z__2.i / alpha;
00248             b.r = z__1.r, b.i = z__1.i;
00249             z__3.r = b.r * 2., z__3.i = b.i * 2.;
00250             z__2.r = z__3.r * 0. - z__3.i * 1., z__2.i = z__3.r * 1. + z__3.i 
00251                     * 0.;
00252             z__1.r = a.r - z__2.r, z__1.i = a.i - z__2.i;
00253             c__.r = z__1.r, c__.i = z__1.i;
00254             z__1.r = c__.r / beta, z__1.i = c__.i / beta;
00255             r__.r = z__1.r, r__.i = z__1.i;
00256             i__2 = jj;
00257             x[i__2].r = a.r, x[i__2].i = a.i;
00258             i__2 = jj + 2;
00259             x[i__2].r = b.r, x[i__2].i = b.i;
00260             jj += *n - j + 1;
00261             i__2 = jj;
00262             zlarnd_(&z__1, &c__2, &iseed[1]);
00263             x[i__2].r = z__1.r, x[i__2].i = z__1.i;
00264             i__2 = jj + 1;
00265             x[i__2].r = r__.r, x[i__2].i = r__.i;
00266             jj += *n - j;
00267             i__2 = jj;
00268             x[i__2].r = c__.r, x[i__2].i = c__.i;
00269             jj += *n - j - 1;
00270             i__2 = jj;
00271             zlarnd_(&z__1, &c__2, &iseed[1]);
00272             x[i__2].r = z__1.r, x[i__2].i = z__1.i;
00273             jj += *n - j - 2;
00274             i__2 = jj;
00275             zlarnd_(&z__1, &c__2, &iseed[1]);
00276             x[i__2].r = z__1.r, x[i__2].i = z__1.i;
00277             if (z_abs(&x[jj - (*n - j - 2)]) > z_abs(&x[jj])) {
00278                 i__2 = jj - (*n - j - 2) + 1;
00279                 i__3 = jj - (*n - j - 2);
00280                 z__1.r = x[i__3].r * 2., z__1.i = x[i__3].i * 2.;
00281                 x[i__2].r = z__1.r, x[i__2].i = z__1.i;
00282             } else {
00283                 i__2 = jj - (*n - j - 2) + 1;
00284                 i__3 = jj;
00285                 z__1.r = x[i__3].r * 2., z__1.i = x[i__3].i * 2.;
00286                 x[i__2].r = z__1.r, x[i__2].i = z__1.i;
00287             }
00288             jj += *n - j - 3;
00289 /* L30: */
00290         }
00291 
00292 /*        Clean-up for N not a multiple of 5. */
00293 
00294         j = n5 + 1;
00295         if (j < *n - 1) {
00296             zlarnd_(&z__2, &c__5, &iseed[1]);
00297             z__1.r = alpha3 * z__2.r, z__1.i = alpha3 * z__2.i;
00298             a.r = z__1.r, a.i = z__1.i;
00299             zlarnd_(&z__2, &c__5, &iseed[1]);
00300             z__1.r = z__2.r / alpha, z__1.i = z__2.i / alpha;
00301             b.r = z__1.r, b.i = z__1.i;
00302             z__3.r = b.r * 2., z__3.i = b.i * 2.;
00303             z__2.r = z__3.r * 0. - z__3.i * 1., z__2.i = z__3.r * 1. + z__3.i 
00304                     * 0.;
00305             z__1.r = a.r - z__2.r, z__1.i = a.i - z__2.i;
00306             c__.r = z__1.r, c__.i = z__1.i;
00307             z__1.r = c__.r / beta, z__1.i = c__.i / beta;
00308             r__.r = z__1.r, r__.i = z__1.i;
00309             i__1 = jj;
00310             x[i__1].r = a.r, x[i__1].i = a.i;
00311             i__1 = jj + 2;
00312             x[i__1].r = b.r, x[i__1].i = b.i;
00313             jj += *n - j + 1;
00314             i__1 = jj;
00315             zlarnd_(&z__1, &c__2, &iseed[1]);
00316             x[i__1].r = z__1.r, x[i__1].i = z__1.i;
00317             i__1 = jj + 1;
00318             x[i__1].r = r__.r, x[i__1].i = r__.i;
00319             jj += *n - j;
00320             i__1 = jj;
00321             x[i__1].r = c__.r, x[i__1].i = c__.i;
00322             jj += *n - j - 1;
00323             j += 3;
00324         }
00325         if (j < *n) {
00326             i__1 = jj;
00327             zlarnd_(&z__1, &c__2, &iseed[1]);
00328             x[i__1].r = z__1.r, x[i__1].i = z__1.i;
00329             i__1 = jj + (*n - j + 1);
00330             zlarnd_(&z__1, &c__2, &iseed[1]);
00331             x[i__1].r = z__1.r, x[i__1].i = z__1.i;
00332             if (z_abs(&x[jj]) > z_abs(&x[jj + (*n - j + 1)])) {
00333                 i__1 = jj + 1;
00334                 i__2 = jj;
00335                 z__1.r = x[i__2].r * 2., z__1.i = x[i__2].i * 2.;
00336                 x[i__1].r = z__1.r, x[i__1].i = z__1.i;
00337             } else {
00338                 i__1 = jj + 1;
00339                 i__2 = jj + (*n - j + 1);
00340                 z__1.r = x[i__2].r * 2., z__1.i = x[i__2].i * 2.;
00341                 x[i__1].r = z__1.r, x[i__1].i = z__1.i;
00342             }
00343             jj = jj + (*n - j + 1) + (*n - j);
00344             j += 2;
00345         } else if (j == *n) {
00346             i__1 = jj;
00347             zlarnd_(&z__1, &c__2, &iseed[1]);
00348             x[i__1].r = z__1.r, x[i__1].i = z__1.i;
00349             jj += *n - j + 1;
00350             ++j;
00351         }
00352     }
00353 
00354     return 0;
00355 
00356 /*     End of ZLATSP */
00357 
00358 } /* zlatsp_ */


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