slatms.c

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

	*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.f) {
	*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_("SLATMS", &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 */

    slatm1_(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 (dabs(d[1]) <= (r__1 = d[mnmin], dabs(r__1))) {
	topdwn = TRUE_;
    } else {
	topdwn = FALSE_;
    }

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

/*        Scale by DMAX */

	temp = dabs(d[1]);
	i__1 = mnmin;
	for (i = 2; i <= i__1; ++i) {
/* Computing MAX */
	    r__2 = temp, r__3 = (r__1 = d[i], dabs(r__1));
	    temp = dmax(r__2,r__3);
/* L20: */
	}

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

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

    }

/*     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, SY, GB, 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;
    slaset_("Full", lda, n, &c_b22, &c_b22, &a[a_offset], lda);

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

    if (llb == 0 && uub == 0) {
	i__1 = ilda + 1;
	scopy_(&mnmin, &d[1], &c__1, &a[1 - iskew + ioffst + a_dim1], &i__1);
	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 = ilda + 1;
	    scopy_(&mnmin, &d[1], &c__1, &a[1 - iskew + ioffst + a_dim1], &
		    i__1);

	    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 = 0.f;
			angle = slarnd_(&c__1, &iseed[1]) * 
				6.2831853071795864769252867663f;
			c = cos(angle);
			s = sin(angle);
/* 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;
			    slarot_(&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) {
				slartg_(&a[ir + 1 - iskew * (ic + 1) + ioffst 
					+ (ic + 1) * a_dim1], &extra, &c, &s, 
					&dummy);
			    }
/* Computing MAX */
			    i__4 = 1, i__5 = jch - jku;
			    irow = max(i__4,i__5);
			    il = ir + 2 - irow;
			    temp = 0.f;
			    iltemp = jch > jku;
			    r__1 = -(doublereal)s;
			    slarot_(&c_false, &iltemp, &c_true, &il, &c, &
				    r__1, &a[irow - iskew * ic + ioffst + ic *
				     a_dim1], &ilda, &temp, &extra);
			    if (iltemp) {
				slartg_(&a[irow + 1 - iskew * (ic + 1) + 
					ioffst + (ic + 1) * a_dim1], &temp, &
					c, &s, &dummy);
/* Computing MAX */
				i__4 = 1, i__5 = jch - jku - jkl;
				icol = max(i__4,i__5);
				il = ic + 2 - icol;
				extra = 0.f;
				L__1 = jch > jku + jkl;
				r__1 = -(doublereal)s;
				slarot_(&c_true, &L__1, &c_true, &il, &c, &
					r__1, &a[irow - iskew * icol + ioffst 
					+ icol * a_dim1], &ilda, &extra, &
					temp);
				ic = icol;
				ir = irow;
			    }
/* L30: */
			}
/* L40: */
		    }
/* L50: */
		}

		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 = 0.f;
			angle = slarnd_(&c__1, &iseed[1]) * 
				6.2831853071795864769252867663f;
			c = cos(angle);
			s = sin(angle);
/* 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;
			    slarot_(&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) {
				slartg_(&a[ir + 1 - iskew * (ic + 1) + ioffst 
					+ (ic + 1) * a_dim1], &extra, &c, &s, 
					&dummy);
			    }
/* Computing MAX */
			    i__4 = 1, i__5 = jch - jkl;
			    icol = max(i__4,i__5);
			    il = ic + 2 - icol;
			    temp = 0.f;
			    iltemp = jch > jkl;
			    r__1 = -(doublereal)s;
			    slarot_(&c_true, &iltemp, &c_true, &il, &c, &r__1,
				     &a[ir - iskew * icol + ioffst + icol * 
				    a_dim1], &ilda, &temp, &extra);
			    if (iltemp) {
				slartg_(&a[ir + 1 - iskew * (icol + 1) + 
					ioffst + (icol + 1) * a_dim1], &temp, 
					&c, &s, &dummy);
/* Computing MAX */
				i__4 = 1, i__5 = jch - jkl - jku;
				irow = max(i__4,i__5);
				il = ir + 2 - irow;
				extra = 0.f;
				L__1 = jch > jkl + jku;
				r__1 = -(doublereal)s;
				slarot_(&c_false, &L__1, &c_true, &il, &c, &
					r__1, &a[irow - iskew * icol + ioffst 
					+ icol * a_dim1], &ilda, &extra, &
					temp);
				ic = icol;
				ir = irow;
			    }
/* L60: */
			}
/* L70: */
		    }
/* L80: */
		}

	    } else {

/*              Bottom-Up -- Start at the bottom right. */

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

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

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

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

/*                    Chase "EXTRA" back down */

			ic = jc;
			i__2 = iendch;
			i__4 = jkl + jku;
			for (jch = jc + jkl; i__4 < 0 ? jch >= i__2 : jch <= 
				i__2; jch += i__4) {
			    ilextr = ic > 0;
			    if (ilextr) {
				slartg_(&a[jch - iskew * ic + ioffst + ic * 
					a_dim1], &extra, &c, &s, &dummy);
			    }
			    ic = max(1,ic);
/* Computing MIN */
			    i__5 = *n - 1, i__6 = jch + jku;
			    icol = min(i__5,i__6);
			    iltemp = jch + jku < *n;
			    temp = 0.f;
			    i__5 = icol + 2 - ic;
			    slarot_(&c_true, &ilextr, &iltemp, &i__5, &c, &s, 
				    &a[jch - iskew * ic + ioffst + ic * 
				    a_dim1], &ilda, &extra, &temp);
			    if (iltemp) {
				slartg_(&a[jch - iskew * icol + ioffst + icol 
					* a_dim1], &temp, &c, &s, &dummy);
/* Computing MIN */
				i__5 = iendch, i__6 = jch + jkl + jku;
				il = min(i__5,i__6) + 2 - jch;
				extra = 0.f;
				L__1 = jch + jkl + jku <= iendch;
				slarot_(&c_false, &c_true, &L__1, &il, &c, &s,
					 &a[jch - iskew * icol + ioffst + 
					icol * a_dim1], &ilda, &temp, &extra);
				ic = icol;
			    }
/* L90: */
			}
/* L100: */
		    }
/* L110: */
		}

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

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

                   First row actually rotated is MIN( N+JK
L, M )   
                   First column actually rotated is N   

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

/*                    Chase "EXTRA" back down */

			ir = jr;
			i__3 = iendch;
			i__2 = jkl + jku;
			for (jch = jr + jku; i__2 < 0 ? jch >= i__3 : jch <= 
				i__3; jch += i__2) {
			    ilextr = ir > 0;
			    if (ilextr) {
				slartg_(&a[ir - iskew * jch + ioffst + jch * 
					a_dim1], &extra, &c, &s, &dummy);
			    }
			    ir = max(1,ir);
/* Computing MIN */
			    i__5 = *m - 1, i__6 = jch + jkl;
			    irow = min(i__5,i__6);
			    iltemp = jch + jkl < *m;
			    temp = 0.f;
			    i__5 = irow + 2 - ir;
			    slarot_(&c_false, &ilextr, &iltemp, &i__5, &c, &s,
				     &a[ir - iskew * jch + ioffst + jch * 
				    a_dim1], &ilda, &extra, &temp);
			    if (iltemp) {
				slartg_(&a[irow - iskew * jch + ioffst + jch *
					 a_dim1], &temp, &c, &s, &dummy);
/* Computing MIN */
				i__5 = iendch, i__6 = jch + jkl + jku;