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

/* Table of constant values */

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

/* Subroutine */ int zlaror_(char *side, char *init, integer *m, integer *n, 
	doublecomplex *a, integer *lda, integer *iseed, doublecomplex *x, 
	integer *info)
{
    /* System generated locals */
    integer a_dim1, a_offset, i__1, i__2, i__3;
    doublecomplex z__1, z__2;

    /* Builtin functions */
    double z_abs(doublecomplex *);
    void d_cnjg(doublecomplex *, doublecomplex *);

    /* Local variables */
    integer j, kbeg, jcol;
    doublereal xabs;
    integer irow;
    extern logical lsame_(char *, char *);
    doublecomplex csign;
    extern /* Subroutine */ int zgerc_(integer *, integer *, doublecomplex *, 
	    doublecomplex *, integer *, doublecomplex *, integer *, 
	    doublecomplex *, integer *), zscal_(integer *, doublecomplex *, 
	    doublecomplex *, integer *);
    integer ixfrm;
    extern /* Subroutine */ int zgemv_(char *, integer *, integer *, 
	    doublecomplex *, doublecomplex *, integer *, doublecomplex *, 
	    integer *, doublecomplex *, doublecomplex *, integer *);
    integer itype, nxfrm;
    doublereal xnorm;
    extern doublereal dznrm2_(integer *, doublecomplex *, integer *);
    extern /* Subroutine */ int xerbla_(char *, integer *);
    doublereal factor;
    extern /* Subroutine */ int zlacgv_(integer *, doublecomplex *, integer *)
	    ;
    extern /* Double Complex */ VOID zlarnd_(doublecomplex *, integer *, 
	    integer *);
    extern /* Subroutine */ int zlaset_(char *, integer *, integer *, 
	    doublecomplex *, doublecomplex *, doublecomplex *, integer *);
    doublecomplex 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 */
/*  ======= */

/*     ZLAROR 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 ZDOTC): */

/*           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*16 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 ZLAROR to continue the same random number */
/*           sequence. */
/*           Modified. */

/*  X      - COMPLEX*16 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 ZLARND 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_("ZLAROR", &i__1);
	return 0;
    }

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

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

    if (lsame_(init, "I")) {
	zlaset_("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., x[i__2].i = 0.;
/* L10: */
    }

    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;
	    zlarnd_(&z__1, &c__3, &iseed[1]);
	    x[i__3].r = z__1.r, x[i__3].i = z__1.i;
/* L20: */
	}

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

	xnorm = dznrm2_(&ixfrm, &x[kbeg], &c__1);
	xabs = z_abs(&x[kbeg]);
	if (xabs != 0.) {
	    i__2 = kbeg;
	    z__1.r = x[i__2].r / xabs, z__1.i = x[i__2].i / xabs;
	    csign.r = z__1.r, csign.i = z__1.i;
	} else {
	    csign.r = 1., csign.i = 0.;
	}
	z__1.r = xnorm * csign.r, z__1.i = xnorm * csign.i;
	xnorms.r = z__1.r, xnorms.i = z__1.i;
	i__2 = nxfrm + kbeg;
	z__1.r = -csign.r, z__1.i = -csign.i;
	x[i__2].r = z__1.r, x[i__2].i = z__1.i;
	factor = xnorm * (xnorm + xabs);
	if (abs(factor) < 1e-20) {
	    *info = 1;
	    i__2 = -(*info);
	    xerbla_("ZLAROR", &i__2);
	    return 0;
	} else {
	    factor = 1. / factor;
	}
	i__2 = kbeg;
	i__3 = kbeg;
	z__1.r = x[i__3].r + xnorms.r, z__1.i = x[i__3].i + xnorms.i;
	x[i__2].r = z__1.r, x[i__2].i = z__1.i;

/*        Apply Householder transformation to A */

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

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

	    zgemv_("C", &ixfrm, n, &c_b2, &a[kbeg + a_dim1], lda, &x[kbeg], &
		    c__1, &c_b1, &x[(nxfrm << 1) + 1], &c__1);
	    z__2.r = factor, z__2.i = 0.;
	    z__1.r = -z__2.r, z__1.i = -z__2.i;
	    zgerc_(&ixfrm, n, &z__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) {
		zlacgv_(&ixfrm, &x[kbeg], &c__1);
	    }

	    zgemv_("N", m, &ixfrm, &c_b2, &a[kbeg * a_dim1 + 1], lda, &x[kbeg]
, &c__1, &c_b1, &x[(nxfrm << 1) + 1], &c__1);
	    z__2.r = factor, z__2.i = 0.;
	    z__1.r = -z__2.r, z__1.i = -z__2.i;
	    zgerc_(m, &ixfrm, &z__1, &x[(nxfrm << 1) + 1], &c__1, &x[kbeg], &
		    c__1, &a[kbeg * a_dim1 + 1], lda);

	}
/* L30: */
    }

    zlarnd_(&z__1, &c__3, &iseed[1]);
    x[1].r = z__1.r, x[1].i = z__1.i;
    xabs = z_abs(&x[1]);
    if (xabs != 0.) {
	z__1.r = x[1].r / xabs, z__1.i = x[1].i / xabs;
	csign.r = z__1.r, csign.i = z__1.i;
    } else {
	csign.r = 1., csign.i = 0.;
    }
    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) {
	    d_cnjg(&z__1, &x[nxfrm + irow]);
	    zscal_(n, &z__1, &a[irow + a_dim1], lda);
/* L40: */
	}
    }

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

    if (itype == 4) {
	i__1 = *n;
	for (jcol = 1; jcol <= i__1; ++jcol) {
	    d_cnjg(&z__1, &x[nxfrm + jcol]);
	    zscal_(m, &z__1, &a[jcol * a_dim1 + 1], &c__1);
/* L60: */
	}
    }
    return 0;

/*     End of ZLAROR */

} /* zlaror_ */