kusano 7d535a
/*  -- translated by f2c (version 19940927).
kusano 7d535a
   You must link the resulting object file with the libraries:
kusano 7d535a
	-lf2c -lm   (in that order)
kusano 7d535a
*/
kusano 7d535a
kusano 7d535a
#include "f2c.h"
kusano 7d535a
kusano 7d535a
/* Table of constant values */
kusano 7d535a
kusano 7d535a
static doublereal c_b9 = 0.;
kusano 7d535a
static doublereal c_b10 = 1.;
kusano 7d535a
static integer c__3 = 3;
kusano 7d535a
static integer c__1 = 1;
kusano 7d535a
kusano 7d535a
/* Subroutine */ int dlaror_(char *side, char *init, integer *m, integer *n, 
kusano 7d535a
	doublereal *a, integer *lda, integer *iseed, doublereal *x, integer *
kusano 7d535a
	info)
kusano 7d535a
{
kusano 7d535a
    /* System generated locals */
kusano 7d535a
    integer a_dim1, a_offset, i__1, i__2;
kusano 7d535a
    doublereal d__1;
kusano 7d535a
kusano 7d535a
    /* Builtin functions */
kusano 7d535a
    double d_sign(doublereal *, doublereal *);
kusano 7d535a
kusano 7d535a
    /* Local variables */
kusano 7d535a
    static integer kbeg;
kusano 7d535a
    extern /* Subroutine */ int dger_(integer *, integer *, doublereal *, 
kusano 7d535a
	    doublereal *, integer *, doublereal *, integer *, doublereal *, 
kusano 7d535a
	    integer *);
kusano 7d535a
    static integer jcol, irow;
kusano 7d535a
    extern doublereal dnrm2_(integer *, doublereal *, integer *);
kusano 7d535a
    static integer j;
kusano 7d535a
    extern /* Subroutine */ int dscal_(integer *, doublereal *, doublereal *, 
kusano 7d535a
	    integer *);
kusano 7d535a
    extern logical lsame_(char *, char *);
kusano 7d535a
    extern /* Subroutine */ int dgemv_(char *, integer *, integer *, 
kusano 7d535a
	    doublereal *, doublereal *, integer *, doublereal *, integer *, 
kusano 7d535a
	    doublereal *, doublereal *, integer *);
kusano 7d535a
    static integer ixfrm, itype, nxfrm;
kusano 7d535a
    static doublereal xnorm;
kusano 7d535a
    extern doublereal dlarnd_(integer *, integer *);
kusano 7d535a
    extern /* Subroutine */ int dlaset_(char *, integer *, integer *, 
kusano 7d535a
	    doublereal *, doublereal *, doublereal *, integer *), 
kusano 7d535a
	    xerbla_(char *, integer *);
kusano 7d535a
    static doublereal factor, xnorms;
kusano 7d535a
kusano 7d535a
kusano 7d535a
/*  -- LAPACK auxiliary test routine (version 2.0) --   
kusano 7d535a
       Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,   
kusano 7d535a
       Courant Institute, Argonne National Lab, and Rice University   
kusano 7d535a
       September 30, 1994   
kusano 7d535a
kusano 7d535a
kusano 7d535a
    Purpose   
kusano 7d535a
    =======   
kusano 7d535a
kusano 7d535a
    DLAROR pre- or post-multiplies an M by N matrix A by a random   
kusano 7d535a
    orthogonal matrix U, overwriting A.  A may optionally be initialized 
kusano 7d535a
  
kusano 7d535a
    to the identity matrix before multiplying by U.  U is generated using 
kusano 7d535a
  
kusano 7d535a
    the method of G.W. Stewart (SIAM J. Numer. Anal. 17, 1980, 403-409). 
kusano 7d535a
  
kusano 7d535a
kusano 7d535a
    Arguments   
kusano 7d535a
    =========   
kusano 7d535a
kusano 7d535a
    SIDE    (input) CHARACTER*1   
kusano 7d535a
            Specifies whether A is multiplied on the left or right by U. 
kusano 7d535a
  
kusano 7d535a
            = 'L':         Multiply A on the left (premultiply) by U   
kusano 7d535a
            = 'R':         Multiply A on the right (postmultiply) by U'   
kusano 7d535a
            = 'C' or 'T':  Multiply A on the left by U and the right   
kusano 7d535a
                            by U' (Here, U' means U-transpose.)   
kusano 7d535a
kusano 7d535a
    INIT    (input) CHARACTER*1   
kusano 7d535a
            Specifies whether or not A should be initialized to the   
kusano 7d535a
            identity matrix.   
kusano 7d535a
            = 'I':  Initialize A to (a section of) the identity matrix   
kusano 7d535a
                     before applying U.   
kusano 7d535a
            = 'N':  No initialization.  Apply U to the input matrix A.   
kusano 7d535a
kusano 7d535a
            INIT = 'I' may be used to generate square or rectangular   
kusano 7d535a
            orthogonal matrices:   
kusano 7d535a
kusano 7d535a
            For M = N and SIDE = 'L' or 'R', the rows will be orthogonal 
kusano 7d535a
  
kusano 7d535a
            to each other, as will the columns.   
kusano 7d535a
kusano 7d535a
            If M < N, SIDE = 'R' produces a dense matrix whose rows are   
kusano 7d535a
            orthogonal and whose columns are not, while SIDE = 'L'   
kusano 7d535a
            produces a matrix whose rows are orthogonal, and whose first 
kusano 7d535a
  
kusano 7d535a
            M columns are orthogonal, and whose remaining columns are   
kusano 7d535a
            zero.   
kusano 7d535a
kusano 7d535a
            If M > N, SIDE = 'L' produces a dense matrix whose columns   
kusano 7d535a
            are orthogonal and whose rows are not, while SIDE = 'R'   
kusano 7d535a
            produces a matrix whose columns are orthogonal, and whose   
kusano 7d535a
            first M rows are orthogonal, and whose remaining rows are   
kusano 7d535a
            zero.   
kusano 7d535a
kusano 7d535a
    M       (input) INTEGER   
kusano 7d535a
            The number of rows of A.   
kusano 7d535a
kusano 7d535a
    N       (input) INTEGER   
kusano 7d535a
            The number of columns of A.   
kusano 7d535a
kusano 7d535a
    A       (input/output) DOUBLE PRECISION array, dimension (LDA, N)   
kusano 7d535a
            On entry, the array A.   
kusano 7d535a
            On exit, overwritten by U A ( if SIDE = 'L' ),   
kusano 7d535a
             or by A U ( if SIDE = 'R' ),   
kusano 7d535a
             or by U A U' ( if SIDE = 'C' or 'T').   
kusano 7d535a
kusano 7d535a
    LDA     (input) INTEGER   
kusano 7d535a
            The leading dimension of the array A.  LDA >= max(1,M).   
kusano 7d535a
kusano 7d535a
    ISEED   (input/output) INTEGER array, dimension (4)   
kusano 7d535a
            On entry ISEED specifies the seed of the random number   
kusano 7d535a
            generator. The array elements should be between 0 and 4095;   
kusano 7d535a
            if not they will be reduced mod 4096.  Also, ISEED(4) must   
kusano 7d535a
            be odd.  The random number generator uses a linear   
kusano 7d535a
            congruential sequence limited to small integers, and so   
kusano 7d535a
            should produce machine independent random numbers. The   
kusano 7d535a
            values of ISEED are changed on exit, and can be used in the   
kusano 7d535a
            next call to DLAROR to continue the same random number   
kusano 7d535a
            sequence.   
kusano 7d535a
kusano 7d535a
    X       (workspace) DOUBLE PRECISION array, dimension (3*MAX( M, N )) 
kusano 7d535a
  
kusano 7d535a
            Workspace of length   
kusano 7d535a
                2*M + N if SIDE = 'L',   
kusano 7d535a
                2*N + M if SIDE = 'R',   
kusano 7d535a
                3*N     if SIDE = 'C' or 'T'.   
kusano 7d535a
kusano 7d535a
    INFO    (output) INTEGER   
kusano 7d535a
            An error flag.  It is set to:   
kusano 7d535a
            = 0:  normal return   
kusano 7d535a
            < 0:  if INFO = -k, the k-th argument had an illegal value   
kusano 7d535a
            = 1:  if the random numbers generated by DLARND are bad.   
kusano 7d535a
kusano 7d535a
    ===================================================================== 
kusano 7d535a
  
kusano 7d535a
kusano 7d535a
kusano 7d535a
       Parameter adjustments */
kusano 7d535a
    a_dim1 = *lda;
kusano 7d535a
    a_offset = a_dim1 + 1;
kusano 7d535a
    a -= a_offset;
kusano 7d535a
    --iseed;
kusano 7d535a
    --x;
kusano 7d535a
kusano 7d535a
    /* Function Body */
kusano 7d535a
    if (*n == 0 || *m == 0) {
kusano 7d535a
	return 0;
kusano 7d535a
    }
kusano 7d535a
kusano 7d535a
    itype = 0;
kusano 7d535a
    if (lsame_(side, "L")) {
kusano 7d535a
	itype = 1;
kusano 7d535a
    } else if (lsame_(side, "R")) {
kusano 7d535a
	itype = 2;
kusano 7d535a
    } else if (lsame_(side, "C") || lsame_(side, "T")) {
kusano 7d535a
	itype = 3;
kusano 7d535a
    }
kusano 7d535a
kusano 7d535a
/*     Check for argument errors. */
kusano 7d535a
kusano 7d535a
    *info = 0;
kusano 7d535a
    if (itype == 0) {
kusano 7d535a
	*info = -1;
kusano 7d535a
    } else if (*m < 0) {
kusano 7d535a
	*info = -3;
kusano 7d535a
    } else if (*n < 0 || itype == 3 && *n != *m) {
kusano 7d535a
	*info = -4;
kusano 7d535a
    } else if (*lda < *m) {
kusano 7d535a
	*info = -6;
kusano 7d535a
    }
kusano 7d535a
    if (*info != 0) {
kusano 7d535a
	i__1 = -(*info);
kusano 7d535a
	xerbla_("DLAROR", &i__1);
kusano 7d535a
	return 0;
kusano 7d535a
    }
kusano 7d535a
kusano 7d535a
    if (itype == 1) {
kusano 7d535a
	nxfrm = *m;
kusano 7d535a
    } else {
kusano 7d535a
	nxfrm = *n;
kusano 7d535a
    }
kusano 7d535a
kusano 7d535a
/*     Initialize A to the identity matrix if desired */
kusano 7d535a
kusano 7d535a
    if (lsame_(init, "I")) {
kusano 7d535a
	dlaset_("Full", m, n, &c_b9, &c_b10, &a[a_offset], lda);
kusano 7d535a
    }
kusano 7d535a
kusano 7d535a
/*     If no rotation possible, multiply by random +/-1   
kusano 7d535a
kusano 7d535a
       Compute rotation by computing Householder transformations   
kusano 7d535a
       H(2), H(3), ..., H(nhouse) */
kusano 7d535a
kusano 7d535a
    i__1 = nxfrm;
kusano 7d535a
    for (j = 1; j <= i__1; ++j) {
kusano 7d535a
	x[j] = 0.;
kusano 7d535a
/* L10: */
kusano 7d535a
    }
kusano 7d535a
kusano 7d535a
    i__1 = nxfrm;
kusano 7d535a
    for (ixfrm = 2; ixfrm <= i__1; ++ixfrm) {
kusano 7d535a
	kbeg = nxfrm - ixfrm + 1;
kusano 7d535a
kusano 7d535a
/*        Generate independent normal( 0, 1 ) random numbers */
kusano 7d535a
kusano 7d535a
	i__2 = nxfrm;
kusano 7d535a
	for (j = kbeg; j <= i__2; ++j) {
kusano 7d535a
	    x[j] = dlarnd_(&c__3, &iseed[1]);
kusano 7d535a
/* L20: */
kusano 7d535a
	}
kusano 7d535a
kusano 7d535a
/*        Generate a Householder transformation from the random vector
kusano 7d535a
 X */
kusano 7d535a
kusano 7d535a
	xnorm = dnrm2_(&ixfrm, &x[kbeg], &c__1);
kusano 7d535a
	xnorms = d_sign(&xnorm, &x[kbeg]);
kusano 7d535a
	d__1 = -x[kbeg];
kusano 7d535a
	x[kbeg + nxfrm] = d_sign(&c_b10, &d__1);
kusano 7d535a
	factor = xnorms * (xnorms + x[kbeg]);
kusano 7d535a
	if (abs(factor) < 1e-20) {
kusano 7d535a
	    *info = 1;
kusano 7d535a
	    xerbla_("DLAROR", info);
kusano 7d535a
	    return 0;
kusano 7d535a
	} else {
kusano 7d535a
	    factor = 1. / factor;
kusano 7d535a
	}
kusano 7d535a
	x[kbeg] += xnorms;
kusano 7d535a
kusano 7d535a
/*        Apply Householder transformation to A */
kusano 7d535a
kusano 7d535a
	if (itype == 1 || itype == 3) {
kusano 7d535a
kusano 7d535a
/*           Apply H(k) from the left. */
kusano 7d535a
kusano 7d535a
	    dgemv_("T", &ixfrm, n, &c_b10, &a[kbeg + a_dim1], lda, &x[kbeg], &
kusano 7d535a
		    c__1, &c_b9, &x[(nxfrm << 1) + 1], &c__1);
kusano 7d535a
	    d__1 = -factor;
kusano 7d535a
	    dger_(&ixfrm, n, &d__1, &x[kbeg], &c__1, &x[(nxfrm << 1) + 1], &
kusano 7d535a
		    c__1, &a[kbeg + a_dim1], lda);
kusano 7d535a
kusano 7d535a
	}
kusano 7d535a
kusano 7d535a
	if (itype == 2 || itype == 3) {
kusano 7d535a
kusano 7d535a
/*           Apply H(k) from the right. */
kusano 7d535a
kusano 7d535a
	    dgemv_("N", m, &ixfrm, &c_b10, &a[kbeg * a_dim1 + 1], lda, &x[
kusano 7d535a
		    kbeg], &c__1, &c_b9, &x[(nxfrm << 1) + 1], &c__1);
kusano 7d535a
	    d__1 = -factor;
kusano 7d535a
	    dger_(m, &ixfrm, &d__1, &x[(nxfrm << 1) + 1], &c__1, &x[kbeg], &
kusano 7d535a
		    c__1, &a[kbeg * a_dim1 + 1], lda);
kusano 7d535a
kusano 7d535a
	}
kusano 7d535a
/* L30: */
kusano 7d535a
    }
kusano 7d535a
kusano 7d535a
    d__1 = dlarnd_(&c__3, &iseed[1]);
kusano 7d535a
    x[nxfrm * 2] = d_sign(&c_b10, &d__1);
kusano 7d535a
kusano 7d535a
/*     Scale the matrix A by D. */
kusano 7d535a
kusano 7d535a
    if (itype == 1 || itype == 3) {
kusano 7d535a
	i__1 = *m;
kusano 7d535a
	for (irow = 1; irow <= i__1; ++irow) {
kusano 7d535a
	    dscal_(n, &x[nxfrm + irow], &a[irow + a_dim1], lda);
kusano 7d535a
/* L40: */
kusano 7d535a
	}
kusano 7d535a
    }
kusano 7d535a
kusano 7d535a
    if (itype == 2 || itype == 3) {
kusano 7d535a
	i__1 = *n;
kusano 7d535a
	for (jcol = 1; jcol <= i__1; ++jcol) {
kusano 7d535a
	    dscal_(m, &x[nxfrm + jcol], &a[jcol * a_dim1 + 1], &c__1);
kusano 7d535a
/* L50: */
kusano 7d535a
	}
kusano 7d535a
    }
kusano 7d535a
    return 0;
kusano 7d535a
kusano 7d535a
/*     End of DLAROR */
kusano 7d535a
kusano 7d535a
} /* dlaror_ */
kusano 7d535a