/* dlaror.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 doublereal c_b9 = 0.;
static doublereal c_b10 = 1.;
static integer c__3 = 3;
static integer c__1 = 1;

/* Subroutine */ int dlaror_(char *side, char *init, integer *m, integer *n, 
	doublereal *a, integer *lda, integer *iseed, doublereal *x, integer *
	info)
{
    /* System generated locals */
    integer a_dim1, a_offset, i__1, i__2;
    doublereal d__1;

    /* Builtin functions */
    double d_sign(doublereal *, doublereal *);

    /* Local variables */
    integer j, kbeg;
    extern /* Subroutine */ int dger_(integer *, integer *, doublereal *, 
	    doublereal *, integer *, doublereal *, integer *, doublereal *, 
	    integer *);
    integer jcol, irow;
    extern doublereal dnrm2_(integer *, doublereal *, integer *);
    extern /* Subroutine */ int dscal_(integer *, doublereal *, doublereal *, 
	    integer *);
    extern logical lsame_(char *, char *);
    extern /* Subroutine */ int dgemv_(char *, integer *, integer *, 
	    doublereal *, doublereal *, integer *, doublereal *, integer *, 
	    doublereal *, doublereal *, integer *);
    integer ixfrm, itype, nxfrm;
    doublereal xnorm;
    extern doublereal dlarnd_(integer *, integer *);
    extern /* Subroutine */ int dlaset_(char *, integer *, integer *, 
	    doublereal *, doublereal *, doublereal *, integer *), 
	    xerbla_(char *, integer *);
    doublereal factor, 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 */
/*  ======= */

/*  DLAROR pre- or post-multiplies an M by N matrix A by a random */
/*  orthogonal 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, 403-409). */

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

/*  SIDE    (input) CHARACTER*1 */
/*          Specifies whether A is multiplied on the left or right by U. */
/*          = 'L':         Multiply A on the left (premultiply) by U */
/*          = 'R':         Multiply A on the right (postmultiply) by U' */
/*          = 'C' or 'T':  Multiply A on the left by U and the right */
/*                          by U' (Here, U' means U-transpose.) */

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

/*          INIT = 'I' may be used to generate square or rectangular */
/*          orthogonal matrices: */

/*          For M = N and SIDE = 'L' or 'R', the rows will be orthogonal */
/*          to each other, as will the columns. */

/*          If M < N, SIDE = 'R' produces a dense matrix whose rows are */
/*          orthogonal and whose columns are not, while SIDE = 'L' */
/*          produces a matrix whose rows are orthogonal, and whose first */
/*          M columns are orthogonal, and whose remaining columns are */
/*          zero. */

/*          If M > N, SIDE = 'L' produces a dense matrix whose columns */
/*          are orthogonal and whose rows are not, while SIDE = 'R' */
/*          produces a matrix whose columns are orthogonal, and whose */
/*          first M rows are orthogonal, and whose remaining rows are */
/*          zero. */

/*  M       (input) INTEGER */
/*          The number of rows of A. */

/*  N       (input) INTEGER */
/*          The number of columns of A. */

/*  A       (input/output) DOUBLE PRECISION array, dimension (LDA, N) */
/*          On entry, the array A. */
/*          On exit, overwritten by U A ( if SIDE = 'L' ), */
/*           or by A U ( if SIDE = 'R' ), */
/*           or by U A U' ( if SIDE = 'C' or 'T'). */

/*  LDA     (input) INTEGER */
/*          The leading dimension of the array A.  LDA >= max(1,M). */

/*  ISEED   (input/output) 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 DLAROR to continue the same random number */
/*          sequence. */

/*  X       (workspace) DOUBLE PRECISION 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'. */

/*  INFO    (output) INTEGER */
/*          An error flag.  It is set to: */
/*          = 0:  normal return */
/*          < 0:  if INFO = -k, the k-th argument had an illegal value */
/*          = 1:  if the random numbers generated by DLARND are bad. */

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

/*     .. 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") || lsame_(side, "T")) {
	itype = 3;
    }

/*     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_("DLAROR", &i__1);
	return 0;
    }

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

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

    if (lsame_(init, "I")) {
	dlaset_("Full", m, n, &c_b9, &c_b10, &a[a_offset], lda);
    }

/*     If no rotation possible, multiply by random +/-1 */

/*     Compute rotation by computing Householder transformations */
/*     H(2), H(3), ..., H(nhouse) */

    i__1 = nxfrm;
    for (j = 1; j <= i__1; ++j) {
	x[j] = 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) {
	    x[j] = dlarnd_(&c__3, &iseed[1]);
/* L20: */
	}

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

	xnorm = dnrm2_(&ixfrm, &x[kbeg], &c__1);
	xnorms = d_sign(&xnorm, &x[kbeg]);
	d__1 = -x[kbeg];
	x[kbeg + nxfrm] = d_sign(&c_b10, &d__1);
	factor = xnorms * (xnorms + x[kbeg]);
	if (abs(factor) < 1e-20) {
	    *info = 1;
	    xerbla_("DLAROR", info);
	    return 0;
	} else {
	    factor = 1. / factor;
	}
	x[kbeg] += xnorms;

/*        Apply Householder transformation to A */

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

/*           Apply H(k) from the left. */

	    dgemv_("T", &ixfrm, n, &c_b10, &a[kbeg + a_dim1], lda, &x[kbeg], &
		    c__1, &c_b9, &x[(nxfrm << 1) + 1], &c__1);
	    d__1 = -factor;
	    dger_(&ixfrm, n, &d__1, &x[kbeg], &c__1, &x[(nxfrm << 1) + 1], &
		    c__1, &a[kbeg + a_dim1], lda);

	}

	if (itype == 2 || itype == 3) {

/*           Apply H(k) from the right. */

	    dgemv_("N", m, &ixfrm, &c_b10, &a[kbeg * a_dim1 + 1], lda, &x[
		    kbeg], &c__1, &c_b9, &x[(nxfrm << 1) + 1], &c__1);
	    d__1 = -factor;
	    dger_(m, &ixfrm, &d__1, &x[(nxfrm << 1) + 1], &c__1, &x[kbeg], &
		    c__1, &a[kbeg * a_dim1 + 1], lda);

	}
/* L30: */
    }

    d__1 = dlarnd_(&c__3, &iseed[1]);
    x[nxfrm * 2] = d_sign(&c_b10, &d__1);

/*     Scale the matrix A by D. */

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

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

/*     End of DLAROR */

} /* dlaror_ */