zlatms.c

来自「LU矩阵分解单机版最新版本」· C语言 代码 · 共 1,649 行 · 第 1/4 页

    } else if (lsame_(pack, "L")) {
	ipack = 2;
	isympk = 1;
    } else if (lsame_(pack, "C")) {
	ipack = 3;
	isympk = 2;
    } else if (lsame_(pack, "R")) {
	ipack = 4;
	isympk = 3;
    } else if (lsame_(pack, "B")) {
	ipack = 5;
	isympk = 3;
    } else if (lsame_(pack, "Q")) {
	ipack = 6;
	isympk = 2;
    } else if (lsame_(pack, "Z")) {
	ipack = 7;
    } else {
	ipack = -1;
    }

/*     Set certain internal parameters */

    mnmin = min(*m,*n);
/* Computing MIN */
    i__1 = *kl, i__2 = *m - 1;
    llb = min(i__1,i__2);
/* Computing MIN */
    i__1 = *ku, i__2 = *n - 1;
    uub = min(i__1,i__2);
/* Computing MIN */
    i__1 = *m, i__2 = *n + llb;
    mr = min(i__1,i__2);
/* Computing MIN */
    i__1 = *n, i__2 = *m + uub;
    nc = min(i__1,i__2);

    if (ipack == 5 || ipack == 6) {
	minlda = uub + 1;
    } else if (ipack == 7) {
	minlda = llb + uub + 1;
    } else {
	minlda = *m;
    }

/*     Use Givens rotation method if bandwidth small enough,   
       or if LDA is too small to store the matrix unpacked. */

    givens = FALSE_;
    if (isym == 1) {
/* Computing MAX */
	i__1 = 1, i__2 = mr + nc;
	if ((doublereal) (llb + uub) < (doublereal) max(i__1,i__2) * .3) {
	    givens = TRUE_;
	}
    } else {
	if (llb << 1 < *m) {
	    givens = TRUE_;
	}
    }
    if (*lda < *m && *lda >= minlda) {
	givens = TRUE_;
    }

/*     Set INFO if an error */

    if (*m < 0) {
	*info = -1;
    } else if (*m != *n && isym != 1) {
	*info = -1;
    } else if (*n < 0) {
	*info = -2;
    } else if (idist == -1) {
	*info = -3;
    } else if (isym == -1) {
	*info = -5;
    } else if (abs(*mode) > 6) {
	*info = -7;
    } else if (*mode != 0 && abs(*mode) != 6 && *cond < 1.) {
	*info = -8;
    } else if (*kl < 0) {
	*info = -10;
    } else if (*ku < 0 || isym != 1 && *kl != *ku) {
	*info = -11;
    } else if (ipack == -1 || isympk == 1 && isym == 1 || isympk == 2 && isym 
	    == 1 && *kl > 0 || isympk == 3 && isym == 1 && *ku > 0 || isympk 
	    != 0 && *m != *n) {
	*info = -12;
    } else if (*lda < max(1,minlda)) {
	*info = -14;
    }

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

/*     Initialize random number generator */

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

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

/*     2)      Set up D  if indicated.   

               Compute D according to COND and MODE */

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

/*     Choose Top-Down if D is (apparently) increasing,   
       Bottom-Up if D is (apparently) decreasing. */

    if (abs(d[1]) <= (d__1 = d[mnmin], abs(d__1))) {
	topdwn = TRUE_;
    } else {
	topdwn = FALSE_;
    }

    if (*mode != 0 && abs(*mode) != 6) {

/*        Scale by DMAX */

	temp = abs(d[1]);
	i__1 = mnmin;
	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);
/* L20: */
	}

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

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

    }

    zlaset_("Full", lda, n, &c_b1, &c_b1, &a[a_offset], lda);

/*     3)      Generate Banded Matrix using Givens rotations.   
               Also the special case of UUB=LLB=0   

                 Compute Addressing constants to cover all   
                 storage formats.  Whether GE, HE, SY, GB, HB, or SB,   
                 upper or lower triangle or both,   
                 the (i,j)-th element is in   
                 A( i - ISKEW*j + IOFFST, j ) */

    if (ipack > 4) {
	ilda = *lda - 1;
	iskew = 1;
	if (ipack > 5) {
	    ioffst = uub + 1;
	} else {
	    ioffst = 1;
	}
    } else {
	ilda = *lda;
	iskew = 0;
	ioffst = 0;
    }

/*     IPACKG is the format that the matrix is generated in. If this is   
       different from IPACK, then the matrix must be repacked at the   
       end.  It also signals how to compute the norm, for scaling. */

    ipackg = 0;

