clatms.c

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

					doublereal)s.i;
				q__2.r = q__3.r * dummy.r - q__3.i * dummy.i, 
					q__2.i = q__3.r * dummy.i + q__3.i * 
					dummy.r;
				r_cnjg(&q__1, &q__2);
				s.r = q__1.r, s.i = q__1.i;
/* Computing MAX */
				i__4 = 1, i__5 = jch - jkl - jku;
				irow = max(i__4,i__5);
				il = ir + 2 - irow;
				extra.r = 0.f, extra.i = 0.f;
				L__1 = jch > jkl + jku;
				clarot_(&c_false, &L__1, &c_true, &il, &c, &s,
					 &a[irow - iskew * icol + ioffst + 
					icol * a_dim1], &ilda, &extra, &ctemp)
					;
				ic = icol;
				ir = irow;
			    }
/* L80: */
			}
/* L90: */
		    }
/* L100: */
		}

	    } 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.r = 0.f, extra.i = 0.f;
			angle = slarnd_(&c__1, &iseed[1]) * 
				6.2831853071795864769252867663f;
			d__1 = cos(angle);
			clarnd_(&q__2, &c__5, &iseed[1]);
			q__1.r = d__1 * q__2.r, q__1.i = d__1 * q__2.i;
			c.r = q__1.r, c.i = q__1.i;
			d__1 = sin(angle);
			clarnd_(&q__2, &c__5, &iseed[1]);
			q__1.r = d__1 * q__2.r, q__1.i = d__1 * q__2.i;
			s.r = q__1.r, s.i = q__1.i;
/* 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;
			    clarot_(&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) {
				clartg_(&a[jch - iskew * ic + ioffst + ic * 
					a_dim1], &extra, &realc, &s, &dummy);
				clarnd_(&q__1, &c__5, &iseed[1]);
				dummy.r = q__1.r, dummy.i = q__1.i;
				q__1.r = realc * dummy.r, q__1.i = realc * 
					dummy.i;
				c.r = q__1.r, c.i = q__1.i;
				q__1.r = s.r * dummy.r - s.i * dummy.i, 
					q__1.i = s.r * dummy.i + s.i * 
					dummy.r;
				s.r = q__1.r, s.i = q__1.i;
			    }
			    ic = max(1,ic);
/* Computing MIN */
			    i__5 = *n - 1, i__6 = jch + jku;
			    icol = min(i__5,i__6);
			    iltemp = jch + jku < *n;
			    ctemp.r = 0.f, ctemp.i = 0.f;
			    i__5 = icol + 2 - ic;
			    clarot_(&c_true, &ilextr, &iltemp, &i__5, &c, &s, 
				    &a[jch - iskew * ic + ioffst + ic * 
				    a_dim1], &ilda, &extra, &ctemp);
			    if (iltemp) {
				clartg_(&a[jch - iskew * icol + ioffst + icol 
					* a_dim1], &ctemp, &realc, &s, &dummy)
					;
				clarnd_(&q__1, &c__5, &iseed[1]);
				dummy.r = q__1.r, dummy.i = q__1.i;
				q__1.r = realc * dummy.r, q__1.i = realc * 
					dummy.i;
				c.r = q__1.r, c.i = q__1.i;
				q__1.r = s.r * dummy.r - s.i * dummy.i, 
					q__1.i = s.r * dummy.i + s.i * 
					dummy.r;
				s.r = q__1.r, s.i = q__1.i;
/* Computing MIN */
				i__5 = iendch, i__6 = jch + jkl + jku;
				il = min(i__5,i__6) + 2 - jch;
				extra.r = 0.f, extra.i = 0.f;
				L__1 = jch + jkl + jku <= iendch;
				clarot_(&c_false, &c_true, &L__1, &il, &c, &s,
					 &a[jch - iskew * icol + ioffst + 
					icol * a_dim1], &ilda, &ctemp, &extra)
					;
				ic = icol;
			    }
/* L110: */
			}
/* L120: */
		    }
