slatms.c

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

				il = min(i__5,i__6) + 2 - jch;
				extra = 0.f;
				L__1 = jch + jkl + jku <= iendch;
				slarot_(&c_true, &c_true, &L__1, &il, &c, &s, 
					&a[irow - iskew * jch + ioffst + jch *
					 a_dim1], &ilda, &temp, &extra);
				ir = irow;
			    }
/* L120: */
			}
/* L130: */
		    }
/* L140: */
		}
	    }

	} else {

/*           Symmetric -- A = U D U' */

	    ipackg = ipack;
	    ioffg = ioffst;

	    if (topdwn) {

/*              Top-Down -- Generate Upper triangle only */

		if (ipack >= 5) {
		    ipackg = 6;
		    ioffg = uub + 1;
		} else {
		    ipackg = 1;
		}
		i__1 = ilda + 1;
		scopy_(&mnmin, &d[1], &c__1, &a[1 - iskew + ioffg + a_dim1], &
			i__1);

		i__1 = uub;
		for (k = 1; k <= i__1; ++k) {
		    i__4 = *n - 1;
		    for (jc = 1; jc <= i__4; ++jc) {
/* Computing MAX */
			i__2 = 1, i__3 = jc - k;
			irow = max(i__2,i__3);
/* Computing MIN */
			i__2 = jc + 1, i__3 = k + 2;
			il = min(i__2,i__3);
			extra = 0.f;
			temp = a[jc - iskew * (jc + 1) + ioffg + (jc + 1) * 
				a_dim1];
			angle = slarnd_(&c__1, &iseed[1]) * 
				6.2831853071795864769252867663f;
			c = cos(angle);
			s = sin(angle);
			L__1 = jc > k;
			slarot_(&c_false, &L__1, &c_true, &il, &c, &s, &a[
				irow - iskew * jc + ioffg + jc * a_dim1], &
				ilda, &extra, &temp);
/* Computing MIN */
			i__3 = k, i__5 = *n - jc;
			i__2 = min(i__3,i__5) + 1;
			slarot_(&c_true, &c_true, &c_false, &i__2, &c, &s, &a[
				(1 - iskew) * jc + ioffg + jc * a_dim1], &
				ilda, &temp, &dummy);

/*                    Chase EXTRA back up the matrix 
*/

			icol = jc;
			i__2 = -k;
			for (jch = jc - k; i__2 < 0 ? jch >= 1 : jch <= 1; 
				jch += i__2) {
			    slartg_(&a[jch + 1 - iskew * (icol + 1) + ioffg + 
				    (icol + 1) * a_dim1], &extra, &c, &s, &
				    dummy);
			    temp = a[jch - iskew * (jch + 1) + ioffg + (jch + 
				    1) * a_dim1];
			    i__3 = k + 2;
			    r__1 = -(doublereal)s;
			    slarot_(&c_true, &c_true, &c_true, &i__3, &c, &
				    r__1, &a[(1 - iskew) * jch + ioffg + jch *
				     a_dim1], &ilda, &temp, &extra);
/* Computing MAX */
			    i__3 = 1, i__5 = jch - k;
			    irow = max(i__3,i__5);
/* Computing MIN */
			    i__3 = jch + 1, i__5 = k + 2;
			    il = min(i__3,i__5);
			    extra = 0.f;
			    L__1 = jch > k;
			    r__1 = -(doublereal)s;
			    slarot_(&c_false, &L__1, &c_true, &il, &c, &r__1, 
				    &a[irow - iskew * jch + ioffg + jch * 
				    a_dim1], &ilda, &extra, &temp);
			    icol = jch;
/* L150: */
			}
/* L160: */
		    }
/* L170: */
		}

/*              If we need lower triangle, copy from upper. No
te that   
                the order of copying is chosen to work for 'q'
 -> 'b' */

		if (ipack != ipackg && ipack != 3) {
		    i__1 = *n;
		    for (jc = 1; jc <= i__1; ++jc) {
			irow = ioffst - iskew * jc;
/* Computing MIN */
			i__2 = *n, i__3 = jc + uub;
			i__4 = min(i__2,i__3);
			for (jr = jc; jr <= i__4; ++jr) {
			    a[jr + irow + jc * a_dim1] = a[jc - iskew * jr + 
				    ioffg + jr * a_dim1];
/* L180: */
			}
/* L190: */
		    }
		    if (ipack == 5) {
			i__1 = *n;
			for (jc = *n - uub + 1; jc <= i__1; ++jc) {
			    i__4 = uub + 1;
			    for (jr = *n + 2 - jc; jr <= i__4; ++jr) {
				a[jr + jc * a_dim1] = 0.f;
/* L200: */
			    }
/* L210: */
			}
		    }
		    if (ipackg == 6) {
			ipackg = ipack;
		    } else {
			ipackg = 0;
		    }
		}
	    } else {

/*              Bottom-Up -- Generate Lower triangle only */

		if (ipack >= 5) {
		    ipackg = 5;
		    if (ipack == 6) {
			ioffg = 1;
		    }
		} else {
		    ipackg = 2;
		}
		i__1 = ilda + 1;
		scopy_(&mnmin, &d[1], &c__1, &a[1 - iskew + ioffg + a_dim1], &
			i__1);

		i__1 = uub;
		for (k = 1; k <= i__1; ++k) {
		    for (jc = *n - 1; jc >= 1; --jc) {
/* Computing MIN */
			i__4 = *n + 1 - jc, i__2 = k + 2;
			il = min(i__4,i__2);
			extra = 0.f;
			temp = a[(1 - iskew) * jc + 1 + ioffg + jc * a_dim1];
			angle = slarnd_(&c__1, &iseed[1]) * 
				6.2831853071795864769252867663f;
			c = cos(angle);
			s = -(doublereal)sin(angle);
			L__1 = *n - jc > k;
			slarot_(&c_false, &c_true, &L__1, &il, &c, &s, &a[(1 
				- iskew) * jc + ioffg + jc * a_dim1], &ilda, &
				temp, &extra);
/* Computing MAX */
			i__4 = 1, i__2 = jc - k + 1;
			icol = max(i__4,i__2);
			i__4 = jc + 2 - icol;
			slarot_(&c_true, &c_false, &c_true, &i__4, &c, &s, &a[
				jc - iskew * icol + ioffg + icol * a_dim1], &
				ilda, &dummy, &temp);

/*                    Chase EXTRA back down the matrix
 */

			icol = jc;
			i__4 = *n - 1;
			i__2 = k;
			for (jch = jc + k; i__2 < 0 ? jch >= i__4 : jch <= 
				i__4; jch += i__2) {
			    slartg_(&a[jch - iskew * icol + ioffg + icol * 
				    a_dim1], &extra, &c, &s, &dummy);
			    temp = a[(1 - iskew) * jch + 1 + ioffg + jch * 
				    a_dim1];
			    i__3 = k + 2;
			    slarot_(&c_true, &c_true, &c_true, &i__3, &c, &s, 
				    &a[jch - iskew * icol + ioffg + icol * 
				    a_dim1], &ilda, &extra, &temp);
/* Computing MIN */
			    i__3 = *n + 1 - jch, i__5 = k + 2;
			    il = min(i__3,i__5);
			    extra = 0.f;
			    L__1 = *n - jch > k;
			    slarot_(&c_false, &c_true, &L__1, &il, &c, &s, &a[
				    (1 - iskew) * jch + ioffg + jch * a_dim1],
				     &ilda, &temp, &extra);
			    icol = jch;
/* L220: */
			}
/* L230: */
		    }
/* L240: */
		}

/*              If we need upper triangle, copy from lower. No
te that   
                the order of copying is chosen to work for 'b'
 -> 'q' */

		if (ipack != ipackg && ipack != 4) {
		    for (jc = *n; jc >= 1; --jc) {
			irow = ioffst - iskew * jc;
/* Computing MAX */
			i__2 = 1, i__4 = jc - uub;
			i__1 = max(i__2,i__4);
			for (jr = jc; jr >= i__1; --jr) {
			    a[jr + irow + jc * a_dim1] = a[jc - iskew * jr + 
				    ioffg + jr * a_dim1];
/* L250: */
			}
/* L260: */
		    }
		    if (ipack == 6) {
			i__1 = uub;
			for (jc = 1; jc <= i__1; ++jc) {
			    i__2 = uub + 1 - jc;
			    for (jr = 1; jr <= i__2; ++jr) {
				a[jr + jc * a_dim1] = 0.f;
/* L270: */
			    }
/* L280: */
			}
		    }
		    if (ipackg == 5) {
			ipackg = ipack;
		    } else {
			ipackg = 0;
		    }
		}
	    }
	}

    } else {

/*        4)      Generate Banded Matrix by first   
                  Rotating by random Unitary matrices,   
                  then reducing the bandwidth using Householder   
                  transformations.   

                  Note: we should get here only if LDA .ge. N */

	if (isym == 1) {

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

	    slagge_(&mr, &nc, &llb, &uub, &d[1], &a[a_offset], lda, &iseed[1],
		     &work[1], &iinfo);
	} else {

/*           Symmetric -- A = U D U' */

	    slagsy_(m, &llb, &d[1], &a[a_offset], lda, &iseed[1], &work[1], &
		    iinfo);

	}
	if (iinfo != 0) {
	    *info = 3;
	    return 0;
	}
    }

