svd.c

大型并行量子化学软件；支持密度泛函（DFT）。可以进行各种量子化学计算。支持CHARMM并行计算。非常具有应用价值。

	if (d[i] != (double)0.) {
	    goto L240;
	}
    }
    b[j] = d[k];
    a[j + k * a_dim1] = (double)0.;
    if (*ip == 0) {
	goto L70;
    }
    i__1 = *n;
    for (i = k; i <= i__1; ++i) {
/* L230: */
	p[i + j * p_dim1] = (double)0.;
    }
    goto L70;
L240:
    t = (r__1 = d[k], fabs(r__1));
    if (t != (double)0.) {
	u = (double)1. / t;
    }
    r = (double)1.;
    i__1 = *n;
    for (l = i; l <= i__1; ++l) {
	s = (r__1 = d[l], fabs(r__1));
	if (s <= t) {
	    goto L250;
	}
	u = (double)1. / s;
/* Computing 2nd power */
	r__1 = t * u;
	r = r * (r__1 * r__1) + (double)1.;
	t = s;
	goto L260;
L250:
/* Computing 2nd power */
	r__1 = s * u;
	r += r__1 * r__1;
L260:
	;
    }
    s = t * sqrt(r);
    r = d[k];
    u = (double)1. / sqrt(s * (s + fabs(r)));
    if (r < (double)0.) {
	s = -(double)s;
    }
    d[k] = u * (r + s);
    i__1 = *n;
    for (i = h; i <= i__1; ++i) {
/* L270: */
	d[i] *= u;
    }
    if (*ip == 0) {
	goto L290;
    }
    i__1 = *n;
    for (i = k; i <= i__1; ++i) {
/* L280: */
	p[i + j * p_dim1] = d[i];
    }
L290:
    b[j] = -(double)s;
    i__1 = *m;
    for (i = k; i <= i__1; ++i) {
/* L300: */
	b[i] = (double)0.;
    }
    i__1 = *n;
    for (l = k; l <= i__1; ++l) {
	t = d[l];
	a[j + l * a_dim1] = t;
	i__2 = *m;
	for (i = k; i <= i__2; ++i) {
/* L310: */
	    b[i] += t * a[i + l * a_dim1];
	}
    }
    i__2 = *n;
    for (l = k; l <= i__2; ++l) {
	t = d[l];
	i__1 = *m;
	for (i = k; i <= i__1; ++i) {
/* L320: */
	    a[i + l * a_dim1] -= t * b[i];
	}
    }
    goto L70;
L330:
    i__1 = *n;
    for (i = k; i <= i__1; ++i) {
/* L340: */
	d[i] = a[k + i * a_dim1];
    }
L350:
    j = k;
    k = h;
    i__1 = *n;
    for (i = k; i <= i__1; ++i) {
/* L360: */
	if (d[i] != (double)0.) {
	    goto L370;
	}
    }
    u = d[j];
    d[j] = (double)0.;
    goto L440;
L370:
    t = (r__1 = d[j], fabs(r__1));
    if (t != (double)0.) {
	u = (double)1. / t;
    }
    r = (double)1.;
    i__1 = *n;
    for (l = i; l <= i__1; ++l) {
	s = (r__1 = d[l], fabs(r__1));
	if (s <= t) {
	    goto L380;
	}
	u = (double)1. / s;
/* Computing 2nd power */
	r__1 = t * u;
	r = r * (r__1 * r__1) + (double)1.;
	t = s;
	goto L390;
L380:
/* Computing 2nd power */
	r__1 = s * u;
	r += r__1 * r__1;
L390:
	;
    }
    s = t * sqrt(r);
    r = d[j];
    u = (double)1. / sqrt(s * (s + fabs(r)));
    if (r < (double)0.) {
	s = -(double)s;
    }
    d[j] = u * (r + s);
    i__1 = *n;
    for (i = k; i <= i__1; ++i) {
/* L400: */
	d[i] *= u;
    }
    u = -(double)s;
    if (k > *m) {
	goto L470;
    }
    i__1 = *m;
    for (i = k; i <= i__1; ++i) {
/* L410: */
	b[i] = (double)0.;
    }
    i__1 = *n;
    for (l = j; l <= i__1; ++l) {
	t = d[l];
	a[j + l * a_dim1] = t;
	i__2 = *m;
	for (i = k; i <= i__2; ++i) {
/* L420: */
	    b[i] += t * a[i + l * a_dim1];
	}
    }
    i__2 = *n;
    for (l = j; l <= i__2; ++l) {
	t = d[l];
	i__1 = *m;
	for (i = k; i <= i__1; ++i) {
/* L430: */
	    a[i + l * a_dim1] -= t * b[i];
	}
    }
L440:
    h = k + 1;
    if (*ip == 0) {
	goto L460;
    }
    i__1 = *n;
    for (i = j; i <= i__1; ++i) {
/* L450: */
	p[i + j * p_dim1] = d[i];
    }
L460:
    d[j] = u;
    if (k < *m) {
	goto L490;
    }
    if (k > *m) {
	goto L620;
    }
    b[j] = a[*m + j * a_dim1];
    goto L330;
L470:
    i__1 = *n;
    for (i = j; i <= i__1; ++i) {
/* L480: */
	a[j + i * a_dim1] = d[i];
    }
    goto L440;
L490:
    i__1 = *m;
    for (i = h; i <= i__1; ++i) {
/* L500: */
	if (a[i + j * a_dim1] != (double)0.) {
	    goto L520;
	}
    }
    b[j] = a[k + j * a_dim1];
    a[k + j * a_dim1] = (double)0.;
    if (*iq == 0) {
	goto L330;
    }
    i__1 = *m;
    for (i = k; i <= i__1; ++i) {
/* L510: */
	q[i + j * q_dim1] = (double)0.;
    }
    goto L330;
L520:
    t = (r__1 = a[k + j * a_dim1], fabs(r__1));
    if (t != (double)0.) {
	u = (double)1. / t;
    }
    r = (double)1.;
    i__1 = *m;
    for (l = i; l <= i__1; ++l) {
	s = (r__1 = a[l + j * a_dim1], fabs(r__1));
	if (s <= t) {
	    goto L530;
	}
	u = (double)1. / s;
/* Computing 2nd power */
	r__1 = t * u;
	r = r * (r__1 * r__1) + (double)1.;
	t = s;
	goto L540;
L530:
/* Computing 2nd power */
	r__1 = s * u;
	r += r__1 * r__1;
L540:
	;
    }
    s = t * sqrt(r);
    r = a[k + j * a_dim1];
    u = (double)1. / sqrt(s * (s + fabs(r)));
    if (r < (double)0.) {
	s = -(double)s;
    }
    b[j] = -(double)s;
    a[k + j * a_dim1] = u * (r + s);
    i__1 = *m;
    for (i = h; i <= i__1; ++i) {
/* L550: */
	a[i + j * a_dim1] *= u;
    }
    if (*iq == 0) {
	goto L570;
    }
    i__1 = *m;
    for (i = k; i <= i__1; ++i) {
/* L560: */
	q[i + j * q_dim1] = a[i + j * a_dim1];
    }
L570:
    i__1 = *n;
    for (l = k; l <= i__1; ++l) {
	t = (double)0.;
	i__2 = *m;
	for (i = k; i <= i__2; ++i) {
/* L580: */
	    t += a[i + j * a_dim1] * a[i + l * a_dim1];
	}
	a[k + l * a_dim1] -= t * a[k + j * a_dim1];
	d[l] = a[k + l * a_dim1];
	i__2 = *m;
	for (i = h; i <= i__2; ++i) {
/* L590: */
	    a[i + l * a_dim1] -= t * a[i + j * a_dim1];
	}
/* L600: */
    }
    goto L350;
L610:
    d[*n] = a[*n + *n * a_dim1];
    b[*n - 1] = a[*n - 1 + *n * a_dim1];
L620:
    if (jq == 0) {
	goto L650;
    }
    if (*n > *m) {
	goto L640;
    }
    if (*n == *m) {
	goto L630;
    }
    if (jq == 1) {
	hsr3_(&q[q_offset], lq, m, n);
    }
    if (jq == 2) {
	hsr4_(&q[q_offset], lq, m, n);
    }
    if (jq == 3) {
	hsr5_(&q[q_offset], lq, m, n);
    }
    goto L650;
L630:
    hsr2_(&q[q_offset], lq, m);
    goto L650;
L640:
    hsr1_(&q[q_offset], lq, m);
L650:
    if (jp == 0) {
	return 0;
    }
    if (*n <= *m) {
	goto L660;
    }
    if (jp == 1) {
	hsr3_(&p[p_offset], lp, n, m);
    }
    if (jp == 2) {
	hsr4_(&p[p_offset], lp, n, m);
    }
    if (jp == 3) {
	hsr5_(&p[p_offset], lp, n, m);
    }
    return 0;
L660:
    hsr1_(&p[p_offset], lp, n);
    return 0;
} /* bidag2_ */


static int
hsr1_(double *a, int *la, int *n)
{
    /* System generated locals */
    integer a_dim1, a_offset, i__1;

    /* Local variables */
    static integer i, j, k, l, m;
    static real s;

    /* Parameter adjustments */
    a_dim1 = *la;
    a_offset = a_dim1 + 1;
    a -= a_offset;

    /* Function Body */
    a[a_dim1 + 1] = (double)1.;
    if (*n == 1) {
	return 0;
    }
    a[(a_dim1 << 1) + 1] = (double)0.;
    if (*n > 2) {
	goto L10;
    }
    a[a_dim1 + 2] = (double)0.;
    a[(a_dim1 << 1) + 2] = (double)1.;
    return 0;
L10:
    l = *n - 2;
    m = *n - 1;
    k = *n;
    s = a[k + l * a_dim1];
    a[*n + *n * a_dim1] = (double)1. - s * a[*n + l * a_dim1];
    a[m + *n * a_dim1] = -(double)s * a[m + l * a_dim1];
L20:
    j = m;
    m = l;
    --l;
    if (l == 0) {
	goto L50;
    }
    s = (double)0.;
    i__1 = *n;
    for (i = j; i <= i__1; ++i) {
/* L30: */
	s += a[i + l * a_dim1] * a[i + k * a_dim1];
    }
    a[m + k * a_dim1] = -(double)s * a[m + l * a_dim1];
    i__1 = *n;
    for (i = j; i <= i__1; ++i) {
/* L40: */
	a[i + k * a_dim1] -= s * a[i + l * a_dim1];
    }
    goto L20;
L50:
    a[k * a_dim1 + 1] = (double)0.;
    --k;
    m = k;
    l = k - 1;
    s = -(double)a[m + l * a_dim1];
    i__1 = *n;
    for (i = k; i <= i__1; ++i) {
/* L60: */
	a[i + k * a_dim1] = s * a[i + l * a_dim1];
    }
    a[k + k * a_dim1] += (double)1.;
    if (l > 1) {
	goto L20;
    }
    i__1 = *n;
    for (i = 2; i <= i__1; ++i) {
/* L70: */
	a[i + a_dim1] = (double)0.;
    }
    return 0;
} /* hsr1_ */

static int
hsr2_(double *a, int *la, int *n)
{
    /* System generated locals */
    integer a_dim1, a_offset, i__1;

    /* Local variables */
    static integer i, j, k, l, m;
    static real s;

    /* Parameter adjustments */
    a_dim1 = *la;
    a_offset = a_dim1 + 1;
    a -= a_offset;

    /* Function Body */
    if (*n > 1) {
	goto L10;
    }
    a[a_dim1 + 1] = (double)1.;
    return 0;
L10:
    l = *n - 2;
    m = *n - 1;
    k = *n;
    s = a[*n + m * a_dim1];
    a[*n + *n * a_dim1] = (double)1. - s * a[*n + m * a_dim1];
    a[m + *n * a_dim1] = -(double)s * a[m + m * a_dim1];
L20:
    j = m;
    --m;
    if (m == 0) {
	goto L50;
    }
    s = (double)0.;
    i__1 = *n;
    for (i = j; i <= i__1; ++i) {
/* L30: */
	s += a[i + m * a_dim1] * a[i + k * a_dim1];
    }
    a[m + k * a_dim1] = -(double)s * a[m + m * a_dim1];
    i__1 = *n;
    for (i = j; i <= i__1; ++i) {
/* L40: */
	a[i + k * a_dim1] -= s * a[i + m * a_dim1];
    }
    goto L20;
L50:
    --k;
    m = k;
    s = -(double)a[k + k * a_dim1];
    i__1 = *n;
    for (i = k; i <= i__1; ++i) {
/* L60: */
	a[i + k * a_dim1] = s * a[i + k * a_dim1];
    }
    a[k + k * a_dim1] += (double)1.;
    if (k > 1) {
	goto L20;
    }
    return 0;
} /* hsr2_ */

