zlatms.c

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

			}
		    }
		    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 = mnmin;
		for (j = 1; j <= i__1; ++j) {
		    i__4 = (1 - iskew) * j + ioffg + j * a_dim1;
		    i__2 = j;
		    z__1.r = d[i__2], z__1.i = 0.;
		    a[i__4].r = z__1.r, a[i__4].i = z__1.i;
/* L260: */
		}

		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.r = 0., extra.i = 0.;
			i__4 = (1 - iskew) * jc + 1 + ioffg + jc * a_dim1;
			ctemp.r = a[i__4].r, ctemp.i = a[i__4].i;
			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;
			if (zsym) {
			    ct.r = c.r, ct.i = c.i;
			    st.r = s.r, st.i = s.i;
			} else {
			    d_cnjg(&z__1, &ctemp);
			    ctemp.r = z__1.r, ctemp.i = z__1.i;
			    d_cnjg(&z__1, &c);
			    ct.r = z__1.r, ct.i = z__1.i;
			    d_cnjg(&z__1, &s);
			    st.r = z__1.r, st.i = z__1.i;
			}
			L__1 = *n - jc > k;
			zlarot_(&c_false, &c_true, &L__1, &il, &c, &s, &a[(1 
				- iskew) * jc + ioffg + jc * a_dim1], &ilda, &
				ctemp, &extra);
/* Computing MAX */
			i__4 = 1, i__2 = jc - k + 1;
			icol = max(i__4,i__2);
			i__4 = jc + 2 - icol;
			zlarot_(&c_true, &c_false, &c_true, &i__4, &ct, &st, &
				a[jc - iskew * icol + ioffg + icol * a_dim1], 
				&ilda, &dummy, &ctemp);

/*                    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) {
			    zlartg_(&a[jch - iskew * icol + ioffg + icol * 
				    a_dim1], &extra, &realc, &s, &dummy);
			    zlarnd_(&z__1, &c__5, &iseed[1]);
			    dummy.r = z__1.r, dummy.i = z__1.i;
			    z__1.r = realc * dummy.r, z__1.i = realc * 
				    dummy.i;
			    c.r = z__1.r, c.i = z__1.i;
			    z__1.r = s.r * dummy.r - s.i * dummy.i, z__1.i = 
				    s.r * dummy.i + s.i * dummy.r;
			    s.r = z__1.r, s.i = z__1.i;
			    i__3 = (1 - iskew) * jch + 1 + ioffg + jch * 
				    a_dim1;
			    ctemp.r = a[i__3].r, ctemp.i = a[i__3].i;
			    if (zsym) {
				ct.r = c.r, ct.i = c.i;
				st.r = s.r, st.i = s.i;
			    } else {
				d_cnjg(&z__1, &ctemp);
				ctemp.r = z__1.r, ctemp.i = z__1.i;
				d_cnjg(&z__1, &c);
				ct.r = z__1.r, ct.i = z__1.i;
				d_cnjg(&z__1, &s);
				st.r = z__1.r, st.i = z__1.i;
			    }
			    i__3 = k + 2;
			    zlarot_(&c_true, &c_true, &c_true, &i__3, &c, &s, 
				    &a[jch - iskew * icol + ioffg + icol * 
				    a_dim1], &ilda, &extra, &ctemp);
/* Computing MIN */
			    i__3 = *n + 1 - jch, i__5 = k + 2;
			    il = min(i__3,i__5);
			    extra.r = 0., extra.i = 0.;
			    L__1 = *n - jch > k;
			    zlarot_(&c_false, &c_true, &L__1, &il, &ct, &st, &
				    a[(1 - iskew) * jch + ioffg + jch * 
				    a_dim1], &ilda, &ctemp, &extra);
			    icol = jch;
/* L270: */
			}
/* L280: */
		    }
/* L290: */
		}

