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


swiftnav
Author(s):
autogenerated on Sat Jun 8 2019 18:55:32