static int
hsr3_(double *a, int *la, int *m, int *n)
{
    /* System generated locals */
    integer a_dim1, a_offset, i__1;

    /* Local variables */
    static integer i, j, k, l;
    static real s;

    /* Parameter adjustments */
    a_dim1 = *la;
    a_offset = a_dim1 + 1;
    a -= a_offset;

    /* Function Body */
    k = *n;
    if (*m >= *n) {
	goto L10;
    }
    fprintf(stderr,"ERROR: ARGUMENT M MUST BE .GE. N IN SUBROUTINE HSR3\n");
    abort();
L10:
    s = -(double)a[k + k * a_dim1];
    i__1 = *m;
    for (i = k; i <= i__1; ++i) {
/* L20: */
	a[i + k * a_dim1] = s * a[i + k * a_dim1];
    }
    a[k + k * a_dim1] += (double)1.;
    if (k == 1) {
	return 0;
    }
    l = k;
L30:
    j = l;
    --l;
    if (l == 0) {
	goto L60;
    }
    s = (double)0.;
    i__1 = *m;
    for (i = j; i <= i__1; ++i) {
/* L40: */
	s += a[i + l * a_dim1] * a[i + k * a_dim1];
    }
    a[l + k * a_dim1] = -(double)s * a[l + l * a_dim1];
    i__1 = *m;
    for (i = j; i <= i__1; ++i) {
/* L50: */
	a[i + k * a_dim1] -= s * a[i + l * a_dim1];
    }
    goto L30;
L60:
    --k;
    goto L10;
} /* hsr3_ */

static int
hsr4_(double *a, int *la, int *m, int *n)
{
    /* System generated locals */
    integer a_dim1, a_offset, i__1, i__2;

    /* Local variables */
    static integer i, j, k, l;
    static real s;

    /* Parameter adjustments */
    a_dim1 = *la;
    a_offset = a_dim1 + 1;
    a -= a_offset;

    /* Function Body */
    k = *m;
    if (*m > *n) {
	goto L10;
    }
    fprintf(stderr,"ERROR: ARGUMENT M MUST BE .GE. N IN SUBROUTINE HSR4\n");
    abort();
L10:
    s = -(double)a[k + *n * a_dim1];
    i__1 = *m;
    for (i = *n; i <= i__1; ++i) {
/* L20: */
	a[i + k * a_dim1] = s * a[i + *n * a_dim1];
    }
    a[k + k * a_dim1] += (double)1.;
    l = *n;
L30:
    j = l;
    --l;
    if (l == 0) {
	goto L60;
    }
    s = (double)0.;
    i__1 = *m;
    for (i = j; i <= i__1; ++i) {
/* L40: */
	s += a[i + l * a_dim1] * a[i + k * a_dim1];
    }
    a[l + k * a_dim1] = -(double)s * a[l + l * a_dim1];
    i__1 = *m;
    for (i = j; i <= i__1; ++i) {
/* L50: */
	a[i + k * a_dim1] -= s * a[i + l * a_dim1];
    }
    goto L30;
L60:
    --k;
    if (k > *n) {
	goto L10;
    }
    k = *m - *n;
    i__1 = k;
    for (j = 1; j <= i__1; ++j) {
	i__2 = *m;
	for (i = 1; i <= i__2; ++i) {
/* L70: */
	    a[i + j * a_dim1] = a[i + (j + *n) * a_dim1];
	}
    }
    return 0;
} /* hsr4_ */

/*      ________________________________________________________ */
/*     |                                                        | */
/*     |    PACKAGE SUBROUTINES: HSR3                           | */