/*     5)      Pack the matrix */

    if (ipack != ipackg) {
	if (ipack == 1) {

/*           'U' -- Upper triangular, not packed */

	    i__1 = *m;
	    for (j = 1; j <= i__1; ++j) {
		i__2 = *m;
		for (i = j + 1; i <= i__2; ++i) {
		    a[i + j * a_dim1] = 0.f;
/* L290: */
		}
/* L300: */
	    }

	} else if (ipack == 2) {

/*           'L' -- Lower triangular, not packed */

	    i__1 = *m;
	    for (j = 2; j <= i__1; ++j) {
		i__2 = j - 1;
		for (i = 1; i <= i__2; ++i) {
		    a[i + j * a_dim1] = 0.f;
/* L310: */
		}
/* L320: */
	    }

	} else if (ipack == 3) {

/*           'C' -- Upper triangle packed Columnwise. */

	    icol = 1;
	    irow = 0;
	    i__1 = *m;
	    for (j = 1; j <= i__1; ++j) {
		i__2 = j;
		for (i = 1; i <= i__2; ++i) {
		    ++irow;
		    if (irow > *lda) {
			irow = 1;
			++icol;
		    }
		    a[irow + icol * a_dim1] = a[i + j * a_dim1];
/* L330: */
		}
/* L340: */
	    }

	} else if (ipack == 4) {

/*           'R' -- Lower triangle packed Columnwise. */

	    icol = 1;
	    irow = 0;
	    i__1 = *m;
	    for (j = 1; j <= i__1; ++j) {
		i__2 = *m;
		for (i = j; i <= i__2; ++i) {
		    ++irow;
		    if (irow > *lda) {
			irow = 1;
			++icol;
		    }
		    a[irow + icol * a_dim1] = a[i + j * a_dim1];
/* L350: */
		}
/* L360: */
	    }

	} else if (ipack >= 5) {

/*           'B' -- The lower triangle is packed as a band matrix.
   
             'Q' -- The upper triangle is packed as a band matrix.
   
             'Z' -- The whole matrix is packed as a band matrix. 
*/

	    if (ipack == 5) {
		uub = 0;
	    }
	    if (ipack == 6) {
		llb = 0;
	    }

	    i__1 = uub;
	    for (j = 1; j <= i__1; ++j) {
/* Computing MIN */
		i__2 = j + llb;
		for (i = min(i__2,*m); i >= 1; --i) {
		    a[i - j + uub + 1 + j * a_dim1] = a[i + j * a_dim1];
/* L370: */
		}
/* L380: */
	    }

	    i__1 = *n;
	    for (j = uub + 2; j <= i__1; ++j) {
/* Computing MIN */
		i__4 = j + llb;
		i__2 = min(i__4,*m);
		for (i = j - uub; i <= i__2; ++i) {
		    a[i - j + uub + 1 + j * a_dim1] = a[i + j * a_dim1];
/* L390: */
		}
/* L400: */
	    }
	}

/*        If packed, zero out extraneous elements.   

          Symmetric/Triangular Packed --   
          zero out everything after A(IROW,ICOL) */

	if (ipack == 3 || ipack == 4) {
	    i__1 = *m;
	    for (jc = icol; jc <= i__1; ++jc) {
		i__2 = *lda;
		for (jr = irow + 1; jr <= i__2; ++jr) {
		    a[jr + jc * a_dim1] = 0.f;
/* L410: */
		}
		irow = 0;
/* L420: */
	    }

	} else if (ipack >= 5) {

/*           Packed Band --   
                1st row is now in A( UUB+2-j, j), zero above it   
                m-th row is now in A( M+UUB-j,j), zero below it   
                last non-zero diagonal is now in A( UUB+LLB+1,j ),
   
                   zero below it, too. */

	    ir1 = uub + llb + 2;
	    ir2 = uub + *m + 2;
	    i__1 = *n;
	    for (jc = 1; jc <= i__1; ++jc) {
		i__2 = uub + 1 - jc;
		for (jr = 1; jr <= i__2; ++jr) {
		    a[jr + jc * a_dim1] = 0.f;
/* L430: */
		}
/* Computing MAX   
   Computing MIN */
		i__3 = ir1, i__5 = ir2 - jc;
		i__2 = 1, i__4 = min(i__3,i__5);
		i__6 = *lda;
		for (jr = max(i__2,i__4); jr <= i__6; ++jr) {
		    a[jr + jc * a_dim1] = 0.f;
/* L440: */
		}
/* L450: */
	    }
	}
    }

    return 0;

/*     End of SLATMS */

} /* slatms_ */