svd.c

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

    i__2 = *n;
    for (j = 1; j <= i__2; ++j) {
	t = f[j];
	r__1 = sqrt((fabs(t)));
	t = r_sign(&r__1, &t);
	i__1 = *mp;
	for (i = 1; i <= i__1; ++i) {
/* L420: */
	    p[i + j * p_dim1] *= t;
	}
    }
/*     -------------------------------- */
/*     |*** REORDER THE EIGENPAIRS ***| */
/*     -------------------------------- */
L430:
    sort2_(&d[1], &e[1], &f[1], n);
    i__1 = *n;
    for (i = 1; i <= i__1; ++i) {
	j = e[i];
/* L440: */
	f[j] = (real) i;
    }
    m = *n + 1;
    i__1 = *n;
    for (j = 1; j <= i__1; ++j) {
	l = m - j;
	k = e[l];
	i = f[j];
	e[i] = (real) k;
	f[k] = (real) i;
	t = d[j];
	d[j] = d[k];
	d[k] = t;
	if (jl == 1) {
	    goto L450;
	}
	if (jr != 1) {
	    goto L480;
	}
L450:
	s = (double)0.;
	i__2 = *mq;
	for (i = 1; i <= i__2; ++i) {
	    t = q[i + k * q_dim1];
	    q[i + k * q_dim1] = q[i + j * q_dim1];
	    q[i + j * q_dim1] = t;
/* L460: */
	    s += t * t;
	}
	s = (double)1. / sqrt(s);
	i__2 = *mq;
	for (i = 1; i <= i__2; ++i) {
/* L470: */
	    q[i + j * q_dim1] = s * q[i + j * q_dim1];
	}
L480:
	if (jr == 2) {
	    goto L490;
	}
	if (jl != 2) {
	    goto L520;
	}
L490:
	s = (double)0.;
	i__2 = *mp;
	for (i = 1; i <= i__2; ++i) {
	    t = p[i + k * p_dim1];
	    p[i + k * p_dim1] = p[i + j * p_dim1];
	    p[i + j * p_dim1] = t;
/* L500: */
	    s += t * t;
	}
	s = (double)1. / sqrt(s);
	i__2 = *mp;
	for (i = 1; i <= i__2; ++i) {
/* L510: */
	    p[i + j * p_dim1] = s * p[i + j * p_dim1];
	}
L520:
	;
    }
    e[1] = (real) (*n);
    return 0;
L530:
    k = *n - k + 1;
    fprintf(stderr,"SINCE THE STOPPING CRITERION NOT SATISFIED\n");
    fprintf(stderr,"AFTER %d ITERATIONS, WE STOP WHILE COMPUTING\n", ns);
    fprintf(stderr,"EIGENVALUE NUMBER %d", k);
    e[1] = (real) k;
    return 0;
} /* singb_ */

/* % */
static int
fgv_(double *x, double *y, double *s, double *p, double *q,
     double *a, double *b)
{
    /* System generated locals */
    real r__1;

    /* Local variables */
    static real c, r, t;

    if (fabs(*p) > fabs(*q)) {
	goto L10;
    }
    if (*q == (double)0.) {
	goto L110;
    }
    r = *a / *b;
    *s = *p / *q;
    t = fabs(r) * *s * *s;
    if (t < (double)1.) {
	goto L70;
    }
    t /= t + (double)1.;
    r = *b / *a;
    *s = *q / *p;
    goto L20;
L10:
    r = *b / *a;
    *s = *q / *p;
    t = fabs(r) * *s * *s;
    if (t > (double)1.) {
	goto L60;
    }
    t = (double)1. / (t + (double)1.);
L20:
    r__1 = *a * t;
    *a = r_sign(&r__1, p);
    if (r_sign(&r, p) == r) {
	goto L40;
    }
    *b = -(double)(r__1 = *b * t, fabs(r__1));
    goto L50;
L40:
    *b = (r__1 = *b * t, fabs(r__1));
L50:
    *y = *s;
    *x = fabs(r) * *s;
    *s = (double)0.;
    return 0;
L60:
    t /= t + (double)1.;
    r = *a / *b;
    *s = *p / *q;
    goto L80;
L70:
    t = (double)1. / (t + (double)1.);
L80:
    c = *a;
    r__1 = *b * t;
    *a = r_sign(&r__1, q);
    if (r_sign(&r, q) == r) {
	goto L90;
    }
    *b = -(double)(r__1 = c * t, fabs(r__1));
    goto L100;
L90:
    *b = (r__1 = c * t, fabs(r__1));
L100:
    *y = *s;
    *x = fabs(r) * *s;
    *s = (double)1.;
    return 0;
L110:
    *x = (double)0.;
    *y = (double)0.;
    *s = (double)0.;
    return 0;
} /* fgv_ */

