svd.c

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

/*     |________________________________________________________| */

static int
hsr5_(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 */
    if (*m >= *n) {
	goto L10;
    }
    fprintf(stderr,"ERROR: ARGUMENT M MUST BE .GE. N IN SUBROUTINE HSR5\n");
    abort();
L10:
    if (*m == *n) {
	goto L90;
    }
    k = *m;
L20:
    s = -(double)a[k + *n * a_dim1];
    i__1 = *m;
    for (i = *n; i <= i__1; ++i) {
/* L30: */
	a[i + k * a_dim1] = s * a[i + *n * a_dim1];
    }
    a[k + k * a_dim1] += (double)1.;
    l = *n;
L40:
    j = l;
    --l;
    if (l == 0) {
	goto L70;
    }
    s = (double)0.;
    i__1 = *m;
    for (i = j; i <= i__1; ++i) {
/* L50: */
	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) {
/* L60: */
	a[i + k * a_dim1] -= s * a[i + l * a_dim1];
    }
    goto L40;
L70:
    --k;
    if (k > *n) {
	goto L20;
    }
L90:
    hsr3_(&a[a_offset], la, m, n);
    return 0;
} /* hsr5_ */

/*      ________________________________________________________ */
/*     |                                                        | */
/*     |  COMPUTE SINGULAR VALUE DECOMPOSITION OF A BIDIAGONAL  | */
/*     | MATRIX:  A = Q TIMES DIAGONAL MATRIX TIMES P TRANSPOSE | */
/*     |                                                        | */
/*     |    INPUT:                                              | */
/*     |                                                        | */
/*     |         D     --DIAGONAL                               | */
/*     |                                                        | */
/*     |         N     --DIMENSION OF BIDIAGONAL MATRIX         | */
/*     |                                                        | */
/*     |         U     --OFF DIAGONAL                           | */
/*     |                                                        | */
/*     |         IU    --= 0 IF U IS SUPERDIAGONAL AND = 1 IF   | */
/*     |                 U IS SUBDIAGONAL                       | */
/*     |                                                        | */
/*     |         Q     --EITHER A SEGMENT OF AN IDENTITY MATRIX | */
/*     |                 OR A SEGMENT OF THE MATRIX USED TO     | */
/*     |                 BIDIAGONALIZE THE COEFFICIENT MATRIX   | */
/*     |                                                        | */
/*     |         LQ    --LEADING (ROW) DIMENSION OF ARRAY Q     | */
/*     |                                                        | */
/*     |         MQ    --NUMBER OF MATRIX ELEMENTS CONTAINED IN | */
/*     |                 EACH COLUMN OF Q                       | */
/*     |                                                        | */
/*     |         IQ    --AN INTEGER WHICH INDICATES WHETHER     | */
/*     |                 COLUMNS OF Q TO BE PROCESSED (= 0 MEANS| */
/*     |                 NO AND = 1 MEANS YES)                  | */
/*     |                                                        | */
/*     |         LP    --LEADING (ROW) DIMENSION OF ARRAY P     | */
/*     |                                                        | */
/*     |         MP    --NUMBER OF MATRIX ELEMENTS CONTAINED IN | */
/*     |                 EACH COLUMN OF P                       | */
/*     |                                                        | */
/*     |         IP    --AN INTEGER DEFINED LIKE IQ             | */
/*     |                                                        | */
/*     |         E     --WORK ARRAY WITH AT LEAST N ELEMENTS    | */
/*     |                                                        | */
/*     |         F     --WORK ARRAY WITH AT LEAST N ELEMENTS    | */
/*     |                                                        | */
/*     |    OUTPUT:                                             | */
/*     |                                                        | */
/*     |         Q     --Q FACTOR IN THE SINGULAR VALUE DECOMP. | */
/*     |                                                        | */
/*     |         D     --SINGULAR VALUES IN DECREASING ORDER    | */
/*     |                                                        | */
/*     |         P     --P FACTOR IN THE SINGULAR VALUE DECOMP. | */
/*     |                                                        | */
/*     |    BUILTIN FUNCTIONS: ABS,AMAX1,SIGN,SQRT              | */
/*     |    PACKAGE SUBROUTINES: EIG3,FGV,HSR1-HSR5,SCL,SFT,    | */
/*     |                         SNG0,SNG1,SORT2                | */
/*     |________________________________________________________| */

