#include "blaswrap.h"
/* 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"

/* 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 */
    static integer j, kbeg;
    extern /* Subroutine */ int dger_(integer *, integer *, doublereal *, 
	    doublereal *, integer *, doublereal *, integer *, doublereal *, 
	    integer *);
    static 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 *);
    static integer ixfrm, itype, nxfrm;
    static doublereal xnorm;
    extern doublereal dlarnd_(integer *, integer *);
    extern /* Subroutine */ int dlaset_(char *, integer *, integer *, 
	    doublereal *, doublereal *, doublereal *, integer *), 
	    xerbla_(char *, integer *);
    static doublereal factor, xnorms;


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


    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.   

    =====================================================================   


       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_ */