/*              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;
			if (zsym) {
/* Computing MAX */
			    i__2 = 1, i__4 = jc - uub;
			    i__1 = max(i__2,i__4);
			    for (jr = jc; jr >= i__1; --jr) {
				i__2 = jr + irow + jc * a_dim1;
				i__4 = jc - iskew * jr + ioffg + jr * a_dim1;
				a[i__2].r = a[i__4].r, a[i__2].i = a[i__4].i;
/* L300: */
			    }
			} else {
/* Computing MAX */
			    i__2 = 1, i__4 = jc - uub;
			    i__1 = max(i__2,i__4);
			    for (jr = jc; jr >= i__1; --jr) {
				i__2 = jr + irow + jc * a_dim1;
				d_cnjg(&z__1, &a[jc - iskew * jr + ioffg + jr 
					* a_dim1]);
				a[i__2].r = z__1.r, a[i__2].i = z__1.i;
/* L310: */
			    }
			}
/* L320: */
		    }
		    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) {
				i__4 = jr + jc * a_dim1;
				a[i__4].r = 0., a[i__4].i = 0.;
/* L330: */
			    }
/* L340: */
			}
		    }
		    if (ipackg == 5) {
			ipackg = ipack;
		    } else {
			ipackg = 0;
		    }
		}
	    }

/*           Ensure that the diagonal is real if Hermitian */

	    if (! zsym) {
		i__1 = *n;
		for (jc = 1; jc <= i__1; ++jc) {
		    irow = ioffst + (1 - iskew) * jc;
		    i__2 = irow + jc * a_dim1;
		    i__4 = irow + jc * a_dim1;
		    d__1 = a[i__4].r;
		    z__1.r = d__1, z__1.i = 0.;
		    a[i__2].r = z__1.r, a[i__2].i = z__1.i;
/* L350: */
		}
	    }

	}

    } 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 */

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

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

	    if (zsym) {
		zlagsy_(m, &llb, &d[1], &a[a_offset], lda, &iseed[1], &work[1]
			, &iinfo);
	    } else {
		zlaghe_(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) {
		    i__4 = i + j * a_dim1;
		    a[i__4].r = 0., a[i__4].i = 0.;
/* L360: */
		}
/* L370: */
	    }

	} 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) {
		    i__4 = i + j * a_dim1;
		    a[i__4].r = 0., a[i__4].i = 0.;
/* L380: */
		}
/* L390: */
	    }

	} 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;
		    }
		    i__4 = irow + icol * a_dim1;
		    i__3 = i + j * a_dim1;
		    a[i__4].r = a[i__3].r, a[i__4].i = a[i__3].i;
/* L400: */
		}
/* L410: */
	    }

	} 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;
		    }
		    i__4 = irow + icol * a_dim1;
		    i__3 = i + j * a_dim1;
		    a[i__4].r = a[i__3].r, a[i__4].i = a[i__3].i;
/* L420: */
		}
/* L430: */
	    }

	} 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) {
		    i__2 = i - j + uub + 1 + j * a_dim1;
		    i__4 = i + j * a_dim1;
		    a[i__2].r = a[i__4].r, a[i__2].i = a[i__4].i;
/* L440: */
		}
/* L450: */
	    }

	    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) {
		    i__4 = i - j + uub + 1 + j * a_dim1;
		    i__3 = i + j * a_dim1;
		    a[i__4].r = a[i__3].r, a[i__4].i = a[i__3].i;
/* L460: */
		}
/* L470: */
	    }
	}

/*        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) {
		    i__4 = jr + jc * a_dim1;
		    a[i__4].r = 0., a[i__4].i = 0.;
/* L480: */
		}
		irow = 0;
/* L490: */
	    }

	} 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) {
		    i__4 = jr + jc * a_dim1;
		    a[i__4].r = 0., a[i__4].i = 0.;
/* L500: */
		}
/* 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) {
		    i__2 = jr + jc * a_dim1;
		    a[i__2].r = 0., a[i__2].i = 0.;
/* L510: */
		}
/* L520: */
	    }
	}
    }

    return 0;

/*     End of ZLATMS */

} /* zlatms_ */