dlatme.c

SuperLU is a general purpose library for the direct solution of large, sparse, nonsymmetric systems

    badei = FALSE_;
    if (lsame_(ei + 1, " ") || *mode != 0) {
	useei = FALSE_;
    } else {
	if (lsame_(ei + 1, "R")) {
	    i__1 = *n;
	    for (j = 2; j <= i__1; ++j) {
		if (lsame_(ei + j, "I")) {
		    if (lsame_(ei + (j - 1), "I")) {
			badei = TRUE_;
		    }
		} else {
		    if (! lsame_(ei + j, "R")) {
			badei = TRUE_;
		    }
		}
/* L10: */
	    }
	} else {
	    badei = TRUE_;
	}
    }

/*     Decode RSIGN */

    if (lsame_(rsign, "T")) {
	irsign = 1;
    } else if (lsame_(rsign, "F")) {
	irsign = 0;
    } else {
	irsign = -1;
    }

/*     Decode UPPER */

    if (lsame_(upper, "T")) {
	iupper = 1;
    } else if (lsame_(upper, "F")) {
	iupper = 0;
    } else {
	iupper = -1;
    }

/*     Decode SIM */

    if (lsame_(sim, "T")) {
	isim = 1;
    } else if (lsame_(sim, "F")) {
	isim = 0;
    } else {
	isim = -1;
    }

/*     Check DS, if MODES=0 and ISIM=1 */

    bads = FALSE_;
    if (*modes == 0 && isim == 1) {
	i__1 = *n;
	for (j = 1; j <= i__1; ++j) {
	    if (ds[j] == 0.) {
		bads = TRUE_;
	    }
/* L20: */
	}
    }

/*     Set INFO if an error */

    if (*n < 0) {
	*info = -1;
    } else if (idist == -1) {
	*info = -2;
    } else if (abs(*mode) > 6) {
	*info = -5;
    } else if (*mode != 0 && abs(*mode) != 6 && *cond < 1.) {
	*info = -6;
    } else if (badei) {
	*info = -8;
    } else if (irsign == -1) {
	*info = -9;
    } else if (iupper == -1) {
	*info = -10;
    } else if (isim == -1) {
	*info = -11;
    } else if (bads) {
	*info = -12;
    } else if (isim == 1 && abs(*modes) > 5) {
	*info = -13;
    } else if (isim == 1 && *modes != 0 && *conds < 1.) {
	*info = -14;
    } else if (*kl < 1) {
	*info = -15;
    } else if (*ku < 1 || *ku < *n - 1 && *kl < *n - 1) {
	*info = -16;
    } else if (*lda < max(1,*n)) {
	*info = -19;
    }

    if (*info != 0) {
	i__1 = -(*info);
	xerbla_("DLATME", &i__1);
	return 0;
    }

/*     Initialize random number generator */

    for (i = 1; i <= 4; ++i) {
	iseed[i] = (i__1 = iseed[i], abs(i__1)) % 4096;
/* L30: */
    }

    if (iseed[4] % 2 != 1) {
	++iseed[4];
    }

/*     2)      Set up diagonal of A   

               Compute D according to COND and MODE */

    dlatm1_(mode, cond, &irsign, &idist, &iseed[1], &d[1], n, &iinfo);
    if (iinfo != 0) {
	*info = 1;
	return 0;
    }
    if (*mode != 0 && abs(*mode) != 6) {

/*        Scale by DMAX */

	temp = abs(d[1]);
	i__1 = *n;
	for (i = 2; i <= i__1; ++i) {
/* Computing MAX */
	    d__2 = temp, d__3 = (d__1 = d[i], abs(d__1));
	    temp = max(d__2,d__3);
/* L40: */
	}

	if (temp > 0.) {
	    alpha = *dmax__ / temp;
	} else if (*dmax__ != 0.) {
	    *info = 2;
	    return 0;
	} else {
	    alpha = 0.;
	}

	dscal_(n, &alpha, &d[1], &c__1);

    }

    dlaset_("Full", n, n, &c_b23, &c_b23, &a[a_offset], lda);
    i__1 = *lda + 1;
    dcopy_(n, &d[1], &c__1, &a[a_offset], &i__1);

/*     Set up complex conjugate pairs */

    if (*mode == 0) {
	if (useei) {
	    i__1 = *n;
	    for (j = 2; j <= i__1; ++j) {
		if (lsame_(ei + j, "I")) {
		    a[j - 1 + j * a_dim1] = a[j + j * a_dim1];
		    a[j + (j - 1) * a_dim1] = -a[j + j * a_dim1];
		    a[j + j * a_dim1] = a[j - 1 + (j - 1) * a_dim1];
		}
/* L50: */
	    }
	}

    } else if (abs(*mode) == 5) {

	i__1 = *n;
	for (j = 2; j <= i__1; j += 2) {
	    if (dlaran_(&iseed[1]) > .5) {
		a[j - 1 + j * a_dim1] = a[j + j * a_dim1];
		a[j + (j - 1) * a_dim1] = -a[j + j * a_dim1];
		a[j + j * a_dim1] = a[j - 1 + (j - 1) * a_dim1];
	    }
/* L60: */
	}
    }

/*     3)      If UPPER='T', set upper triangle of A to random numbers.   
               (but don't modify the corners of 2x2 blocks.) */

    if (iupper != 0) {
	i__1 = *n;
	for (jc = 2; jc <= i__1; ++jc) {
	    if (a[jc - 1 + jc * a_dim1] != 0.) {
		jr = jc - 2;
	    } else {
		jr = jc - 1;
	    }
	    dlarnv_(&idist, &iseed[1], &jr, &a[jc * a_dim1 + 1]);
/* L70: */
	}
    }