static int
singb_(double *d, int *n, double *u, int *iu,
       double *q, int *lq, int *mq, int *iq,
       double *p, int *lp, int *mp, int *ip, double *e, double *f)
{
    int one = 1;

    /* System generated locals */
    integer q_dim1, q_offset, p_dim1, p_offset, i__1, i__2;
    real r__1, r__2;

    /* Local variables */
    static real b, c;
    static integer g, h, i, j, k, l, m;
    static real r, s, t, v, w, x, y, z;
    static integer j0, k2, l1;
    static real t0, t1, t2, t3;
    static integer id, jl, ll, jr, ns;

    /* Parameter adjustments */
    --f;
    --e;
    p_dim1 = *lp;
    p_offset = p_dim1 + 1;
    p -= p_offset;
    q_dim1 = *lq;
    q_offset = q_dim1 + 1;
    q -= q_offset;
    --u;
    --d;

    /* Function Body */
    if (*n > 1) {
	goto L10;
    }
    if (*iq == 1) {
	q[q_dim1 + 1] = (double)1.;
    }
    if (*ip == 1) {
	p[p_dim1 + 1] = (double)1.;
    }
    if (d[1] >= (double)0.) {
	return 0;
    }
    d[1] = -(double)d[1];
    if (*iq == 1) {
	q[q_dim1 + 1] = (double)-1.;
    }
    return 0;
L10:
    jl = *iq;
    if (jl == 0) {
	jl = 3;
    }
    if (jl != 3) {
	jl = 1;
    }
    jr = *ip;
    if (jr == 0) {
	jr = 3;
    }
    if (jr != 3) {
	jr = 2;
    }
    if (*iu == 0) {
	goto L20;
    }
    i = jl;
    jl = jr;
    jr = i;
L20:
    j = 0;
    l = *n - 1;
    k2 = *n - 2;
    i__1 = l;
    for (i = 1; i <= i__1; ++i) {
	e[i] = (double)1.;
	f[i] = (double)1.;
	if (u[i] == (double)0.) {
	    j = i;
	}
/* L30: */
	if (d[i] == (double)0.) {
	    j = i;
	}
    }
    e[*n] = (double)1.;
    f[*n] = (double)1.;
    b = (double)3.5527136788005009e-15;
    t = (double)1.;
L40:
    t *= (double).5;
    s = t + (double)1.;
    if (s > (double)1.) {
	goto L40;
    }
    t0 = (double)1. / (t + t);
/* Computing 2nd power */
    r__1 = t + t;
    t2 = r__1 * r__1;
    ns = *n * 50;
    l1 = 0;
    k = *n;
    ll = 0;
    goto L70;
L50:
    j = 0;
    i__1 = l;
    for (i = 1; i <= i__1; ++i) {
	if (u[i] == (double)0.) {
	    j = i;
	}
/* L60: */
	if (d[i] == (double)0.) {
	    j = i;
	}
    }
L70:
    if (j == 0) {
	goto L140;
    }
    if (u[j] == (double)0.) {
	goto L140;
    }
/*     ------------------------------- */
/*     |*** ZERO DIAGONAL ELEMENT ***| */
/*     ------------------------------- */
    i = j;
    v = u[j];
    u[j] = (double)0.;
    s = -(double)v * (r__1 = e[j], fabs(r__1));
/*     --------------------------- */
/*     |*** PROCESS LEFT SIDE ***| */
/*     --------------------------- */
L80:
    h = i;
    ++i;
    r = (r__1 = e[i], fabs(r__1)) * d[i];
    fgv_(&x, &y, &t, &r, &s, &e[j], &e[i]);
    if (t == (double)1.) {
	goto L90;
    }
    d[i] -= y * v;
    sng0_(&q[q_offset], lq, mq, &p[p_offset], lp, mp, &jl, &j, &i, &x, &y);
    if (i == k) {
	goto L100;
    }
    v = x * u[i];
    s = -(double)v * (r__1 = e[j], fabs(r__1));
    goto L80;
L90:
    d[i] = d[i] * y - v;
    sng1_(&q[q_offset], lq, mq, &p[p_offset], lp, mp, &jl, &j, &i, &x, &y);
    if (i == k) {
	goto L100;
    }
    v = u[i];
    u[i] = v * y;
    s = -(double)v * (r__1 = e[j], fabs(r__1));
    goto L80;
L100:
    if (j == 1) {
	goto L130;
    }
    i = j - 1;
    s = u[i];
    u[i] = (double)0.;
/*     ---------------------------- */
/*     |*** PROCESS RIGHT SIDE ***| */
/*     ---------------------------- */
L110:
    h = i;
    --i;
    r = d[h];
    fgv_(&x, &y, &t, &r, &s, &f[h], &f[j]);
    if (t == (double)1.) {
	goto L120;
    }
    d[h] = r + x * s;
    sng0_(&q[q_offset], lq, mq, &p[p_offset], lp, mp, &jr, &h, &j, &x, &y);
    if (h == 1) {
	goto L130;
    }
    s = -(double)y * u[i];
    goto L110;
L120:
    d[h] = x * r + s;
    sng1_(&q[q_offset], lq, mq, &p[p_offset], lp, mp, &jr, &h, &j, &x, &y);
    if (h == 1) {
	goto L130;
    }
    s = -(double)u[i];
    u[i] = x * u[i];
    goto L110;
L130:
    scl_(&d[1], &u[1], n, &q[q_offset], lq, mq, &p[p_offset], lp, mp, &e[1], &
	    f[1], &b, &one, &k, &jl, &jr);
L140:
    ++j;
    if (j == k) {
	goto L320;
    }
    s = (double)0.;
    t = (double)0.;
/*     ----------------------------- */
/*     |*** SET ERROR TOLERANCE ***| */
/*     ----------------------------- */
    i__1 = l;
    for (i = j; i <= i__1; ++i) {
	x = (r__1 = d[i] * e[i] * (d[i] * f[i]), fabs(r__1));
	y = (r__1 = u[i] * e[i] * (u[i] * f[i + 1]), fabs(r__1));
/* Computing MAX */
	r__1 = max(s,x);
	s = dmax(r__1,y);
/* L150: */
	t = t + x + y;
    }
    x = (r__1 = e[k] * d[k] * (f[k] * d[k]), fabs(r__1));
    s = dmax(s,x);
    t3 = s * t2;
    t += x;
    if (t == (double)0.) {
	t = (double)1.;
    }
    t1 = t0 / t;
    goto L280;
L160:
    ++ll;
    if (ll > ns) {
	goto L530;
    }
    if (l > j) {
	goto L170;
    }
    s = (double)0.;
    t = (double)0.;
    goto L180;
L170:
    s = u[k2];
    t = e[k2];
/*     ----------------------- */
/*     |*** COMPUTE SHIFT ***| */
/*     ----------------------- */
L180:
    sft_(&c, &d[k], &d[l], &u[l], &s, &e[k], &e[l], &t, &f[k], &f[l]);
    v = (r__1 = e[j], fabs(r__1));
    w = (r__1 = f[j], fabs(r__1));
    t = d[j] * w;
    r = t * (d[j] * v) - c;
    s = t * (u[j] * v);
    id = 1;
    j0 = j;
    i = j;
    h = j - 1;
L190:
    g = h;
    h = i;
    ++i;
    z = (r__1 = f[i], fabs(r__1));
/*     ---------------------------- */
/*     |*** PROCESS RIGHT SIDE ***| */
/*     ---------------------------- */
    fgv_(&x, &y, &t, &r, &s, &f[h], &f[i]);
    v = d[h];
    if (t == (double)1.) {
	goto L210;
    }
    if (h == j) {
	goto L200;
    }
    t = u[g] + x * s;
    u[g] = t;
    if ((r__1 = t * e[g] * (t * f[h]), fabs(r__1)) > t3) {
	goto L200;
    }
    j = h;
    id = 1;
L200:
    r = v + x * u[h];
    s = x * d[i];
    u[h] -= y * v;
    sng0_(&q[q_offset], lq, mq, &p[p_offset], lp, mp, &jr, &h, &i, &x, &y);
    goto L230;
L210:
    if (h == j) {
	goto L220;
    }
    t = x * u[g] + s;
    u[g] = t;
    if ((r__1 = t * e[g] * (t * f[h]), fabs(r__1)) > t3) {
	goto L220;
    }
    j = h;
    id = 1;
L220:
    r = x * v + u[h];
    s = d[i];
    u[h] = y * u[h] - v;
    d[i] = y * s;
    sng1_(&q[q_offset], lq, mq, &p[p_offset], lp, mp, &jr, &h, &i, &x, &y);
/*     --------------------------- */
/*     |*** PROCESS LEFT SIDE ***| */
/*     --------------------------- */
L230:
    fgv_(&x, &y, &t, &r, &s, &e[h], &e[i]);
    if (t == (double)1.) {
	goto L250;
    }
    t = r + x * s;
    d[h] = t;
    if ((r__1 = t * e[h] * (t * f[h]), fabs(r__1)) > t3) {
	goto L240;
    }
    id = 0;
    j = h;
L240:
    r = u[h] + x * d[i];
    d[i] -= y * u[h];
    u[h] = r;
    sng0_(&q[q_offset], lq, mq, &p[p_offset], lp, mp, &jl, &h, &i, &x, &y);
    if (i == k) {
	goto L270;
    }
    s = x * u[i];
    goto L190;
L250:
    t = s + x * r;
    d[h] = t;
    if ((r__1 = t * e[h] * (t * f[h]), fabs(r__1)) > t3) {
	goto L260;
    }
    id = 0;
    j = h;
L260:
    r = d[i] + x * u[h];
    d[i] = y * d[i] - u[h];
    u[h] = r;
    sng1_(&q[q_offset], lq, mq, &p[p_offset], lp, mp, &jl, &h, &i, &x, &y);
    if (i == k) {
	goto L270;
    }
    s = u[i];
    u[i] = s * y;
    goto L190;
L270:
    scl_(&d[1], &u[1], n, &q[q_offset], lq, mq, &p[p_offset], lp, mp, &e[1], &
	    f[1], &b, &j0, &k, &jl, &jr);
    if (id == 0) {
	goto L70;
    }
L280:
    w = e[l];
    x = e[k];
    y = f[l];
    z = f[k];
    r = (r__1 = x * d[k] * (z * d[k]), fabs(r__1));
    s = (r__1 = w * d[l] * (y * d[l]), fabs(r__1));
    t = (r__1 = x * u[l] * (y * u[l]), fabs(r__1));
/*     ------------------------------ */
/*     |*** TEST FOR CONVERGENCE ***| */
/*     ------------------------------ */
    if (s * t1 * (t * t1) > (double)1.) {
	goto L160;
    }
    ++l1;
    if (l1 > 40) {
	goto L290;
    }
    if (t == (double)0.) {
	goto L290;
    }
    r += t;
    if (s / r * (t / r) > t2) {
	goto L160;
    }
L290:
    l1 = 0;
    if (s > r) {
	goto L310;
    }
    r = -(double)d[k] * fabs(x);
    s = u[l] * fabs(w);
    fgv_(&w, &y, &t, &r, &s, &e[l], &e[k]);
    x = e[k];
    if (t == (double)1.) {
	goto L300;
    }
    d[k] -= y * u[l];
    sng0_(&q[q_offset], lq, mq, &p[p_offset], lp, mp, &jl, &l, &k, &w, &y);
    goto L310;
L300:
    d[l] *= w;
    d[k] = y * d[k] - u[l];
    sng1_(&q[q_offset], lq, mq, &p[p_offset], lp, mp, &jl, &l, &k, &w, &y);
L310:
    r__1 = sqrt((fabs(x)));
    r__2 = sqrt((fabs(z)));
    t = r_sign(&r__1, &x) * d[k] * r_sign(&r__2, &z);
    if (t < (double)0.) {
	e[k] = -(double)e[k];
    }
    d[k] = fabs(t);
    k = l;
    l = k2;
    --k2;
    if (k > j) {
	goto L280;
    }
L320:
    x = e[k];
    z = f[k];
    r__1 = sqrt((fabs(x)));
    r__2 = sqrt((fabs(z)));
    t = r_sign(&r__1, &x) * d[k] * r_sign(&r__2, &z);
    if (t < (double)0.) {
	e[k] = -(double)e[k];
    }
    d[k] = fabs(t);
    k = l;
    l = k2;
    --k2;
    if (k > 1) {
	goto L50;
    }
    if (k == 0) {
	goto L330;
    }
    goto L320;
/*     ------------------------- */
/*     |*** RESCALE P AND Q ***| */
/*     ------------------------- */
L330:
    switch ((int)jl) {
	case 1:  goto L340;
	case 2:  goto L360;
	case 3:  goto L380;
    }
L340:
    i__1 = *n;
    for (j = 1; j <= i__1; ++j) {
	t = e[j];
	r__1 = sqrt((fabs(t)));
	t = r_sign(&r__1, &t);
	i__2 = *mq;
	for (i = 1; i <= i__2; ++i) {
/* L350: */
	    q[i + j * q_dim1] *= t;
	}
    }
    goto L380;
L360:
    i__2 = *n;
    for (j = 1; j <= i__2; ++j) {
	t = e[j];
	r__1 = sqrt((fabs(t)));
	t = r_sign(&r__1, &t);
	i__1 = *mp;
	for (i = 1; i <= i__1; ++i) {
/* L370: */
	    p[i + j * p_dim1] *= t;
	}
    }
L380:
    switch ((int)jr) {
	case 1:  goto L390;
	case 2:  goto L410;
	case 3:  goto L430;
    }
L390:
    i__1 = *n;
    for (j = 1; j <= i__1; ++j) {
	t = f[j];
	r__1 = sqrt((fabs(t)));
	t = r_sign(&r__1, &t);
	i__2 = *mq;
	for (i = 1; i <= i__2; ++i) {
/* L400: */
	    q[i + j * q_dim1] *= t;
	}
    }
    goto L430;
L410: