zlatrs.c

来自「算断裂的」· C语言 代码 · 共 1,187 行 · 第 1/3 页

product is 1,   
                   call ZDOTU to perform the dot product. 
*/

		    if (upper) {
			i__3 = j - 1;
			zdotu_(&z__1, &i__3, &A(1,j), &c__1, &X(1),
				 &c__1);
			csumj.r = z__1.r, csumj.i = z__1.i;
		    } else if (j < *n) {
			i__3 = *n - j;
			zdotu_(&z__1, &i__3, &A(j+1,j), &c__1, &
				X(j + 1), &c__1);
			csumj.r = z__1.r, csumj.i = z__1.i;
		    }
		} else {

/*                 Otherwise, use in-line code for the dot
 product. */

		    if (upper) {
			i__3 = j - 1;
			for (i = 1; i <= j-1; ++i) {
			    i__4 = i + j * a_dim1;
			    z__3.r = A(i,j).r * uscal.r - A(i,j).i * 
				    uscal.i, z__3.i = A(i,j).r * uscal.i + A(i,j).i * uscal.r;
			    i__5 = i;
			    z__2.r = z__3.r * X(i).r - z__3.i * X(i).i, 
				    z__2.i = z__3.r * X(i).i + z__3.i * X(
				    i).r;
			    z__1.r = csumj.r + z__2.r, z__1.i = csumj.i + 
				    z__2.i;
			    csumj.r = z__1.r, csumj.i = z__1.i;
/* L130: */
			}
		    } else if (j < *n) {
			i__3 = *n;
			for (i = j + 1; i <= *n; ++i) {
			    i__4 = i + j * a_dim1;
			    z__3.r = A(i,j).r * uscal.r - A(i,j).i * 
				    uscal.i, z__3.i = A(i,j).r * uscal.i + A(i,j).i * uscal.r;
			    i__5 = i;
			    z__2.r = z__3.r * X(i).r - z__3.i * X(i).i, 
				    z__2.i = z__3.r * X(i).i + z__3.i * X(
				    i).r;
			    z__1.r = csumj.r + z__2.r, z__1.i = csumj.i + 
				    z__2.i;
			    csumj.r = z__1.r, csumj.i = z__1.i;
/* L140: */
			}
		    }
		}

		z__1.r = tscal, z__1.i = 0.;
		if (uscal.r == z__1.r && uscal.i == z__1.i) {

/*                 Compute x(j) := ( x(j) - CSUMJ ) / A(j,
j) if 1/A(j,j)   
                   was not used to scale the dotproduct. 
*/

		    i__3 = j;
		    i__4 = j;
		    z__1.r = X(j).r - csumj.r, z__1.i = X(j).i - 
			    csumj.i;
		    X(j).r = z__1.r, X(j).i = z__1.i;
		    i__3 = j;
		    xj = (d__1 = X(j).r, abs(d__1)) + (d__2 = d_imag(&X(j))
			    , abs(d__2));
		    if (nounit) {
			i__3 = j + j * a_dim1;
			z__1.r = tscal * A(j,j).r, z__1.i = tscal * A(j,j)
				.i;
			tjjs.r = z__1.r, tjjs.i = z__1.i;
		    } else {
			tjjs.r = tscal, tjjs.i = 0.;
			if (tscal == 1.) {
			    goto L160;
			}
		    }

/*                    Compute x(j) = x(j) / A(j,j), scalin
g if necessary. */

		    tjj = (d__1 = tjjs.r, abs(d__1)) + (d__2 = d_imag(&tjjs), 
			    abs(d__2));
		    if (tjj > smlnum) {

/*                       abs(A(j,j)) > SMLNUM: */

			if (tjj < 1.) {
			    if (xj > tjj * bignum) {

/*                             Scale X by 1/ab
s(x(j)). */

				rec = 1. / xj;
				zdscal_(n, &rec, &X(1), &c__1);
				*scale *= rec;
				xmax *= rec;
			    }
			}
			i__3 = j;
			zladiv_(&z__1, &X(j), &tjjs);
			X(j).r = z__1.r, X(j).i = z__1.i;
		    } else if (tjj > 0.) {

/*                       0 < abs(A(j,j)) <= SMLNUM: */

			if (xj > tjj * bignum) {

/*                          Scale x by (1/abs(x(j)
))*abs(A(j,j))*BIGNUM. */

			    rec = tjj * bignum / xj;
			    zdscal_(n, &rec, &X(1), &c__1);
			    *scale *= rec;
			    xmax *= rec;
			}
			i__3 = j;
			zladiv_(&z__1, &X(j), &tjjs);
			X(j).r = z__1.r, X(j).i = z__1.i;
		    } else {

/*                       A(j,j) = 0:  Set x(1:n) = 0, 
x(j) = 1, and   
                         scale = 0 and compute a solut
ion to A**T *x = 0. */

			i__3 = *n;
			for (i = 1; i <= *n; ++i) {
			    i__4 = i;
			    X(i).r = 0., X(i).i = 0.;
/* L150: */
			}
			i__3 = j;
			X(j).r = 1., X(j).i = 0.;
			*scale = 0.;
			xmax = 0.;
		    }
L160:
		    ;
		} else {

/*                 Compute x(j) := x(j) / A(j,j) - CSUMJ i
f the dot   
                   product has already been divided by 1/A
(j,j). */

		    i__3 = j;
		    zladiv_(&z__2, &X(j), &tjjs);
		    z__1.r = z__2.r - csumj.r, z__1.i = z__2.i - csumj.i;
		    X(j).r = z__1.r, X(j).i = z__1.i;
		}