/* % */
static int
sng0_(double *q, int *lq, int *m, double *p, int *lp, int *n,
      int *l, int *j, int *k, double *x, double *y)
{
    /* System generated locals */
    integer q_dim1, q_offset, p_dim1, p_offset, i__1;

    /* Local variables */
    static integer i;
    static real s, t;

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

    /* Function Body */
    switch ((int)*l) {
	case 1:  goto L10;
	case 2:  goto L30;
	case 3:  goto L50;
    }
L10:
    i__1 = *m;
    for (i = 1; i <= i__1; ++i) {
	t = q[i + *j * q_dim1];
	s = q[i + *k * q_dim1];
	q[i + *j * q_dim1] = t + *x * s;
/* L20: */
	q[i + *k * q_dim1] = s - *y * t;
    }
    return 0;
L30:
    i__1 = *n;
    for (i = 1; i <= i__1; ++i) {
	t = p[i + *j * p_dim1];
	s = p[i + *k * p_dim1];
	p[i + *j * p_dim1] = t + *x * s;
/* L40: */
	p[i + *k * p_dim1] = s - *y * t;
    }
L50:
    return 0;
} /* sng0_ */

/* % */
static int
sng1_(double *q, int *lq, int *m, double *p,
      int *lp, int *n, int *l, int *j, int *k, double *x, double *y)
{
    /* System generated locals */
    integer q_dim1, q_offset, p_dim1, p_offset, i__1;

    /* Local variables */
    static integer i;
    static real s, t;

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

    /* Function Body */
    switch ((int)*l) {
	case 1:  goto L10;
	case 2:  goto L30;
	case 3:  goto L50;
    }
L10:
    i__1 = *m;
    for (i = 1; i <= i__1; ++i) {
	t = q[i + *j * q_dim1];
	s = q[i + *k * q_dim1];
	q[i + *j * q_dim1] = *x * t + s;
/* L20: */
	q[i + *k * q_dim1] = *y * s - t;
    }
    return 0;
L30:
    i__1 = *n;
    for (i = 1; i <= i__1; ++i) {
	t = p[i + *j * p_dim1];
	s = p[i + *k * p_dim1];
	p[i + *j * p_dim1] = *x * t + s;
/* L40: */
	p[i + *k * p_dim1] = *y * s - t;
    }
L50:
    return 0;
} /* sng1_ */

/* % */
static int
sft_(double *s, double *a, double *b, double *c,
     double *d, double *e2, double *e1, double *e0, double *f2, double *f1)
{
    static real w, x, y, z, g0, g1, g2, h1, h2;

    g0 = fabs(*e0);
    g1 = fabs(*e1);
    g2 = fabs(*e2);
    h1 = fabs(*f1);
    h2 = fabs(*f2);
    w = *a * g2 * (*a * h2) + *c * g1 * (*c * h2);
    x = *b * g1 * (*b * h1) + *d * g0 * (*d * h1);
    y = *b * g1 * (*c * h1);
    z = *b * g1 * (*c * h2);
    eig3_(s, s, &x, &w, &y, &z);
    return 0;
} /* sft_ */

/* % */
static int
scl_(double *d, double *u, int *n, double *q,
     int *lq, int *mq, double *p, int *lp, int *mp,
     double *e, double *f, double *b,
     int *j, int *k, int *jl, 
     int *jr)
{
    /* System generated locals */
    integer p_dim1, p_offset, q_dim1, q_offset, i__1, i__2;
    real r__1, r__2;

    /* Local variables */
    static integer h, i, l, m;
    static real r, s, t, v;

    /* 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 */
    t = (double)1.;
    i__1 = *k;
    for (i = *j; i <= i__1; ++i) {
	if ((r__1 = e[i], fabs(r__1)) < t) {
	    t = (r__2 = e[i], fabs(r__2));
	}
/* L10: */
	if ((r__1 = f[i], fabs(r__1)) < t) {
	    t = (r__2 = f[i], fabs(r__2));
	}
    }
    if (t > *b) {
	return 0;
    }
/*     ---------------------------- */
/*     |*** RESCALE THE MATRIX ***| */
/*     ---------------------------- */
    r = e[*j];
    r__1 = sqrt((fabs(r)));
    r = r_sign(&r__1, &r);
    e[*j] = (double)1.;
    s = f[*j];
    r__1 = sqrt((fabs(s)));
    s = r_sign(&r__1, &s);
    f[*j] = (double)1.;
    d[*j] = r * d[*j] * s;
    t = r;
    if (*jl == 1) {
	goto L20;
    }
    if (*jr != 1) {
	goto L40;
    }
    t = s;
L20:
    i__1 = *mq;
    for (i = 1; i <= i__1; ++i) {
/* L30: */
	q[i + *j * q_dim1] *= t;
    }
L40:
    t = s;
    if (*jr == 2) {
	goto L50;
    }
    if (*jl != 2) {
	goto L70;
    }
    t = r;
L50:
    i__1 = *mp;
    for (i = 1; i <= i__1; ++i) {
/* L60: */
	p[i + *j * p_dim1] *= t;
    }
L70:
    l = *j + 1;
    h = *j;
    i__1 = *k;
    for (m = l; m <= i__1; ++m) {
	t = e[m];
	r__1 = sqrt((fabs(t)));
	t = r_sign(&r__1, &t);
	e[m] = (double)1.;
	s = f[m];
	r__1 = sqrt((fabs(s)));
	s = r_sign(&r__1, &s);
	f[m] = (double)1.;
	d[m] = s * d[m] * t;
	u[h] = r * u[h] * s;
	h = m;
	r = t;
	v = t;
	if (*jl == 1) {
	    goto L80;
	}
	if (*jr != 1) {
	    goto L100;
	}
	v = s;
L80:
	i__2 = *mq;
	for (i = 1; i <= i__2; ++i) {
/* L90: */
	    q[i + m * q_dim1] *= v;
	}
L100:
	v = s;
	if (*jr == 2) {
	    goto L110;
	}
	if (*jl != 2) {
	    goto L130;
	}
	v = t;
L110:
	i__2 = *mp;
	for (i = 1; i <= i__2; ++i) {
/* L120: */
	    p[i + m * p_dim1] *= v;
	}
L130:
	;
    }
    return 0;
} /* scl_ */