/*     4)      If SIM='T', apply similarity transformation.   

                                  -1   
               Transform is  X A X  , where X = U S V, thus   

               it is  U S V A V' (1/S) U' */

    if (isim != 0) {

/*        Compute S (singular values of the eigenvector matrix)   
          according to CONDS and MODES */

	dlatm1_(modes, conds, &c__0, &c__0, &iseed[1], &ds[1], n, &iinfo);
	if (iinfo != 0) {
	    *info = 3;
	    return 0;
	}

/*        Multiply by V and V' */

	dlarge_(n, &a[a_offset], lda, &iseed[1], &work[1], &iinfo);
	if (iinfo != 0) {
	    *info = 4;
	    return 0;
	}

/*        Multiply by S and (1/S) */

	i__1 = *n;
	for (j = 1; j <= i__1; ++j) {
	    dscal_(n, &ds[j], &a[j + a_dim1], lda);
	    if (ds[j] != 0.) {
		d__1 = 1. / ds[j];
		dscal_(n, &d__1, &a[j * a_dim1 + 1], &c__1);
	    } else {
		*info = 5;
		return 0;
	    }
/* L80: */
	}

/*        Multiply by U and U' */

	dlarge_(n, &a[a_offset], lda, &iseed[1], &work[1], &iinfo);
	if (iinfo != 0) {
	    *info = 4;
	    return 0;
	}
    }

/*     5)      Reduce the bandwidth. */

    if (*kl < *n - 1) {

/*        Reduce bandwidth -- kill column */

	i__1 = *n - 1;
	for (jcr = *kl + 1; jcr <= i__1; ++jcr) {
	    ic = jcr - *kl;
	    irows = *n + 1 - jcr;
	    icols = *n + *kl - jcr;

	    dcopy_(&irows, &a[jcr + ic * a_dim1], &c__1, &work[1], &c__1);
	    xnorms = work[1];
	    dlarfg_(&irows, &xnorms, &work[2], &c__1, &tau);
	    work[1] = 1.;

	    dgemv_("T", &irows, &icols, &c_b39, &a[jcr + (ic + 1) * a_dim1], 
		    lda, &work[1], &c__1, &c_b23, &work[irows + 1], &c__1)
		    ;
	    d__1 = -tau;
	    dger_(&irows, &icols, &d__1, &work[1], &c__1, &work[irows + 1], &
		    c__1, &a[jcr + (ic + 1) * a_dim1], lda);

	    dgemv_("N", n, &irows, &c_b39, &a[jcr * a_dim1 + 1], lda, &work[1]
		    , &c__1, &c_b23, &work[irows + 1], &c__1);
	    d__1 = -tau;
	    dger_(n, &irows, &d__1, &work[irows + 1], &c__1, &work[1], &c__1, 
		    &a[jcr * a_dim1 + 1], lda);

	    a[jcr + ic * a_dim1] = xnorms;
	    i__2 = irows - 1;
	    dlaset_("Full", &i__2, &c__1, &c_b23, &c_b23, &a[jcr + 1 + ic * 
		    a_dim1], lda);
/* L90: */
	}
    } else if (*ku < *n - 1) {

/*        Reduce upper bandwidth -- kill a row at a time. */

	i__1 = *n - 1;
	for (jcr = *ku + 1; jcr <= i__1; ++jcr) {
	    ir = jcr - *ku;
	    irows = *n + *ku - jcr;
	    icols = *n + 1 - jcr;

	    dcopy_(&icols, &a[ir + jcr * a_dim1], lda, &work[1], &c__1);
	    xnorms = work[1];
	    dlarfg_(&icols, &xnorms, &work[2], &c__1, &tau);
	    work[1] = 1.;

	    dgemv_("N", &irows, &icols, &c_b39, &a[ir + 1 + jcr * a_dim1], 
		    lda, &work[1], &c__1, &c_b23, &work[icols + 1], &c__1)
		    ;
	    d__1 = -tau;
	    dger_(&irows, &icols, &d__1, &work[icols + 1], &c__1, &work[1], &
		    c__1, &a[ir + 1 + jcr * a_dim1], lda);

	    dgemv_("C", &icols, n, &c_b39, &a[jcr + a_dim1], lda, &work[1], &
		    c__1, &c_b23, &work[icols + 1], &c__1);
	    d__1 = -tau;
	    dger_(&icols, n, &d__1, &work[1], &c__1, &work[icols + 1], &c__1, 
		    &a[jcr + a_dim1], lda);

	    a[ir + jcr * a_dim1] = xnorms;
	    i__2 = icols - 1;
	    dlaset_("Full", &c__1, &i__2, &c_b23, &c_b23, &a[ir + (jcr + 1) * 
		    a_dim1], lda);
/* L100: */
	}
    }

/*     Scale the matrix to have norm ANORM */

    if (*anorm >= 0.) {
	temp = dlange_("M", n, n, &a[a_offset], lda, tempa);
	if (temp > 0.) {
	    alpha = *anorm / temp;
	    i__1 = *n;
	    for (j = 1; j <= i__1; ++j) {
		dscal_(n, &alpha, &a[j * a_dim1 + 1], &c__1);
/* L110: */
	    }
	}
    }

    return 0;

/*     End of DLATME */

} /* dlatme_ */