claror.c

Go to the documentation of this file.

/* claror.f -- translated by f2c (version 20061008).
   You must link the resulting object file with libf2c:
        on Microsoft Windows system, link with libf2c.lib;
        on Linux or Unix systems, link with .../path/to/libf2c.a -lm
        or, if you install libf2c.a in a standard place, with -lf2c -lm
        -- in that order, at the end of the command line, as in
                cc *.o -lf2c -lm
        Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,

                http://www.netlib.org/f2c/libf2c.zip
*/

#include "f2c.h"
#include "blaswrap.h"

/* Table of constant values */

static complex c_b1 = {0.f,0.f};
static complex c_b2 = {1.f,0.f};
static integer c__3 = 3;
static integer c__1 = 1;

/* Subroutine */ int claror_(char *side, char *init, integer *m, integer *n, 
        complex *a, integer *lda, integer *iseed, complex *x, integer *info)
{
    /* System generated locals */
    integer a_dim1, a_offset, i__1, i__2, i__3;
    complex q__1, q__2;

    /* Builtin functions */
    double c_abs(complex *);
    void r_cnjg(complex *, complex *);

    /* Local variables */
    integer j, kbeg, jcol;
    real xabs;
    integer irow;
    extern /* Subroutine */ int cgerc_(integer *, integer *, complex *, 
            complex *, integer *, complex *, integer *, complex *, integer *),
             cscal_(integer *, complex *, complex *, integer *);
    extern logical lsame_(char *, char *);
    extern /* Subroutine */ int cgemv_(char *, integer *, integer *, complex *
, complex *, integer *, complex *, integer *, complex *, complex *
, integer *);
    complex csign;
    integer ixfrm, itype, nxfrm;
    real xnorm;
    extern doublereal scnrm2_(integer *, complex *, integer *);
    extern /* Subroutine */ int clacgv_(integer *, complex *, integer *);
    extern /* Complex */ VOID clarnd_(complex *, integer *, integer *);
    extern /* Subroutine */ int claset_(char *, integer *, integer *, complex 
            *, complex *, complex *, integer *), xerbla_(char *, 
            integer *);
    real factor;
    complex xnorms;


/*  -- LAPACK auxiliary test routine (version 3.1) -- */
/*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
/*     November 2006 */

/*     .. Scalar Arguments .. */
/*     .. */
/*     .. Array Arguments .. */
/*     .. */

/*  Purpose */
/*  ======= */

/*     CLAROR pre- or post-multiplies an M by N matrix A by a random */
/*     unitary matrix U, overwriting A. A may optionally be */
/*     initialized to the identity matrix before multiplying by U. */
/*     U is generated using the method of G.W. Stewart */
/*     ( SIAM J. Numer. Anal. 17, 1980, pp. 403-409 ). */
/*     (BLAS-2 version) */

/*  Arguments */
/*  ========= */

/*  SIDE   - CHARACTER*1 */
/*           SIDE specifies whether A is multiplied on the left or right */
/*           by U. */
/*       SIDE = 'L'   Multiply A on the left (premultiply) by U */
/*       SIDE = 'R'   Multiply A on the right (postmultiply) by U* */
/*       SIDE = 'C'   Multiply A on the left by U and the right by U* */
/*       SIDE = 'T'   Multiply A on the left by U and the right by U' */
/*           Not modified. */

/*  INIT   - CHARACTER*1 */
/*           INIT specifies whether or not A should be initialized to */
/*           the identity matrix. */
/*              INIT = 'I'   Initialize A to (a section of) the */
/*                           identity matrix before applying U. */
/*              INIT = 'N'   No initialization.  Apply U to the */
/*                           input matrix A. */

/*           INIT = 'I' may be used to generate square (i.e., unitary) */
/*           or rectangular orthogonal matrices (orthogonality being */
/*           in the sense of CDOTC): */

/*           For square matrices, M=N, and SIDE many be either 'L' or */
/*           'R'; the rows will be orthogonal to each other, as will the */
/*           columns. */
/*           For rectangular matrices where M < N, SIDE = 'R' will */
/*           produce a dense matrix whose rows will be orthogonal and */
/*           whose columns will not, while SIDE = 'L' will produce a */
/*           matrix whose rows will be orthogonal, and whose first M */
/*           columns will be orthogonal, the remaining columns being */
/*           zero. */
/*           For matrices where M > N, just use the previous */
/*           explaination, interchanging 'L' and 'R' and "rows" and */
/*           "columns". */

/*           Not modified. */

/*  M      - INTEGER */
/*           Number of rows of A. Not modified. */

/*  N      - INTEGER */
/*           Number of columns of A. Not modified. */

/*  A      - COMPLEX array, dimension ( LDA, N ) */
/*           Input and output array. Overwritten by U A ( if SIDE = 'L' ) */
/*           or by A U ( if SIDE = 'R' ) */
/*           or by U A U* ( if SIDE = 'C') */
/*           or by U A U' ( if SIDE = 'T') on exit. */

/*  LDA    - INTEGER */
/*           Leading dimension of A. Must be at least MAX ( 1, M ). */
/*           Not modified. */

/*  ISEED  - INTEGER array, dimension ( 4 ) */
/*           On entry ISEED specifies the seed of the random number */
/*           generator. The array elements should be between 0 and 4095; */
/*           if not they will be reduced mod 4096.  Also, ISEED(4) must */
/*           be odd.  The random number generator uses a linear */
/*           congruential sequence limited to small integers, and so */
/*           should produce machine independent random numbers. The */
/*           values of ISEED are changed on exit, and can be used in the */
/*           next call to CLAROR to continue the same random number */
/*           sequence. */
/*           Modified. */

/*  X      - COMPLEX array, dimension ( 3*MAX( M, N ) ) */
/*           Workspace. Of length: */
/*               2*M + N if SIDE = 'L', */
/*               2*N + M if SIDE = 'R', */
/*               3*N     if SIDE = 'C' or 'T'. */
/*           Modified. */

/*  INFO   - INTEGER */
/*           An error flag.  It is set to: */
/*            0  if no error. */
/*            1  if CLARND returned a bad random number (installation */
/*               problem) */
/*           -1  if SIDE is not L, R, C, or T. */
/*           -3  if M is negative. */
/*           -4  if N is negative or if SIDE is C or T and N is not equal */
/*               to M. */
/*           -6  if LDA is less than M. */

/*  ===================================================================== */

/*     .. Parameters .. */
/*     .. */
/*     .. Local Scalars .. */
/*     .. */
/*     .. External Functions .. */
/*     .. */
/*     .. External Subroutines .. */
/*     .. */
/*     .. Intrinsic Functions .. */
/*     .. */
/*     .. Executable Statements .. */

    /* Parameter adjustments */
    a_dim1 = *lda;
    a_offset = 1 + a_dim1;
    a -= a_offset;
    --iseed;
    --x;

    /* Function Body */
    if (*n == 0 || *m == 0) {
        return 0;
    }

    itype = 0;
    if (lsame_(side, "L")) {
        itype = 1;
    } else if (lsame_(side, "R")) {
        itype = 2;
    } else if (lsame_(side, "C")) {
        itype = 3;
    } else if (lsame_(side, "T")) {
        itype = 4;
    }

/*     Check for argument errors. */

    *info = 0;
    if (itype == 0) {
        *info = -1;
    } else if (*m < 0) {
        *info = -3;
    } else if (*n < 0 || itype == 3 && *n != *m) {
        *info = -4;
    } else if (*lda < *m) {
        *info = -6;
    }
    if (*info != 0) {
        i__1 = -(*info);
        xerbla_("CLAROR", &i__1);
        return 0;
    }

    if (itype == 1) {
        nxfrm = *m;
    } else {
        nxfrm = *n;
    }

/*     Initialize A to the identity matrix if desired */

    if (lsame_(init, "I")) {
        claset_("Full", m, n, &c_b1, &c_b2, &a[a_offset], lda);
    }

/*     If no rotation possible, still multiply by */
/*     a random complex number from the circle |x| = 1 */

/*      2)      Compute Rotation by computing Householder */
/*              Transformations H(2), H(3), ..., H(n).  Note that the */
/*              order in which they are computed is irrelevant. */

    i__1 = nxfrm;
    for (j = 1; j <= i__1; ++j) {
        i__2 = j;
        x[i__2].r = 0.f, x[i__2].i = 0.f;
/* L40: */
    }

    i__1 = nxfrm;
    for (ixfrm = 2; ixfrm <= i__1; ++ixfrm) {
        kbeg = nxfrm - ixfrm + 1;

/*        Generate independent normal( 0, 1 ) random numbers */

        i__2 = nxfrm;
        for (j = kbeg; j <= i__2; ++j) {
            i__3 = j;
            clarnd_(&q__1, &c__3, &iseed[1]);
            x[i__3].r = q__1.r, x[i__3].i = q__1.i;
/* L50: */
        }

/*        Generate a Householder transformation from the random vector X */

        xnorm = scnrm2_(&ixfrm, &x[kbeg], &c__1);
        xabs = c_abs(&x[kbeg]);
        if (xabs != 0.f) {
            i__2 = kbeg;
            q__1.r = x[i__2].r / xabs, q__1.i = x[i__2].i / xabs;
            csign.r = q__1.r, csign.i = q__1.i;
        } else {
            csign.r = 1.f, csign.i = 0.f;
        }
        q__1.r = xnorm * csign.r, q__1.i = xnorm * csign.i;
        xnorms.r = q__1.r, xnorms.i = q__1.i;
        i__2 = nxfrm + kbeg;
        q__1.r = -csign.r, q__1.i = -csign.i;
        x[i__2].r = q__1.r, x[i__2].i = q__1.i;
        factor = xnorm * (xnorm + xabs);
        if (dabs(factor) < 1e-20f) {
            *info = 1;
            i__2 = -(*info);
            xerbla_("CLAROR", &i__2);
            return 0;
        } else {
            factor = 1.f / factor;
        }
        i__2 = kbeg;
        i__3 = kbeg;
        q__1.r = x[i__3].r + xnorms.r, q__1.i = x[i__3].i + xnorms.i;
        x[i__2].r = q__1.r, x[i__2].i = q__1.i;

/*        Apply Householder transformation to A */

        if (itype == 1 || itype == 3 || itype == 4) {

/*           Apply H(k) on the left of A */

            cgemv_("C", &ixfrm, n, &c_b2, &a[kbeg + a_dim1], lda, &x[kbeg], &
                    c__1, &c_b1, &x[(nxfrm << 1) + 1], &c__1);
            q__2.r = factor, q__2.i = 0.f;
            q__1.r = -q__2.r, q__1.i = -q__2.i;
            cgerc_(&ixfrm, n, &q__1, &x[kbeg], &c__1, &x[(nxfrm << 1) + 1], &
                    c__1, &a[kbeg + a_dim1], lda);

        }

        if (itype >= 2 && itype <= 4) {

/*           Apply H(k)* (or H(k)') on the right of A */

            if (itype == 4) {
                clacgv_(&ixfrm, &x[kbeg], &c__1);
            }

            cgemv_("N", m, &ixfrm, &c_b2, &a[kbeg * a_dim1 + 1], lda, &x[kbeg]
, &c__1, &c_b1, &x[(nxfrm << 1) + 1], &c__1);
            q__2.r = factor, q__2.i = 0.f;
            q__1.r = -q__2.r, q__1.i = -q__2.i;
            cgerc_(m, &ixfrm, &q__1, &x[(nxfrm << 1) + 1], &c__1, &x[kbeg], &
                    c__1, &a[kbeg * a_dim1 + 1], lda);

        }
/* L60: */
    }

    clarnd_(&q__1, &c__3, &iseed[1]);
    x[1].r = q__1.r, x[1].i = q__1.i;
    xabs = c_abs(&x[1]);
    if (xabs != 0.f) {
        q__1.r = x[1].r / xabs, q__1.i = x[1].i / xabs;
        csign.r = q__1.r, csign.i = q__1.i;
    } else {
        csign.r = 1.f, csign.i = 0.f;
    }
    i__1 = nxfrm << 1;
    x[i__1].r = csign.r, x[i__1].i = csign.i;

/*     Scale the matrix A by D. */

    if (itype == 1 || itype == 3 || itype == 4) {
        i__1 = *m;
        for (irow = 1; irow <= i__1; ++irow) {
            r_cnjg(&q__1, &x[nxfrm + irow]);
            cscal_(n, &q__1, &a[irow + a_dim1], lda);
/* L70: */
        }
    }

    if (itype == 2 || itype == 3) {
        i__1 = *n;
        for (jcol = 1; jcol <= i__1; ++jcol) {
            cscal_(m, &x[nxfrm + jcol], &a[jcol * a_dim1 + 1], &c__1);
/* L80: */
        }
    }

    if (itype == 4) {
        i__1 = *n;
        for (jcol = 1; jcol <= i__1; ++jcol) {
            r_cnjg(&q__1, &x[nxfrm + jcol]);
            cscal_(m, &q__1, &a[jcol * a_dim1 + 1], &c__1);
/* L90: */
        }
    }
    return 0;

/*     End of CLAROR */

} /* claror_ */