linpackd.c

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

/*     SCALES A VECTOR BY A CONSTANT. */
/*     USES UNROLLED LOOPS FOR INCREMENT EQUAL TO ONE. */
/*     JACK DONGARRA, LINPACK, 3/11/78. */


    /* Parameter adjustments */
    --dx;

    /* Function Body */
    if (*n <= 0) {
	return 0;
    }
    if (*incx == 1) {
	goto L20;
    }

/*        CODE FOR INCREMENT NOT EQUAL TO 1 */

    nincx = *n * *incx;
    i__1 = nincx;
    i__2 = *incx;
    for (i = 1; i__2 < 0 ? i >= i__1 : i <= i__1; i += i__2) {
	dx[i] = *da * dx[i];
/* L10: */
    }
    return 0;

/*        CODE FOR INCREMENT EQUAL TO 1 */


/*        CLEAN-UP LOOP */

L20:
    m = *n % 5;
    if (m == 0) {
	goto L40;
    }
    i__2 = m;
    for (i = 1; i <= i__2; ++i) {
	dx[i] = *da * dx[i];
/* L30: */
    }
    if (*n < 5) {
	return 0;
    }
L40:
    mp1 = m + 1;
    i__2 = *n;
    for (i = mp1; i <= i__2; i += 5) {
	dx[i] = *da * dx[i];
	dx[i + 1] = *da * dx[i + 1];
	dx[i + 2] = *da * dx[i + 2];
	dx[i + 3] = *da * dx[i + 3];
	dx[i + 4] = *da * dx[i + 4];
/* L50: */
    }
    return 0;
} /* dscal_ */

/* Subroutine */ int dswap_(n, dx, incx, dy, incy)
integer *n;
doublereal *dx;
integer *incx;
doublereal *dy;
integer *incy;
{
    /* System generated locals */
    integer i__1;

    /* Local variables */
    static integer i, m;
    static doublereal dtemp;
    static integer ix, iy, mp1;


/*     INTERCHANGES TWO VECTORS. */
/*     USES UNROLLED LOOPS FOR INCREMENTS EQUAL ONE. */
/*     JACK DONGARRA, LINPACK, 3/11/78. */


    /* Parameter adjustments */
    --dy;
    --dx;

    /* Function Body */
    if (*n <= 0) {
	return 0;
    }
    if (*incx == 1 && *incy == 1) {
	goto L20;
    }

/*       CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS NOT EQUAL */
/*         TO 1 */

    ix = 1;
    iy = 1;
    if (*incx < 0) {
	ix = (-(*n) + 1) * *incx + 1;
    }
    if (*incy < 0) {
	iy = (-(*n) + 1) * *incy + 1;
    }
    i__1 = *n;
    for (i = 1; i <= i__1; ++i) {
	dtemp = dx[ix];
	dx[ix] = dy[iy];
	dy[iy] = dtemp;
	ix += *incx;
	iy += *incy;
/* L10: */
    }
    return 0;

/*       CODE FOR BOTH INCREMENTS EQUAL TO 1 */


/*       CLEAN-UP LOOP */

L20:
    m = *n % 3;
    if (m == 0) {
	goto L40;
    }
    i__1 = m;
    for (i = 1; i <= i__1; ++i) {
	dtemp = dx[i];
	dx[i] = dy[i];
	dy[i] = dtemp;
/* L30: */
    }
    if (*n < 3) {
	return 0;
    }
L40:
    mp1 = m + 1;
    i__1 = *n;
    for (i = mp1; i <= i__1; i += 3) {
	dtemp = dx[i];
	dx[i] = dy[i];
	dy[i] = dtemp;
	dtemp = dx[i + 1];
	dx[i + 1] = dy[i + 1];
	dy[i + 1] = dtemp;
	dtemp = dx[i + 2];
	dx[i + 2] = dy[i + 2];
	dy[i + 2] = dtemp;
/* L50: */
    }
    return 0;
} /* dswap_ */

/* Subroutine */ int dsvdc_(x, ldx, n, p, s, e, u, ldu, v, ldv, work, job, 
	info)