static int
eig3_(double *ea, double *eb, double *a, double *b, double *y, double *z)
{
    /* Local variables */
    static real c, s, t;

    t = (*b - *a) * (double).5;
    c = sqrt((fabs(*y))) * sqrt((fabs(*z)));
    if (fabs(t) > fabs(c)) {
	goto L30;
    }
    if (c != (double)0.) {
	goto L10;
    }
    *ea = *a;
    *eb = *b;
    return 0;
L10:
    t /= fabs(c);
    s = fabs(c) / (fabs(t) + sqrt(t * t + (double)1.));
    if (t < (double)0.) {
	goto L20;
    }
    *ea = *a - s;
    *eb = *b + s;
    return 0;
L20:
    *ea = *a + s;
    *eb = *b - s;
    return 0;
L30:
    t = fabs(c) / t;
    s = t * fabs(c) / (sqrt(t * t + (double)1.) + (double)1.);
    *ea = *a - s;
    *eb = *b + s;
    return 0;
} /* eig3_ */

/*      ________________________________________________________ */
/*     |                                                        | */
/*     |            SORT AN ARRAY IN INCREASING ORDER           | */
/*     |                                                        | */
/*     |    INPUT:                                              | */
/*     |                                                        | */
/*     |         X     --ARRAY OF NUMBERS                       | */
/*     |                                                        | */
/*     |         W     --WORKING ARRAY (LENGTH  AT LEAST N)     | */
/*     |                                                        | */
/*     |         N     --NUMBER OF ARRAY ELEMENTS TO SORT       | */
/*     |                                                        | */
/*     |    OUTPUT:                                             | */
/*     |                                                        | */
/*     |         X     --ORIGINAL ARRAY                         | */
/*     |                                                        | */
/*     |         Y     --INDICES OF X GIVING INCREASING ORDER   | */
/*     |________________________________________________________| */

static int
sort2_(double *x, double *y, double *w, int *n)
{
    static integer i, j, k, l, m, p, q;
    static real s, t;

    /* Parameter adjustments */
    --w;
    --y;
    --x;

    /* Function Body */
    i = 1;
L10:
    k = i;
L20:
    j = i;
    y[i] = (real) i;
    ++i;
    if (j == *n) {
	goto L30;
    }
    if (x[i] >= x[j]) {
	goto L20;
    }
    w[k] = (real) i;
    goto L10;
L30:
    if (k == 1) {
	return 0;
    }
    w[k] = (real) (*n + 1);
L40:
    m = 1;
    l = 1;
L50:
    i = l;
    if (i > *n) {
	goto L120;
    }
    p = y[i];
    s = x[p];
    j = w[i];
    k = j;
    if (j > *n) {
	goto L100;
    }
    q = y[j];
    t = x[q];
    l = w[j];
    y[i] = (real) l;
L60:
    if (s > t) {
	goto L70;
    }
    w[m] = (real) p;
    ++m;
    ++i;
    if (i == k) {
	goto L80;
    }
    p = y[i];
    s = x[p];
    goto L60;
L70:
    w[m] = (real) q;
    ++m;
    ++j;
    if (j == l) {
	goto L110;
    }
    q = y[j];
    t = x[q];
    goto L60;
L80:
    w[m] = (real) q;
    k = m + l - j;
    i = j - m;
L90:
    ++m;
    if (m == k) {
	goto L50;
    }
    w[m] = y[m + i];
    goto L90;
L100:
    y[i] = (real) j;
    l = j;
L110:
    w[m] = (real) p;
    k = m + k - i;
    i -= m;
    goto L90;
L120:
    i = 1;
L130:
    k = i;
    j = y[i];
L140:
    y[i] = w[i];
    ++i;
    if (i < j) {
	goto L140;
    }
    w[k] = (real) i;
    if (i <= *n) {
	goto L130;
    }
    if (k > 1) {
	goto L40;
    }
    return 0;
} /* sort2_ */