/*     Diagonal Matrix -- We are done, unless it   
       is to be stored HP/SP/PP/TP (PACK='R' or 'C') */

    if (llb == 0 && uub == 0) {
	i__1 = mnmin;
	for (j = 1; j <= i__1; ++j) {
	    i__2 = (1 - iskew) * j + ioffst + j * a_dim1;
	    i__3 = j;
	    z__1.r = d[i__3], z__1.i = 0.;
	    a[i__2].r = z__1.r, a[i__2].i = z__1.i;
/* L30: */
	}

	if (ipack <= 2 || ipack >= 5) {
	    ipackg = ipack;
	}

    } else if (givens) {

/*        Check whether to use Givens rotations,   
          Householder transformations, or nothing. */

	if (isym == 1) {

/*           Non-symmetric -- A = U D V */

	    if (ipack > 4) {
		ipackg = ipack;
	    } else {
		ipackg = 0;
	    }

	    i__1 = mnmin;
	    for (j = 1; j <= i__1; ++j) {
		i__2 = (1 - iskew) * j + ioffst + j * a_dim1;
		i__3 = j;
		z__1.r = d[i__3], z__1.i = 0.;
		a[i__2].r = z__1.r, a[i__2].i = z__1.i;
/* L40: */
	    }

	    if (topdwn) {
		jkl = 0;
		i__1 = uub;
		for (jku = 1; jku <= i__1; ++jku) {

/*                 Transform from bandwidth JKL, JKU-1 to 
JKL, JKU   

                   Last row actually rotated is M   
                   Last column actually rotated is MIN( M+
JKU, N )   

   Computing MIN */
		    i__3 = *m + jku;
		    i__2 = min(i__3,*n) + jkl - 1;
		    for (jr = 1; jr <= i__2; ++jr) {
			extra.r = 0., extra.i = 0.;
			angle = dlarnd_(&c__1, &iseed[1]) * 
				6.2831853071795864769252867663;
			d__1 = cos(angle);
			zlarnd_(&z__2, &c__5, &iseed[1]);
			z__1.r = d__1 * z__2.r, z__1.i = d__1 * z__2.i;
			c.r = z__1.r, c.i = z__1.i;
			d__1 = sin(angle);
			zlarnd_(&z__2, &c__5, &iseed[1]);
			z__1.r = d__1 * z__2.r, z__1.i = d__1 * z__2.i;
			s.r = z__1.r, s.i = z__1.i;
/* Computing MAX */
			i__3 = 1, i__4 = jr - jkl;
			icol = max(i__3,i__4);
			if (jr < *m) {
/* Computing MIN */
			    i__3 = *n, i__4 = jr + jku;
			    il = min(i__3,i__4) + 1 - icol;
			    L__1 = jr > jkl;
			    zlarot_(&c_true, &L__1, &c_false, &il, &c, &s, &a[
				    jr - iskew * icol + ioffst + icol * 
				    a_dim1], &ilda, &extra, &dummy);
			}

/*                    Chase "EXTRA" back up */

			ir = jr;
			ic = icol;
			i__3 = -jkl - jku;
			for (jch = jr - jkl; i__3 < 0 ? jch >= 1 : jch <= 1; 
				jch += i__3) {
			    if (ir < *m) {
				zlartg_(&a[ir + 1 - iskew * (ic + 1) + ioffst 
					+ (ic + 1) * a_dim1], &extra, &realc, 
					&s, &dummy);
				zlarnd_(&z__1, &c__5, &iseed[1]);
				dummy.r = z__1.r, dummy.i = z__1.i;
				z__2.r = realc * dummy.r, z__2.i = realc * 
					dummy.i;
				d_cnjg(&z__1, &z__2);
				c.r = z__1.r, c.i = z__1.i;
				z__3.r = -s.r, z__3.i = -s.i;
				z__2.r = z__3.r * dummy.r - z__3.i * dummy.i, 
					z__2.i = z__3.r * dummy.i + z__3.i * 
					dummy.r;
				d_cnjg(&z__1, &z__2);
				s.r = z__1.r, s.i = z__1.i;
			    }
/* Computing MAX */
			    i__4 = 1, i__5 = jch - jku;
			    irow = max(i__4,i__5);
			    il = ir + 2 - irow;
			    ctemp.r = 0., ctemp.i = 0.;
			    iltemp = jch > jku;
			    zlarot_(&c_false, &iltemp, &c_true, &il, &c, &s, &
				    a[irow - iskew * ic + ioffst + ic * 
				    a_dim1], &ilda, &ctemp, &extra);
			    if (iltemp) {
				zlartg_(&a[irow + 1 - iskew * (ic + 1) + 
					ioffst + (ic + 1) * a_dim1], &ctemp, &
					realc, &s, &dummy);
				zlarnd_(&z__1, &c__5, &iseed[1]);
				dummy.r = z__1.r, dummy.i = z__1.i;
				z__2.r = realc * dummy.r, z__2.i = realc * 
					dummy.i;
				d_cnjg(&z__1, &z__2);
				c.r = z__1.r, c.i = z__1.i;
				z__3.r = -s.r, z__3.i = -s.i;
				z__2.r = z__3.r * dummy.r - z__3.i * dummy.i, 
					z__2.i = z__3.r * dummy.i + z__3.i * 
					dummy.r;
				d_cnjg(&z__1, &z__2);
				s.r = z__1.r, s.i = z__1.i;

/* Computing MAX */
				i__4 = 1, i__5 = jch - jku - jkl;
				icol = max(i__4,i__5);
				il = ic + 2 - icol;
				extra.r = 0., extra.i = 0.;
				L__1 = jch > jku + jkl;
				zlarot_(&c_true, &L__1, &c_true, &il, &c, &s, 
					&a[irow - iskew * icol + ioffst + 
					icol * a_dim1], &ilda, &extra, &ctemp)
					;
				ic = icol;
				ir = irow;
			    }
/* L50: */
			}
/* L60: */
		    }
/* L70: */
		}

		jku = uub;
		i__1 = llb;
		for (jkl = 1; jkl <= i__1; ++jkl) {

/*                 Transform from bandwidth JKL-1, JKU to 
JKL, JKU   

   Computing MIN */
		    i__3 = *n + jkl;
		    i__2 = min(i__3,*m) + jku - 1;
		    for (jc = 1; jc <= i__2; ++jc) {
			extra.r = 0., extra.i = 0.;
			angle = dlarnd_(&c__1, &iseed[1]) * 
				6.2831853071795864769252867663;
			d__1 = cos(angle);
			zlarnd_(&z__2, &c__5, &iseed[1]);
			z__1.r = d__1 * z__2.r, z__1.i = d__1 * z__2.i;
			c.r = z__1.r, c.i = z__1.i;
			d__1 = sin(angle);
			zlarnd_(&z__2, &c__5, &iseed[1]);
			z__1.r = d__1 * z__2.r, z__1.i = d__1 * z__2.i;
			s.r = z__1.r, s.i = z__1.i;
/* Computing MAX */
			i__3 = 1, i__4 = jc - jku;
			irow = max(i__3,i__4);
			if (jc < *n) {
/* Computing MIN */
			    i__3 = *m, i__4 = jc + jkl;
			    il = min(i__3,i__4) + 1 - irow;
			    L__1 = jc > jku;
			    zlarot_(&c_false, &L__1, &c_false, &il, &c, &s, &
				    a[irow - iskew * jc + ioffst + jc * 
				    a_dim1], &ilda, &extra, &dummy);
			}

/*                    Chase "EXTRA" back up */

			ic = jc;
			ir = irow;
			i__3 = -jkl - jku;
			for (jch = jc - jku; i__3 < 0 ? jch >= 1 : jch <= 1; 
				jch += i__3) {
			    if (ic < *n) {
				zlartg_(&a[ir + 1 - iskew * (ic + 1) + ioffst 
					+ (ic + 1) * a_dim1], &extra, &realc, 
					&s, &dummy);
				zlarnd_(&z__1, &c__5, &iseed[1]);
				dummy.r = z__1.r, dummy.i = z__1.i;
				z__2.r = realc * dummy.r, z__2.i = realc * 
					dummy.i;
				d_cnjg(&z__1, &z__2);
				c.r = z__1.r, c.i = z__1.i;
				z__3.r = -s.r, z__3.i = -s.i;
				z__2.r = z__3.r * dummy.r - z__3.i * dummy.i, 
					z__2.i = z__3.r * dummy.i + z__3.i * 
					dummy.r;
				d_cnjg(&z__1, &z__2);
				s.r = z__1.r, s.i = z__1.i;
			    }
/* Computing MAX */
			    i__4 = 1, i__5 = jch - jkl;
			    icol = max(i__4,i__5);
			    il = ic + 2 - icol;
			    ctemp.r = 0., ctemp.i = 0.;
			    iltemp = jch > jkl;
			    zlarot_(&c_true, &iltemp, &c_true, &il, &c, &s, &
				    a[ir - iskew * icol + ioffst + icol * 
				    a_dim1], &ilda, &ctemp, &extra);
			    if (iltemp) {
				zlartg_(&a[ir + 1 - iskew * (icol + 1) + 
					ioffst + (icol + 1) * a_dim1], &ctemp,
					 &realc, &s, &dummy);
				zlarnd_(&z__1, &c__5, &iseed[1]);
				dummy.r = z__1.r, dummy.i = z__1.i;
				z__2.r = realc * dummy.r, z__2.i = realc * 
					dummy.i;
				d_cnjg(&z__1, &z__2);
				c.r = z__1.r, c.i = z__1.i;
				z__3.r = -s.r, z__3.i = -s.i;