doublereal *x;
integer *ldx, *n, *p;
doublereal *s, *e, *u;
integer *ldu;
doublereal *v;
integer *ldv;
doublereal *work;
integer *job, *info;
{
    /* System generated locals */
    integer x_dim1, x_offset, u_dim1, u_offset, v_dim1, v_offset, i__1, i__2, 
	    i__3;
    doublereal d__1, d__2, d__3, d__4, d__5, d__6, d__7;

    /* Local variables */
    static integer kase;
    extern doublereal ddot_();
    static integer jobu, iter;
    extern /* Subroutine */ int drot_();
    static doublereal test;
    extern doublereal dnrm2_();
    static integer nctp1;
    static doublereal b, c;
    static integer nrtp1;
    static doublereal f, g;
    static integer i, j, k, l, m;
    static doublereal t, scale;
    extern /* Subroutine */ int dscal_();
    static doublereal shift;
    extern /* Subroutine */ int dswap_(), drotg_();
    static integer maxit;
    extern /* Subroutine */ int daxpy_();
    static logical wantu, wantv;
    static doublereal t1, ztest, el;
    static integer kk;
    static doublereal cs;
    static integer ll, mm, ls;
    static doublereal sl;
    static integer lu;
    static doublereal sm, sn;
    static integer lm1, mm1, lp1, mp1, nct, ncu, lls, nrt;
    static doublereal emm1, smm1;



/*     DSVDC IS A SUBROUTINE TO REDUCE A DOUBLE PRECISION NXP MATRIX X */
/*     BY ORTHOGONAL TRANSFORMATIONS U AND V TO DIAGONAL FORM.  THE */
/*     DIAGONAL ELEMENTS S(I) ARE THE SINGULAR VALUES OF X.  THE */
/*     COLUMNS OF U ARE THE CORRESPONDING LEFT SINGULAR VECTORS, */
/*     AND THE COLUMNS OF V THE RIGHT SINGULAR VECTORS. */

/*     ON ENTRY */

/*         X         DOUBLE PRECISION(LDX,P), WHERE LDX.GE.N. */
/*                   X CONTAINS THE MATRIX WHOSE SINGULAR VALUE */
/*                   DECOMPOSITION IS TO BE COMPUTED.  X IS */
/*                   DESTROYED BY DSVDC. */

/*         LDX       INTEGER. */
/*                   LDX IS THE LEADING DIMENSION OF THE ARRAY X. */

/*         N         INTEGER. */
/*                   N IS THE NUMBER OF COLUMNS OF THE MATRIX X. */

/*         P         INTEGER. */
/*                   P IS THE NUMBER OF ROWS OF THE MATRIX X. */

/*         LDU       INTEGER. */
/*                   LDU IS THE LEADING DIMENSION OF THE ARRAY U. */
/*                   (SEE BELOW). */

/*         LDV       INTEGER. */
/*                   LDV IS THE LEADING DIMENSION OF THE ARRAY V. */
/*                   (SEE BELOW). */

/*         WORK      DOUBLE PRECISION(N). */
/*                   WORK IS A SCRATCH ARRAY. */

/*         JOB       INTEGER. */
/*                   JOB CONTROLS THE COMPUTATION OF THE SINGULAR */
/*                   VECTORS.  IT HAS THE DECIMAL EXPANSION AB */
/*                   WITH THE FOLLOWING MEANING */

/*                        A.EQ.0    DO NOT COMPUTE THE LEFT SINGULAR */
/*                                  VECTORS. */
/*                        A.EQ.1    RETURN THE N LEFT SINGULAR VECTORS */
/*                                  IN U. */
/*                        A.GE.2    RETURN THE FIRST MIN(N,P) SINGULAR */
/*                                  VECTORS IN U. */
/*                        B.EQ.0    DO NOT COMPUTE THE RIGHT SINGULAR */
/*                                  VECTORS. */
/*                        B.EQ.1    RETURN THE RIGHT SINGULAR VECTORS */
/*                                  IN V. */

/*     ON RETURN */

/*         S         DOUBLE PRECISION(MM), WHERE MM=MIN(N+1,P). */
/*                   THE FIRST MIN(N,P) ENTRIES OF S CONTAIN THE */
/*                   SINGULAR VALUES OF X ARRANGED IN DESCENDING */
/*                   ORDER OF MAGNITUDE. */

/*         E         DOUBLE PRECISION(P). */
/*                   E ORDINARILY CONTAINS ZEROS.  HOWEVER SEE THE */
/*                   DISCUSSION OF INFO FOR EXCEPTIONS. */

/*         U         DOUBLE PRECISION(LDU,K), WHERE LDU.GE.N.  IF */
/*                                   JOBA.EQ.1 THEN K.EQ.N, IF JOBA.GE.2 
*/
/*                                   THEN K.EQ.MIN(N,P). */
/*                   U CONTAINS THE MATRIX OF RIGHT SINGULAR VECTORS. */
/*                   U IS NOT REFERENCED IF JOBA.EQ.0.  IF N.LE.P */
/*                   OR IF JOBA.EQ.2, THEN U MAY BE IDENTIFIED WITH X */
/*                   IN THE SUBROUTINE CALL. */

/*         V         DOUBLE PRECISION(LDV,P), WHERE LDV.GE.P. */
/*                   V CONTAINS THE MATRIX OF RIGHT SINGULAR VECTORS. */
/*                   V IS NOT REFERENCED IF JOB.EQ.0.  IF P.LE.N, */
/*                   THEN V MAY BE IDENTIFIED WITH X IN THE */
/*                   SUBROUTINE CALL. */

/*         INFO      INTEGER. */
/*                   THE SINGULAR VALUES (AND THEIR CORRESPONDING */
/*                   SINGULAR VECTORS) S(INFO+1),S(INFO+2),...,S(M) */
/*                   ARE CORRECT (HERE M=MIN(N,P)).  THUS IF */
/*                   INFO.EQ.0, ALL THE SINGULAR VALUES AND THEIR */
/*                   VECTORS ARE CORRECT.  IN ANY EVENT, THE MATRIX */
/*                   B = TRANS(U)*X*V IS THE BIDIAGONAL MATRIX */
/*                   WITH THE ELEMENTS OF S ON ITS DIAGONAL AND THE */
/*                   ELEMENTS OF E ON ITS SUPER-DIAGONAL (TRANS(U) */
/*                   IS THE TRANSPOSE OF U).  THUS THE SINGULAR */
/*                   VALUES OF X AND B ARE THE SAME. */

/*     LINPACK. THIS VERSION DATED 03/19/79 . */
/*     G.W. STEWART, UNIVERSITY OF MARYLAND, ARGONNE NATIONAL LAB. */

/*     DSVDC USES THE FOLLOWING FUNCTIONS AND SUBPROGRAMS. */

/*     EXTERNAL DROT */
/*     BLAS DAXPY,DDOT,DSCAL,DSWAP,DNRM2,DROTG */
/*     FORTRAN DABS,DMAX1,MAX0,MIN0,MOD,DSQRT */

/*     INTERNAL VARIABLES */



/*     SET THE MAXIMUM NUMBER OF ITERATIONS. */

    /* Parameter adjustments */
    x_dim1 = *ldx;
    x_offset = x_dim1 + 1;
    x -= x_offset;
    --s;
    --e;
    u_dim1 = *ldu;
    u_offset = u_dim1 + 1;
    u -= u_offset;
    v_dim1 = *ldv;
    v_offset = v_dim1 + 1;
    v -= v_offset;
    --work;

    /* Function Body */
    maxit = 30;

/*     DETERMINE WHAT IS TO BE COMPUTED. */

    wantu = FALSE_;
    wantv = FALSE_;
    jobu = *job % 100 / 10;
    ncu = *n;
    if (jobu > 1) {
	ncu = min(*n,*p);
    }
    if (jobu != 0) {
	wantu = TRUE_;
    }
    if (*job % 10 != 0) {
	wantv = TRUE_;
    }

/*     REDUCE X TO BIDIAGONAL FORM, STORING THE DIAGONAL ELEMENTS */
/*     IN S AND THE SUPER-DIAGONAL ELEMENTS IN E. */

    *info = 0;
/* Computing MIN */
    i__1 = *n - 1;
    nct = min(i__1,*p);
/* Computing MAX */
/* Computing MIN */
    i__3 = *p - 2;
    i__1 = 0, i__2 = min(i__3,*n);
    nrt = max(i__1,i__2);
    lu = max(nct,nrt);
    if (lu < 1) {
	goto L170;
    }
    i__1 = lu;
    for (l = 1; l <= i__1; ++l) {
	lp1 = l + 1;
	if (l > nct) {
	    goto L20;
	}

/*           COMPUTE THE TRANSFORMATION FOR THE L-TH COLUMN AND */
/*           PLACE THE L-TH DIAGONAL IN S(L). */

	i__2 = *n - l + 1;
	s[l] = dnrm2_(&i__2, &x[l + l * x_dim1], &c__1);
	if (s[l] == 0.) {
	    goto L10;
	}
	if (x[l + l * x_dim1] != 0.) {
	    s[l] = d_sign(&s[l], &x[l + l * x_dim1]);
	}
	i__2 = *n - l + 1;
	d__1 = 1. / s[l];
	dscal_(&i__2, &d__1, &x[l + l * x_dim1], &c__1);
	x[l + l * x_dim1] += 1.;
L10:
	s[l] = -s[l];
L20:
	if (*p < lp1) {
	    goto L50;
	}
	i__2 = *p;
	for (j = lp1; j <= i__2; ++j) {
	    if (l > nct) {
		goto L30;
	    }
	    if (s[l] == 0.) {
		goto L30;
	    }

/*              APPLY THE TRANSFORMATION. */

	    i__3 = *n - l + 1;
	    t = -ddot_(&i__3, &x[l + l * x_dim1], &c__1, &x[l + j * x_dim1], &
		    c__1) / x[l + l * x_dim1];
	    i__3 = *n - l + 1;
	    daxpy_(&i__3, &t, &x[l + l * x_dim1], &c__1, &x[l + j * x_dim1], &
		    c__1);
L30:

/*           PLACE THE L-TH ROW OF X INTO  E FOR THE */
/*           SUBSEQUENT CALCULATION OF THE ROW TRANSFORMATION. */

	    e[j] = x[l + j * x_dim1];
/* L40: */
	}
L50:
	if (! wantu || l > nct) {
	    goto L70;
	}

/*           PLACE THE TRANSFORMATION IN U FOR SUBSEQUENT BACK */
/*           MULTIPLICATION. */

	i__2 = *n;
	for (i = l; i <= i__2; ++i) {
	    u[i + l * u_dim1] = x[i + l * x_dim1];
/* L60: */
	}
L70:
	if (l > nrt) {
	    goto L150;
	}

/*           COMPUTE THE L-TH ROW TRANSFORMATION AND PLACE THE */
/*           L-TH SUPER-DIAGONAL IN E(L). */

	i__2 = *p - l;
	e[l] = dnrm2_(&i__2, &e[lp1], &c__1);
	if (e[l] == 0.) {
	    goto L80;
	}
	if (e[lp1] != 0.) {
	    e[l] = d_sign(&e[l], &e[lp1]);
	}
	i__2 = *p - l;
	d__1 = 1. / e[l];
	dscal_(&i__2, &d__1, &e[lp1], &c__1);
	e[lp1] += 1.;
L80:
	e[l] = -e[l];
	if (lp1 > *n || e[l] == 0.) {
	    goto L120;
	}

/*              APPLY THE TRANSFORMATION. */

	i__2 = *n;
	for (i = lp1; i <= i__2; ++i) {
	    work[i] = 0.;
/* L90: */
	}
	i__2 = *p;
	for (j = lp1; j <= i__2; ++j) {
	    i__3 = *n - l;
	    daxpy_(&i__3, &e[j], &x[lp1 + j * x_dim1], &c__1, &work[lp1], &
		    c__1);
/* L100: */
	}
	i__2 = *p;
	for (j = lp1; j <= i__2; ++j) {
	    i__3 = *n - l;
	    d__1 = -e[j] / e[lp1];
	    daxpy_(&i__3, &d__1, &work[lp1], &c__1, &x[lp1 + j * x_dim1], &
		    c__1);
/* L110: */
	}
L120:
	if (! wantv) {
	    goto L140;
	}

/*              PLACE THE TRANSFORMATION IN V FOR SUBSEQUENT */
/*              BACK MULTIPLICATION. */

	i__2 = *p;
	for (i = lp1; i <= i__2; ++i) {
	    v[i + l * v_dim1] = e[i];
/* L130: */
	}
L140:
L150:
/* L160: */
	;
    }
L170:

/*     SET UP THE FINAL BIDIAGONAL MATRIX OR ORDER M. */

/* Computing MIN */
    i__1 = *p, i__2 = *n + 1;
    m = min(i__1,i__2);
    nctp1 = nct + 1;
    nrtp1 = nrt + 1;
    if (nct < *p) {
	s[nctp1] = x[nctp1 + nctp1 * x_dim1];
    }
    if (*n < m) {
	s[m] = 0.;
    }
    if (nrtp1 < m) {
	e[nrtp1] = x[nrtp1 + m * x_dim1];
    }
    e[m] = 0.;