/* L130: */
		}

		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.r = 0.f, extra.i = 0.f;
			angle = slarnd_(&c__1, &iseed[1]) * 
				6.2831853071795864769252867663f;
			d__1 = cos(angle);
			clarnd_(&q__2, &c__5, &iseed[1]);
			q__1.r = d__1 * q__2.r, q__1.i = d__1 * q__2.i;
			c.r = q__1.r, c.i = q__1.i;
			d__1 = sin(angle);
			clarnd_(&q__2, &c__5, &iseed[1]);
			q__1.r = d__1 * q__2.r, q__1.i = d__1 * q__2.i;
			s.r = q__1.r, s.i = q__1.i;
/* 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;
			    clarot_(&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) {
				clartg_(&a[ir - iskew * jch + ioffst + jch * 
					a_dim1], &extra, &realc, &s, &dummy);
				clarnd_(&q__1, &c__5, &iseed[1]);
				dummy.r = q__1.r, dummy.i = q__1.i;
				q__1.r = realc * dummy.r, q__1.i = realc * 
					dummy.i;
				c.r = q__1.r, c.i = q__1.i;
				q__1.r = s.r * dummy.r - s.i * dummy.i, 
					q__1.i = s.r * dummy.i + s.i * 
					dummy.r;
				s.r = q__1.r, s.i = q__1.i;
			    }
			    ir = max(1,ir);
/* Computing MIN */
			    i__5 = *m - 1, i__6 = jch + jkl;
			    irow = min(i__5,i__6);
			    iltemp = jch + jkl < *m;
			    ctemp.r = 0.f, ctemp.i = 0.f;
			    i__5 = irow + 2 - ir;
			    clarot_(&c_false, &ilextr, &iltemp, &i__5, &c, &s,
				     &a[ir - iskew * jch + ioffst + jch * 
				    a_dim1], &ilda, &extra, &ctemp);
			    if (iltemp) {
				clartg_(&a[irow - iskew * jch + ioffst + jch *
					 a_dim1], &ctemp, &realc, &s, &dummy);
				clarnd_(&q__1, &c__5, &iseed[1]);
				dummy.r = q__1.r, dummy.i = q__1.i;
				q__1.r = realc * dummy.r, q__1.i = realc * 
					dummy.i;
				c.r = q__1.r, c.i = q__1.i;
				q__1.r = s.r * dummy.r - s.i * dummy.i, 
					q__1.i = s.r * dummy.i + s.i * 
					dummy.r;
				s.r = q__1.r, s.i = q__1.i;
/* Computing MIN */
				i__5 = iendch, i__6 = jch + jkl + jku;
				il = min(i__5,i__6) + 2 - jch;
				extra.r = 0.f, extra.i = 0.f;
				L__1 = jch + jkl + jku <= iendch;
				clarot_(&c_true, &c_true, &L__1, &il, &c, &s, 
					&a[irow - iskew * jch + ioffst + jch *
					 a_dim1], &ilda, &ctemp, &extra);
				ir = irow;
			    }
/* L140: */
			}
/* L150: */
		    }
/* L160: */
		}

	    }

	} else {

/*           Symmetric -- A = U D U'   
             Hermitian -- 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 = mnmin;
		for (j = 1; j <= i__1; ++j) {
		    i__4 = (1 - iskew) * j + ioffg + j * a_dim1;
		    i__2 = j;
		    q__1.r = d[i__2], q__1.i = 0.f;
		    a[i__4].r = q__1.r, a[i__4].i = q__1.i;
/* L170: */
		}

		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.r = 0.f, extra.i = 0.f;
			i__2 = jc - iskew * (jc + 1) + ioffg + (jc + 1) * 
				a_dim1;
			ctemp.r = a[i__2].r, ctemp.i = a[i__2].i;
			angle = slarnd_(&c__1, &iseed[1]) * 
				6.2831853071795864769252867663f;
			d__1 = cos(angle);
			clarnd_(&q__2, &c__5, &iseed[1]);
			q__1.r = d__1 * q__2.r, q__1.i = d__1 * q__2.i;
			c.r = q__1.r, c.i = q__1.i;
			d__1 = sin(angle);
			clarnd_(&q__2, &c__5, &iseed[1]);
			q__1.r = d__1 * q__2.r, q__1.i = d__1 * q__2.i;
			s.r = q__1.r, s.i = q__1.i;
			if (csym) {
			    ct.r = c.r, ct.i = c.i;
			    st.r = s.r, st.i = s.i;
			} else {
			    r_cnjg(&q__1, &ctemp);
			    ctemp.r = q__1.r, ctemp.i = q__1.i;
			    r_cnjg(&q__1, &c);
			    ct.r = q__1.r, ct.i = q__1.i;
			    r_cnjg(&q__1, &s);
			    st.r = q__1.r, st.i = q__1.i;
			}
			L__1 = jc > k;
			clarot_(&c_false, &L__1, &c_true, &il, &c, &s, &a[
				irow - iskew * jc + ioffg + jc * a_dim1], &
				ilda, &extra, &ctemp);
