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 integer c__3 = 3;
kusano 7d535a
static integer c__1 = 1;
kusano 7d535a
static real c_b12 = 0.f;
kusano 7d535a
static real c_b19 = -1.f;
kusano 7d535a
static real c_b26 = 1.f;
kusano 7d535a
kusano 7d535a
/* Subroutine */ int slagsy_(integer *n, integer *k, real *d, real *a, 
kusano 7d535a
	integer *lda, integer *iseed, real *work, integer *info)
kusano 7d535a
{
kusano 7d535a
    /* System generated locals */
kusano 7d535a
    integer a_dim1, a_offset, i__1, i__2, i__3;
kusano 7d535a
    real r__1;
kusano 7d535a
kusano 7d535a
    /* Builtin functions */
kusano 7d535a
    double r_sign(real *, real *);
kusano 7d535a
kusano 7d535a
    /* Local variables */
kusano 7d535a
    extern /* Subroutine */ int sger_(integer *, integer *, real *, real *, 
kusano 7d535a
	    integer *, real *, integer *, real *, integer *);
kusano 7d535a
    extern real sdot_(integer *, real *, integer *, real *, integer *), 
kusano 7d535a
	    snrm2_(integer *, real *, integer *);
kusano 7d535a
    static integer i, j;
kusano 7d535a
    extern /* Subroutine */ int ssyr2_(char *, integer *, real *, real *, 
kusano 7d535a
	    integer *, real *, integer *, real *, integer *);
kusano 7d535a
    static real alpha;
kusano 7d535a
    extern /* Subroutine */ int sscal_(integer *, real *, real *, integer *), 
kusano 7d535a
	    sgemv_(char *, integer *, integer *, real *, real *, integer *, 
kusano 7d535a
	    real *, integer *, real *, real *, integer *), saxpy_(
kusano 7d535a
	    integer *, real *, real *, integer *, real *, integer *), ssymv_(
kusano 7d535a
	    char *, integer *, real *, real *, integer *, real *, integer *, 
kusano 7d535a
	    real *, real *, integer *);
kusano 7d535a
    static real wa, wb, wn;
kusano 7d535a
    extern /* Subroutine */ int xerbla_(char *, integer *), slarnv_(
kusano 7d535a
	    integer *, integer *, integer *, real *);
kusano 7d535a
    static real tau;
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
       February 29, 1992   
kusano 7d535a
kusano 7d535a
kusano 7d535a
    Purpose   
kusano 7d535a
    =======   
kusano 7d535a
kusano 7d535a
    SLAGSY generates a real symmetric matrix A, by pre- and post-   
kusano 7d535a
    multiplying a real diagonal matrix D with a random orthogonal matrix: 
kusano 7d535a
  
kusano 7d535a
    A = U*D*U'. The semi-bandwidth may then be reduced to k by additional 
kusano 7d535a
  
kusano 7d535a
    orthogonal transformations.   
kusano 7d535a
kusano 7d535a
    Arguments   
kusano 7d535a
    =========   
kusano 7d535a
kusano 7d535a
    N       (input) INTEGER   
kusano 7d535a
            The order of the matrix A.  N >= 0.   
kusano 7d535a
kusano 7d535a
    K       (input) INTEGER   
kusano 7d535a
            The number of nonzero subdiagonals within the band of A.   
kusano 7d535a
            0 <= K <= N-1.   
kusano 7d535a
kusano 7d535a
    D       (input) REAL array, dimension (N)   
kusano 7d535a
            The diagonal elements of the diagonal matrix D.   
kusano 7d535a
kusano 7d535a
    A       (output) REAL array, dimension (LDA,N)   
kusano 7d535a
            The generated n by n symmetric matrix A (the full matrix is   
kusano 7d535a
            stored).   
kusano 7d535a
kusano 7d535a
    LDA     (input) INTEGER   
kusano 7d535a
            The leading dimension of the array A.  LDA >= N.   
kusano 7d535a
kusano 7d535a
    ISEED   (input/output) INTEGER array, dimension (4)   
kusano 7d535a
            On entry, the seed of the random number generator; the array 
kusano 7d535a
  
kusano 7d535a
            elements must be between 0 and 4095, and ISEED(4) must be   
kusano 7d535a
            odd.   
kusano 7d535a
            On exit, the seed is updated.   
kusano 7d535a
kusano 7d535a
    WORK    (workspace) REAL array, dimension (2*N)   
kusano 7d535a
kusano 7d535a
    INFO    (output) INTEGER   
kusano 7d535a
            = 0: successful exit   
kusano 7d535a
            < 0: if INFO = -i, the i-th argument had an illegal value   
kusano 7d535a
kusano 7d535a
    ===================================================================== 
kusano 7d535a
  
kusano 7d535a
kusano 7d535a
kusano 7d535a
       Test the input arguments   
kusano 7d535a
kusano 7d535a
       Parameter adjustments */
kusano 7d535a
    --d;
kusano 7d535a
    a_dim1 = *lda;
kusano 7d535a
    a_offset = a_dim1 + 1;
kusano 7d535a
    a -= a_offset;
kusano 7d535a
    --iseed;
kusano 7d535a
    --work;
kusano 7d535a
kusano 7d535a
    /* Function Body */
kusano 7d535a
    *info = 0;
kusano 7d535a
    if (*n < 0) {
kusano 7d535a
	*info = -1;
kusano 7d535a
    } else if (*k < 0 || *k > *n - 1) {
kusano 7d535a
	*info = -2;
kusano 7d535a
    } else if (*lda < max(1,*n)) {
kusano 7d535a
	*info = -5;
kusano 7d535a
    }
kusano 7d535a
    if (*info < 0) {
kusano 7d535a
	i__1 = -(*info);
kusano 7d535a
	xerbla_("SLAGSY", &i__1);
kusano 7d535a
	return 0;
kusano 7d535a
    }
kusano 7d535a
kusano 7d535a
/*     initialize lower triangle of A to diagonal matrix */
kusano 7d535a
kusano 7d535a
    i__1 = *n;
kusano 7d535a
    for (j = 1; j <= i__1; ++j) {
kusano 7d535a
	i__2 = *n;
kusano 7d535a
	for (i = j + 1; i <= i__2; ++i) {
kusano 7d535a
	    a[i + j * a_dim1] = 0.f;
kusano 7d535a
/* L10: */
kusano 7d535a
	}
kusano 7d535a
/* L20: */
kusano 7d535a
    }
kusano 7d535a
    i__1 = *n;
kusano 7d535a
    for (i = 1; i <= i__1; ++i) {
kusano 7d535a
	a[i + i * a_dim1] = d[i];
kusano 7d535a
/* L30: */
kusano 7d535a
    }
kusano 7d535a
kusano 7d535a
/*     Generate lower triangle of symmetric matrix */
kusano 7d535a
kusano 7d535a
    for (i = *n - 1; i >= 1; --i) {
kusano 7d535a
kusano 7d535a
/*        generate random reflection */
kusano 7d535a
kusano 7d535a
	i__1 = *n - i + 1;
kusano 7d535a
	slarnv_(&c__3, &iseed[1], &i__1, &work[1]);
kusano 7d535a
	i__1 = *n - i + 1;
kusano 7d535a
	wn = snrm2_(&i__1, &work[1], &c__1);
kusano 7d535a
	wa = r_sign(&wn, &work[1]);
kusano 7d535a
	if (wn == 0.f) {
kusano 7d535a
	    tau = 0.f;
kusano 7d535a
	} else {
kusano 7d535a
	    wb = work[1] + wa;
kusano 7d535a
	    i__1 = *n - i;
kusano 7d535a
	    r__1 = 1.f / wb;
kusano 7d535a
	    sscal_(&i__1, &r__1, &work[2], &c__1);
kusano 7d535a
	    work[1] = 1.f;
kusano 7d535a
	    tau = wb / wa;
kusano 7d535a
	}
kusano 7d535a
kusano 7d535a
/*        apply random reflection to A(i:n,i:n) from the left   
kusano 7d535a
          and the right   
kusano 7d535a
kusano 7d535a
          compute  y := tau * A * u */
kusano 7d535a
kusano 7d535a
	i__1 = *n - i + 1;
kusano 7d535a
	ssymv_("Lower", &i__1, &tau, &a[i + i * a_dim1], lda, &work[1], &c__1,
kusano 7d535a
		 &c_b12, &work[*n + 1], &c__1);
kusano 7d535a
kusano 7d535a
/*        compute  v := y - 1/2 * tau * ( y, u ) * u */
kusano 7d535a
kusano 7d535a
	i__1 = *n - i + 1;
kusano 7d535a
	alpha = tau * -.5f * sdot_(&i__1, &work[*n + 1], &c__1, &work[1], &
kusano 7d535a
		c__1);
kusano 7d535a
	i__1 = *n - i + 1;
kusano 7d535a
	saxpy_(&i__1, &alpha, &work[1], &c__1, &work[*n + 1], &c__1);
kusano 7d535a
kusano 7d535a
/*        apply the transformation as a rank-2 update to A(i:n,i:n) */
kusano 7d535a
kusano 7d535a
	i__1 = *n - i + 1;
kusano 7d535a
	ssyr2_("Lower", &i__1, &c_b19, &work[1], &c__1, &work[*n + 1], &c__1, 
kusano 7d535a
		&a[i + i * a_dim1], lda);
kusano 7d535a
/* L40: */
kusano 7d535a
    }
kusano 7d535a
kusano 7d535a
/*     Reduce number of subdiagonals to K */
kusano 7d535a
kusano 7d535a
    i__1 = *n - 1 - *k;
kusano 7d535a
    for (i = 1; i <= i__1; ++i) {
kusano 7d535a
kusano 7d535a
/*        generate reflection to annihilate A(k+i+1:n,i) */
kusano 7d535a
kusano 7d535a
	i__2 = *n - *k - i + 1;
kusano 7d535a
	wn = snrm2_(&i__2, &a[*k + i + i * a_dim1], &c__1);
kusano 7d535a
	wa = r_sign(&wn, &a[*k + i + i * a_dim1]);
kusano 7d535a
	if (wn == 0.f) {
kusano 7d535a
	    tau = 0.f;
kusano 7d535a
	} else {
kusano 7d535a
	    wb = a[*k + i + i * a_dim1] + wa;
kusano 7d535a
	    i__2 = *n - *k - i;
kusano 7d535a
	    r__1 = 1.f / wb;
kusano 7d535a
	    sscal_(&i__2, &r__1, &a[*k + i + 1 + i * a_dim1], &c__1);
kusano 7d535a
	    a[*k + i + i * a_dim1] = 1.f;
kusano 7d535a
	    tau = wb / wa;
kusano 7d535a
	}
kusano 7d535a
kusano 7d535a
/*        apply reflection to A(k+i:n,i+1:k+i-1) from the left */
kusano 7d535a
kusano 7d535a
	i__2 = *n - *k - i + 1;
kusano 7d535a
	i__3 = *k - 1;
kusano 7d535a
	sgemv_("Transpose", &i__2, &i__3, &c_b26, &a[*k + i + (i + 1) * 
kusano 7d535a
		a_dim1], lda, &a[*k + i + i * a_dim1], &c__1, &c_b12, &work[1]
kusano 7d535a
		, &c__1);
kusano 7d535a
	i__2 = *n - *k - i + 1;
kusano 7d535a
	i__3 = *k - 1;
kusano 7d535a
	r__1 = -(doublereal)tau;
kusano 7d535a
	sger_(&i__2, &i__3, &r__1, &a[*k + i + i * a_dim1], &c__1, &work[1], &
kusano 7d535a
		c__1, &a[*k + i + (i + 1) * a_dim1], lda);
kusano 7d535a
kusano 7d535a
/*        apply reflection to A(k+i:n,k+i:n) from the left and the rig
kusano 7d535a
ht   
kusano 7d535a
kusano 7d535a
          compute  y := tau * A * u */
kusano 7d535a
kusano 7d535a
	i__2 = *n - *k - i + 1;
kusano 7d535a
	ssymv_("Lower", &i__2, &tau, &a[*k + i + (*k + i) * a_dim1], lda, &a[*
kusano 7d535a
		k + i + i * a_dim1], &c__1, &c_b12, &work[1], &c__1);
kusano 7d535a
kusano 7d535a
/*        compute  v := y - 1/2 * tau * ( y, u ) * u */
kusano 7d535a
kusano 7d535a
	i__2 = *n - *k - i + 1;
kusano 7d535a
	alpha = tau * -.5f * sdot_(&i__2, &work[1], &c__1, &a[*k + i + i * 
kusano 7d535a
		a_dim1], &c__1);
kusano 7d535a
	i__2 = *n - *k - i + 1;
kusano 7d535a
	saxpy_(&i__2, &alpha, &a[*k + i + i * a_dim1], &c__1, &work[1], &c__1)
kusano 7d535a
		;
kusano 7d535a
kusano 7d535a
/*        apply symmetric rank-2 update to A(k+i:n,k+i:n) */
kusano 7d535a
kusano 7d535a
	i__2 = *n - *k - i + 1;
kusano 7d535a
	ssyr2_("Lower", &i__2, &c_b19, &a[*k + i + i * a_dim1], &c__1, &work[
kusano 7d535a
		1], &c__1, &a[*k + i + (*k + i) * a_dim1], lda);
kusano 7d535a
kusano 7d535a
	a[*k + i + i * a_dim1] = -(doublereal)wa;
kusano 7d535a
	i__2 = *n;
kusano 7d535a
	for (j = *k + i + 1; j <= i__2; ++j) {
kusano 7d535a
	    a[j + i * a_dim1] = 0.f;
kusano 7d535a
/* L50: */
kusano 7d535a
	}
kusano 7d535a
/* L60: */
kusano 7d535a
    }
kusano 7d535a
kusano 7d535a
/*     Store full symmetric matrix */
kusano 7d535a
kusano 7d535a
    i__1 = *n;
kusano 7d535a
    for (j = 1; j <= i__1; ++j) {
kusano 7d535a
	i__2 = *n;
kusano 7d535a
	for (i = j + 1; i <= i__2; ++i) {
kusano 7d535a
	    a[j + i * a_dim1] = a[i + j * a_dim1];
kusano 7d535a
/* L70: */
kusano 7d535a
	}
kusano 7d535a
/* L80: */
kusano 7d535a
    }
kusano 7d535a
    return 0;
kusano 7d535a
kusano 7d535a
/*     End of SLAGSY */
kusano 7d535a
kusano 7d535a
} /* slagsy_ */
kusano 7d535a