/* Computing MAX */
		i__3 = j;
		d__3 = xmax, d__4 = (d__1 = X(j).r, abs(d__1)) + (d__2 = 
			d_imag(&X(j)), abs(d__2));
		xmax = max(d__3,d__4);
/* L170: */
	    }

	} else {

/*           Solve A**H * x = b */

	    i__1 = jlast;
	    i__2 = jinc;
	    for (j = jfirst; jinc < 0 ? j >= jlast : j <= jlast; j += jinc) {

/*              Compute x(j) = b(j) - sum A(k,j)*x(k).   
                                      k<>j */

		i__3 = j;
		xj = (d__1 = X(j).r, abs(d__1)) + (d__2 = d_imag(&X(j)), 
			abs(d__2));
		uscal.r = tscal, uscal.i = 0.;
		rec = 1. / max(xmax,1.);
		if (CNORM(j) > (bignum - xj) * rec) {

/*                 If x(j) could overflow, scale x by 1/(2
*XMAX). */

		    rec *= .5;
		    if (nounit) {
			d_cnjg(&z__2, &A(j,j));
			z__1.r = tscal * z__2.r, z__1.i = tscal * z__2.i;
			tjjs.r = z__1.r, tjjs.i = z__1.i;
		    } else {
			tjjs.r = tscal, tjjs.i = 0.;
		    }
		    tjj = (d__1 = tjjs.r, abs(d__1)) + (d__2 = d_imag(&tjjs), 
			    abs(d__2));
		    if (tjj > 1.) {

/*                       Divide by A(j,j) when scaling
 x if A(j,j) > 1.   

   Computing MIN */
			d__1 = 1., d__2 = rec * tjj;
			rec = min(d__1,d__2);
			zladiv_(&z__1, &uscal, &tjjs);
			uscal.r = z__1.r, uscal.i = z__1.i;
		    }
		    if (rec < 1.) {
			zdscal_(n, &rec, &X(1), &c__1);
			*scale *= rec;
			xmax *= rec;
		    }
		}

		csumj.r = 0., csumj.i = 0.;
		if (uscal.r == 1. && uscal.i == 0.) {

/*                 If the scaling needed for A in the dot 
product is 1,   
                   call ZDOTC to perform the dot product. 
*/

		    if (upper) {
			i__3 = j - 1;
			zdotc_(&z__1, &i__3, &A(1,j), &c__1, &X(1),
				 &c__1);
			csumj.r = z__1.r, csumj.i = z__1.i;
		    } else if (j < *n) {
			i__3 = *n - j;
			zdotc_(&z__1, &i__3, &A(j+1,j), &c__1, &
				X(j + 1), &c__1);
			csumj.r = z__1.r, csumj.i = z__1.i;
		    }
		} else {

/*                 Otherwise, use in-line code for the dot
 product. */

		    if (upper) {
			i__3 = j - 1;
			for (i = 1; i <= j-1; ++i) {
			    d_cnjg(&z__4, &A(i,j));
			    z__3.r = z__4.r * uscal.r - z__4.i * uscal.i, 
				    z__3.i = z__4.r * uscal.i + z__4.i * 
				    uscal.r;
			    i__4 = i;
			    z__2.r = z__3.r * X(i).r - z__3.i * X(i).i, 
				    z__2.i = z__3.r * X(i).i + z__3.i * X(
				    i).r;
			    z__1.r = csumj.r + z__2.r, z__1.i = csumj.i + 
				    z__2.i;
			    csumj.r = z__1.r, csumj.i = z__1.i;
/* L180: */
			}
		    } else if (j < *n) {
			i__3 = *n;
			for (i = j + 1; i <= *n; ++i) {
			    d_cnjg(&z__4, &A(i,j));
			    z__3.r = z__4.r * uscal.r - z__4.i * uscal.i, 
				    z__3.i = z__4.r * uscal.i + z__4.i * 
				    uscal.r;
			    i__4 = i;
			    z__2.r = z__3.r * X(i).r - z__3.i * X(i).i, 
				    z__2.i = z__3.r * X(i).i + z__3.i * X(
				    i).r;
			    z__1.r = csumj.r + z__2.r, z__1.i = csumj.i + 
				    z__2.i;
			    csumj.r = z__1.r, csumj.i = z__1.i;
/* L190: */
			}
		    }
		}

		z__1.r = tscal, z__1.i = 0.;
		if (uscal.r == z__1.r && uscal.i == z__1.i) {

/*                 Compute x(j) := ( x(j) - CSUMJ ) / A(j,
j) if 1/A(j,j)   
                   was not used to scale the dotproduct. 
*/

		    i__3 = j;
		    i__4 = j;
		    z__1.r = X(j).r - csumj.r, z__1.i = X(j).i - 
			    csumj.i;
		    X(j).r = z__1.r, X(j).i = z__1.i;
		    i__3 = j;
		    xj = (d__1 = X(j).r, abs(d__1)) + (d__2 = d_imag(&X(j))
			    , abs(d__2));
		    if (nounit) {
			d_cnjg(&z__2, &A(j,j));
			z__1.r = tscal * z__2.r, z__1.i = tscal * z__2.i;
			tjjs.r = z__1.r, tjjs.i = z__1.i;
		    } else {
			tjjs.r = tscal, tjjs.i = 0.;
			if (tscal == 1.) {
			    goto L210;
			}
		    }

/*                    Compute x(j) = x(j) / A(j,j), scalin
g if necessary. */

		    tjj = (d__1 = tjjs.r, abs(d__1)) + (d__2 = d_imag(&tjjs), 
			    abs(d__2));
		    if (tjj > smlnum) {

/*                       abs(A(j,j)) > SMLNUM: */

			if (tjj < 1.) {
			    if (xj > tjj * bignum) {

/*                             Scale X by 1/ab
s(x(j)). */

				rec = 1. / xj;
				zdscal_(n, &rec, &X(1), &c__1);
				*scale *= rec;
				xmax *= rec;
			    }
			}
			i__3 = j;
			zladiv_(&z__1, &X(j), &tjjs);
			X(j).r = z__1.r, X(j).i = z__1.i;
		    } else if (tjj > 0.) {

/*                       0 < abs(A(j,j)) <= SMLNUM: */

			if (xj > tjj * bignum) {

/*                          Scale x by (1/abs(x(j)
))*abs(A(j,j))*BIGNUM. */

			    rec = tjj * bignum / xj;
			    zdscal_(n, &rec, &X(1), &c__1);
			    *scale *= rec;
			    xmax *= rec;
			}
			i__3 = j;
			zladiv_(&z__1, &X(j), &tjjs);
			X(j).r = z__1.r, X(j).i = z__1.i;
		    } else {

/*                       A(j,j) = 0:  Set x(1:n) = 0, 
x(j) = 1, and   
                         scale = 0 and compute a solut
ion to A**H *x = 0. */

			i__3 = *n;
			for (i = 1; i <= *n; ++i) {
			    i__4 = i;
			    X(i).r = 0., X(i).i = 0.;
/* L200: */
			}
			i__3 = j;
			X(j).r = 1., X(j).i = 0.;
			*scale = 0.;
			xmax = 0.;
		    }
L210:
		    ;
		} else {

/*                 Compute x(j) := x(j) / A(j,j) - CSUMJ i
f the dot   
                   product has already been divided by 1/A
(j,j). */

		    i__3 = j;
		    zladiv_(&z__2, &X(j), &tjjs);
		    z__1.r = z__2.r - csumj.r, z__1.i = z__2.i - csumj.i;
		    X(j).r = z__1.r, X(j).i = z__1.i;
		}
/* Computing MAX */
		i__3 = j;
		d__3 = xmax, d__4 = (d__1 = X(j).r, abs(d__1)) + (d__2 = 
			d_imag(&X(j)), abs(d__2));
		xmax = max(d__3,d__4);
/* L220: */
	    }
	}
	*scale /= tscal;
    }

/*     Scale the column norms by 1/TSCAL for return. */

    if (tscal != 1.) {
	d__1 = 1. / tscal;
	dscal_(n, &d__1, &CNORM(1), &c__1);
    }

    return 0;

/*     End of ZLATRS */

} /* zlatrs_ */