/* Computing MIN */
			i__3 = k, i__5 = *n - jc;
			i__2 = min(i__3,i__5) + 1;
			clarot_(&c_true, &c_true, &c_false, &i__2, &ct, &st, &
				a[(1 - iskew) * jc + ioffg + jc * a_dim1], &
				ilda, &ctemp, &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) {
			    clartg_(&a[jch + 1 - iskew * (icol + 1) + ioffg + 
				    (icol + 1) * a_dim1], &extra, &realc, &s, 
				    &dummy);
			    clarnd_(&q__1, &c__5, &iseed[1]);
			    dummy.r = q__1.r, dummy.i = q__1.i;
			    q__2.r = realc * dummy.r, q__2.i = realc * 
				    dummy.i;
			    r_cnjg(&q__1, &q__2);
			    c.r = q__1.r, c.i = q__1.i;
			    q__3.r = -(doublereal)s.r, q__3.i = -(doublereal)
				    s.i;
			    q__2.r = q__3.r * dummy.r - q__3.i * dummy.i, 
				    q__2.i = q__3.r * dummy.i + q__3.i * 
				    dummy.r;
			    r_cnjg(&q__1, &q__2);
			    s.r = q__1.r, s.i = q__1.i;
			    i__3 = jch - iskew * (jch + 1) + ioffg + (jch + 1)
				     * a_dim1;
			    ctemp.r = a[i__3].r, ctemp.i = a[i__3].i;
			    if (csym) {
				ct.r = c.r, ct.i = c.i;
				st.r = s.r, st.i = s.i;
			    } else {
				r_cnjg(&q__1, &ctemp);
				ctemp.r = q__1.r, ctemp.i = q__1.i;
				r_cnjg(&q__1, &c);
				ct.r = q__1.r, ct.i = q__1.i;
				r_cnjg(&q__1, &s);
				st.r = q__1.r, st.i = q__1.i;
			    }
			    i__3 = k + 2;
			    clarot_(&c_true, &c_true, &c_true, &i__3, &c, &s, 
				    &a[(1 - iskew) * jch + ioffg + jch * 
				    a_dim1], &ilda, &ctemp, &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.r = 0.f, extra.i = 0.f;
			    L__1 = jch > k;
			    clarot_(&c_false, &L__1, &c_true, &il, &ct, &st, &
				    a[irow - iskew * jch + ioffg + jch * 
				    a_dim1], &ilda, &extra, &ctemp);
			    icol = jch;
/* L180: */
			}
/* L190: */
		    }
/* L200: */
		}

/*              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;
			if (csym) {
/* Computing MIN */
			    i__2 = *n, i__3 = jc + uub;
			    i__4 = min(i__2,i__3);
			    for (jr = jc; jr <= i__4; ++jr) {
				i__2 = jr + irow + jc * a_dim1;
				i__3 = jc - iskew * jr + ioffg + jr * a_dim1;
				a[i__2].r = a[i__3].r, a[i__2].i = a[i__3].i;
/* L210: */
			    }
			} else {
/* Computing MIN */
			    i__2 = *n, i__3 = jc + uub;
			    i__4 = min(i__2,i__3);
			    for (jr = jc; jr <= i__4; ++jr) {
				i__2 = jr + irow + jc * a_dim1;
				r_cnjg(&q__1, &a[jc - iskew * jr + ioffg + jr 
					* a_dim1]);
				a[i__2].r = q__1.r, a[i__2].i = q__1.i;
/* L220: */
			    }
			}
/* L230: */
		    }
		    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) {
				i__2 = jr + jc * a_dim1;
				a[i__2].r = 0.f, a[i__2].i = 0.f;
/* L240: */