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 complex c_b1 = {0.f,0.f};
kusano 7d535a
static complex c_b2 = {1.f,0.f};
kusano 7d535a
static integer c__3 = 3;
kusano 7d535a
static integer c__1 = 1;
kusano 7d535a
kusano 7d535a
/* Subroutine */ int clagsy_(integer *n, integer *k, real *d, complex *a, 
kusano 7d535a
	integer *lda, integer *iseed, complex *work, integer *info)
kusano 7d535a
{
kusano 7d535a
    /* System generated locals */
kusano 7d535a
    integer a_dim1, a_offset, i__1, i__2, i__3, i__4, i__5, i__6, i__7, i__8, 
kusano 7d535a
	    i__9;
kusano 7d535a
    doublereal d__1;
kusano 7d535a
    complex q__1, q__2, q__3, q__4;
kusano 7d535a
kusano 7d535a
    /* Builtin functions */
kusano 7d535a
    double c_abs(complex *);
kusano 7d535a
    void c_div(complex *, complex *, complex *);
kusano 7d535a
kusano 7d535a
    /* Local variables */
kusano 7d535a
    static integer i, j;
kusano 7d535a
    extern /* Subroutine */ int cgerc_(integer *, integer *, complex *, 
kusano 7d535a
	    complex *, integer *, complex *, integer *, complex *, integer *);
kusano 7d535a
    static complex alpha;
kusano 7d535a
    extern /* Subroutine */ int cscal_(integer *, complex *, complex *, 
kusano 7d535a
	    integer *);
kusano 7d535a
    extern /* Complex */ VOID cdotc_(complex *, integer *, complex *, integer 
kusano 7d535a
	    *, complex *, integer *);
kusano 7d535a
    extern /* Subroutine */ int cgemv_(char *, integer *, integer *, complex *
kusano 7d535a
	    , complex *, integer *, complex *, integer *, complex *, complex *
kusano 7d535a
	    , integer *), caxpy_(integer *, complex *, complex *, 
kusano 7d535a
	    integer *, complex *, integer *), csymv_(char *, integer *, 
kusano 7d535a
	    complex *, complex *, integer *, complex *, integer *, complex *, 
kusano 7d535a
	    complex *, integer *);
kusano 7d535a
    extern real scnrm2_(integer *, complex *, integer *);
kusano 7d535a
    static integer ii, jj;
kusano 7d535a
    static complex wa, wb;
kusano 7d535a
    extern /* Subroutine */ int clacgv_(integer *, complex *, integer *);
kusano 7d535a
    static real wn;
kusano 7d535a
    extern /* Subroutine */ int xerbla_(char *, integer *), clarnv_(
kusano 7d535a
	    integer *, integer *, integer *, complex *);
kusano 7d535a
    static complex 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
       September 30, 1994   
kusano 7d535a
kusano 7d535a
kusano 7d535a
    Purpose   
kusano 7d535a
    =======   
kusano 7d535a
kusano 7d535a
    CLAGSY generates a complex symmetric matrix A, by pre- and post-   
kusano 7d535a
    multiplying a real diagonal matrix D with a random unitary matrix:   
kusano 7d535a
    A = U*D*U**T. The semi-bandwidth may then be reduced to k by   
kusano 7d535a
    additional unitary 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) COMPLEX 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) COMPLEX 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_("CLAGSY", &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
	    i__3 = i + j * a_dim1;
kusano 7d535a
	    a[i__3].r = 0.f, a[i__3].i = 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
	i__2 = i + i * a_dim1;
kusano 7d535a
	i__3 = i;
kusano 7d535a
	a[i__2].r = d[i__3], a[i__2].i = 0.f;
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
	clarnv_(&c__3, &iseed[1], &i__1, &work[1]);
kusano 7d535a
	i__1 = *n - i + 1;
kusano 7d535a
	wn = scnrm2_(&i__1, &work[1], &c__1);
kusano 7d535a
	d__1 = wn / c_abs(&work[1]);
kusano 7d535a
	q__1.r = d__1 * work[1].r, q__1.i = d__1 * work[1].i;
kusano 7d535a
	wa.r = q__1.r, wa.i = q__1.i;
kusano 7d535a
	if (wn == 0.f) {
kusano 7d535a
	    tau.r = 0.f, tau.i = 0.f;
kusano 7d535a
	} else {
kusano 7d535a
	    q__1.r = work[1].r + wa.r, q__1.i = work[1].i + wa.i;
kusano 7d535a
	    wb.r = q__1.r, wb.i = q__1.i;
kusano 7d535a
	    i__1 = *n - i;
kusano 7d535a
	    c_div(&q__1, &c_b2, &wb);
kusano 7d535a
	    cscal_(&i__1, &q__1, &work[2], &c__1);
kusano 7d535a
	    work[1].r = 1.f, work[1].i = 0.f;
kusano 7d535a
	    c_div(&q__1, &wb, &wa);
kusano 7d535a
	    d__1 = q__1.r;
kusano 7d535a
	    tau.r = d__1, tau.i = 0.f;
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 * conjg(u) */
kusano 7d535a
kusano 7d535a
	i__1 = *n - i + 1;
kusano 7d535a
	clacgv_(&i__1, &work[1], &c__1);
kusano 7d535a
	i__1 = *n - i + 1;
kusano 7d535a
	csymv_("Lower", &i__1, &tau, &a[i + i * a_dim1], lda, &work[1], &c__1,
kusano 7d535a
		 &c_b1, &work[*n + 1], &c__1);
kusano 7d535a
	i__1 = *n - i + 1;
kusano 7d535a
	clacgv_(&i__1, &work[1], &c__1);
kusano 7d535a
kusano 7d535a
/*        compute  v := y - 1/2 * tau * ( u, y ) * u */
kusano 7d535a
kusano 7d535a
	q__3.r = -.5f, q__3.i = 0.f;
kusano 7d535a
	q__2.r = q__3.r * tau.r - q__3.i * tau.i, q__2.i = q__3.r * tau.i + 
kusano 7d535a
		q__3.i * tau.r;
kusano 7d535a
	i__1 = *n - i + 1;
kusano 7d535a
	cdotc_(&q__4, &i__1, &work[1], &c__1, &work[*n + 1], &c__1);
kusano 7d535a
	q__1.r = q__2.r * q__4.r - q__2.i * q__4.i, q__1.i = q__2.r * q__4.i 
kusano 7d535a
		+ q__2.i * q__4.r;
kusano 7d535a
	alpha.r = q__1.r, alpha.i = q__1.i;
kusano 7d535a
	i__1 = *n - i + 1;
kusano 7d535a
	caxpy_(&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
          CALL CSYR2( 'Lower', N-I+1, -ONE, WORK, 1, WORK( N+1 ), 1, 
kusano 7d535a
  
kusano 7d535a
          $               A( I, I ), LDA ) */
kusano 7d535a
kusano 7d535a
	i__1 = *n;
kusano 7d535a
	for (jj = i; jj <= i__1; ++jj) {
kusano 7d535a
	    i__2 = *n;
kusano 7d535a
	    for (ii = jj; ii <= i__2; ++ii) {
kusano 7d535a
		i__3 = ii + jj * a_dim1;
kusano 7d535a
		i__4 = ii + jj * a_dim1;
kusano 7d535a
		i__5 = ii - i + 1;
kusano 7d535a
		i__6 = *n + jj - i + 1;
kusano 7d535a
		q__3.r = work[i__5].r * work[i__6].r - work[i__5].i * work[
kusano 7d535a
			i__6].i, q__3.i = work[i__5].r * work[i__6].i + work[
kusano 7d535a
			i__5].i * work[i__6].r;
kusano 7d535a
		q__2.r = a[i__4].r - q__3.r, q__2.i = a[i__4].i - q__3.i;
kusano 7d535a
		i__7 = *n + ii - i + 1;
kusano 7d535a
		i__8 = jj - i + 1;
kusano 7d535a
		q__4.r = work[i__7].r * work[i__8].r - work[i__7].i * work[
kusano 7d535a
			i__8].i, q__4.i = work[i__7].r * work[i__8].i + work[
kusano 7d535a
			i__7].i * work[i__8].r;
kusano 7d535a
		q__1.r = q__2.r - q__4.r, q__1.i = q__2.i - q__4.i;
kusano 7d535a
		a[i__3].r = q__1.r, a[i__3].i = q__1.i;
kusano 7d535a
/* L40: */
kusano 7d535a
	    }
kusano 7d535a
/* L50: */
kusano 7d535a
	}
kusano 7d535a
/* L60: */
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 = scnrm2_(&i__2, &a[*k + i + i * a_dim1], &c__1);
kusano 7d535a
	d__1 = wn / c_abs(&a[*k + i + i * a_dim1]);
kusano 7d535a
	i__2 = *k + i + i * a_dim1;
kusano 7d535a
	q__1.r = d__1 * a[i__2].r, q__1.i = d__1 * a[i__2].i;
kusano 7d535a
	wa.r = q__1.r, wa.i = q__1.i;
kusano 7d535a
	if (wn == 0.f) {
kusano 7d535a
	    tau.r = 0.f, tau.i = 0.f;
kusano 7d535a
	} else {
kusano 7d535a
	    i__2 = *k + i + i * a_dim1;
kusano 7d535a
	    q__1.r = a[i__2].r + wa.r, q__1.i = a[i__2].i + wa.i;
kusano 7d535a
	    wb.r = q__1.r, wb.i = q__1.i;
kusano 7d535a
	    i__2 = *n - *k - i;
kusano 7d535a
	    c_div(&q__1, &c_b2, &wb);
kusano 7d535a
	    cscal_(&i__2, &q__1, &a[*k + i + 1 + i * a_dim1], &c__1);
kusano 7d535a
	    i__2 = *k + i + i * a_dim1;
kusano 7d535a
	    a[i__2].r = 1.f, a[i__2].i = 0.f;
kusano 7d535a
	    c_div(&q__1, &wb, &wa);
kusano 7d535a
	    d__1 = q__1.r;
kusano 7d535a
	    tau.r = d__1, tau.i = 0.f;
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
	cgemv_("Conjugate transpose", &i__2, &i__3, &c_b2, &a[*k + i + (i + 1)
kusano 7d535a
		 * a_dim1], lda, &a[*k + i + i * a_dim1], &c__1, &c_b1, &work[
kusano 7d535a
		1], &c__1);
kusano 7d535a
	i__2 = *n - *k - i + 1;
kusano 7d535a
	i__3 = *k - 1;
kusano 7d535a
	q__1.r = -(doublereal)tau.r, q__1.i = -(doublereal)tau.i;
kusano 7d535a
	cgerc_(&i__2, &i__3, &q__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 * conjg(u) */
kusano 7d535a
kusano 7d535a
	i__2 = *n - *k - i + 1;
kusano 7d535a
	clacgv_(&i__2, &a[*k + i + i * a_dim1], &c__1);
kusano 7d535a
	i__2 = *n - *k - i + 1;
kusano 7d535a
	csymv_("Lower", &i__2, &tau, &a[*k + i + (*k + i) * a_dim1], lda, &a[*
kusano 7d535a
		k + i + i * a_dim1], &c__1, &c_b1, &work[1], &c__1);
kusano 7d535a
	i__2 = *n - *k - i + 1;
kusano 7d535a
	clacgv_(&i__2, &a[*k + i + i * a_dim1], &c__1);
kusano 7d535a
kusano 7d535a
/*        compute  v := y - 1/2 * tau * ( u, y ) * u */
kusano 7d535a
kusano 7d535a
	q__3.r = -.5f, q__3.i = 0.f;
kusano 7d535a
	q__2.r = q__3.r * tau.r - q__3.i * tau.i, q__2.i = q__3.r * tau.i + 
kusano 7d535a
		q__3.i * tau.r;
kusano 7d535a
	i__2 = *n - *k - i + 1;
kusano 7d535a
	cdotc_(&q__4, &i__2, &a[*k + i + i * a_dim1], &c__1, &work[1], &c__1);
kusano 7d535a
	q__1.r = q__2.r * q__4.r - q__2.i * q__4.i, q__1.i = q__2.r * q__4.i 
kusano 7d535a
		+ q__2.i * q__4.r;
kusano 7d535a
	alpha.r = q__1.r, alpha.i = q__1.i;
kusano 7d535a
	i__2 = *n - *k - i + 1;
kusano 7d535a
	caxpy_(&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
          CALL CSYR2( 'Lower', N-K-I+1, -ONE, A( K+I, I ), 1, WORK, 1,
kusano 7d535a
   
kusano 7d535a
          $               A( K+I, K+I ), LDA ) */
kusano 7d535a
kusano 7d535a
	i__2 = *n;
kusano 7d535a
	for (jj = *k + i; jj <= i__2; ++jj) {
kusano 7d535a
	    i__3 = *n;
kusano 7d535a
	    for (ii = jj; ii <= i__3; ++ii) {
kusano 7d535a
		i__4 = ii + jj * a_dim1;
kusano 7d535a
		i__5 = ii + jj * a_dim1;
kusano 7d535a
		i__6 = ii + i * a_dim1;
kusano 7d535a
		i__7 = jj - *k - i + 1;
kusano 7d535a
		q__3.r = a[i__6].r * work[i__7].r - a[i__6].i * work[i__7].i, 
kusano 7d535a
			q__3.i = a[i__6].r * work[i__7].i + a[i__6].i * work[
kusano 7d535a
			i__7].r;
kusano 7d535a
		q__2.r = a[i__5].r - q__3.r, q__2.i = a[i__5].i - q__3.i;
kusano 7d535a
		i__8 = ii - *k - i + 1;
kusano 7d535a
		i__9 = jj + i * a_dim1;
kusano 7d535a
		q__4.r = work[i__8].r * a[i__9].r - work[i__8].i * a[i__9].i, 
kusano 7d535a
			q__4.i = work[i__8].r * a[i__9].i + work[i__8].i * a[
kusano 7d535a
			i__9].r;
kusano 7d535a
		q__1.r = q__2.r - q__4.r, q__1.i = q__2.i - q__4.i;
kusano 7d535a
		a[i__4].r = q__1.r, a[i__4].i = q__1.i;
kusano 7d535a
/* L70: */
kusano 7d535a
	    }
kusano 7d535a
/* L80: */
kusano 7d535a
	}
kusano 7d535a
kusano 7d535a
	i__2 = *k + i + i * a_dim1;
kusano 7d535a
	q__1.r = -(doublereal)wa.r, q__1.i = -(doublereal)wa.i;
kusano 7d535a
	a[i__2].r = q__1.r, a[i__2].i = q__1.i;
kusano 7d535a
	i__2 = *n;
kusano 7d535a
	for (j = *k + i + 1; j <= i__2; ++j) {
kusano 7d535a
	    i__3 = j + i * a_dim1;
kusano 7d535a
	    a[i__3].r = 0.f, a[i__3].i = 0.f;
kusano 7d535a
/* L90: */
kusano 7d535a
	}
kusano 7d535a
/* L100: */
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
	    i__3 = j + i * a_dim1;
kusano 7d535a
	    i__4 = i + j * a_dim1;
kusano 7d535a
	    a[i__3].r = a[i__4].r, a[i__3].i = a[i__4].i;
kusano 7d535a
/* L110: */
kusano 7d535a
	}
kusano 7d535a
/* L120: */
kusano 7d535a
    }
kusano 7d535a
    return 0;
kusano 7d535a
kusano 7d535a
/*     End of CLAGSY */
kusano 7d535a
kusano 7d535a
} /* clagsy_ */
kusano 7d535a