svd.c

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

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

typedef double real;
typedef int integer;

#define r_sign(a,b) ((*b<0.0)?-fabs(*a):fabs(*a))
#define min(a,b) (((a)<(b))?(a):(b))
#define max(a,b) (((a)>(b))?(a):(b))
#define dmax(a,b) (((a)>(b))?(a):(b))

int
sing_(double *q, int *lq, int *iq, double *s, double *p,
     int *lp, int *ip, double *a, int *la, int *m, int *n, double *w);

static int
bidag2_(double *d, double *b, double *q, int *lq, int *iq,
        double *p, int *lp, int *ip, double *a, int *la, int *m, int *n);
static int
hsr1_(double *a, int *la, int *n);
static int
hsr2_(double *a, int *la, int *n);
static int
hsr3_(double *a, int *la, int *m, int *n);
static int
hsr4_(double *a, int *la, int *m, int *n);
static int
hsr5_(double *a, int *la, int *m, int *n);
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);
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);
static int
sft_(double *s, double *a, double *b, double *c,
     double *d, double *e2, double *e1, double *e0, double *f2, double *f1);
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);
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);
static int
eig3_(double *ea, double *eb, double *a, double *b, double *y, double *z);
static int
sort2_(double *x, double *y, double *w, int *n);
static int
fgv_(double *x, double *y, double *s, double *p, double *q,
     double *a, double *b);

#ifdef TEST

main(int argc, char**argv)
{
  int i,j;
  double *tmp;

  int m = atoi(argv[1]);
  int n = atoi(argv[2]);

  int l = ((m<n)?m:n);
  
  double *q = (double*)malloc(sizeof(double)*m*m);
  double *s = (double*)malloc(sizeof(double)*n);
  double *p = (double*)malloc(sizeof(double)*n*n);
  double *a = (double*)malloc(sizeof(double)*n*m);
  double *w = (double*)malloc(sizeof(double)*3*m);

  int lq = m;
  int lp = n;
  int iq = 3;
  int ip = 3;
  int la = lq;

  tmp = a;
  for (i=0; i<m; i++) {
      for (j=0; j<n; j++) {
          *tmp++ = drand48() * ((drand48()<0.5)?1.0:-1.0);
          /* *tmp++ = 1.0; */
        }
    }

  printf("A:\n");
  for (i=0; i<m; i++) {
      for (j=0; j<n; j++) {
          printf(" % 6.4f", a[j*m + i]);
        }
      printf("\n");
    }

  sing_(q, &lq, &iq, s, p,
        &lp, &ip, a, &la, &m, &n, w);

  printf("Q:\n");
  for (i=0; i<m; i++) {
      for (j=0; j<m; j++) {
          printf(" % 6.4f",q[j*m+i]);
        }
      printf("\n");
    }

  printf("P:\n");
  for (i=0; i<n; i++) {
      for (j=0; j<n; j++) {
          printf(" % 6.4f",p[j*n+i]);
        }
      printf("\n");
    }

  printf("S:\n");
  tmp = s;
  for (j=0; j<l; j++) {
      printf(" % 6.4f",*tmp++);
    }
  printf("\n");

  printf("QQt:\n");
  for (i=0; i<m; i++) {
      for (j=0; j<m; j++) {
          int k;
          double tmp = 0.0;
          for (k=0; k<m; k++) {
              tmp += q[k*m+i]*q[k*m+j];
            }
          printf(" % 6.4f",tmp);
        }
      printf("\n");
    }

  printf("QtQ:\n");
  for (i=0; i<m; i++) {
      for (j=0; j<m; j++) {
          int k;
          double tmp = 0.0;
          for (k=0; k<m; k++) {
              tmp += q[i*m+k]*q[j*m+k];
            }
          printf(" % 6.4f",tmp);
        }
      printf("\n");
    }

  printf("PPt:\n");
  for (i=0; i<n; i++) {
      for (j=0; j<n; j++) {
          int k;
          double tmp = 0.0;
          for (k=0; k<l; k++) {
              tmp += p[k*n+i]*p[k*n+j];
            }
          printf(" % 6.4f",tmp);
        }
      printf("\n");
    }

  printf("PtP:\n");
  for (i=0; i<n; i++) {
      for (j=0; j<n; j++) {
          int k;
          double tmp = 0.0;
          for (k=0; k<l; k++) {
              tmp += p[i*n+k]*p[j*n+k];
            }
          printf(" % 6.4f",tmp);
        }
      printf("\n");
    }

  printf("QSPt:\n");
  for (i=0; i<m; i++) {
      for (j=0; j<n; j++) {
          int k;
          double tmp = 0.0;
          for (k=0; k<l; k++) {
              tmp += q[k*m+i]*s[k]*p[k*n+j];
            }
          printf(" % 6.4f",tmp);
        }
      printf("\n");
    }
  

  return 0;
}

#endif /* TEST */

/*      ________________________________________________________ */
/*     |                                                        | */
/*     |COMPUTE SINGULAR VALUE DECOMPOSITION OF A GENERAL MATRIX| */
/*     |      A = Q TIMES DIAGONAL MATRIX TIMES P TRANSPOSE     | */
/*     |                                                        | */
/*     |    INPUT:                                              | */
/*     |                                                        | */
/*     |         S     --ARRAY WITH AT LEAST N ELEMENTS         | */
/*     |                                                        | */
/*     |         LQ    --LEADING (ROW) DIMENSION OF ARRAY Q     | */
/*     |                                                        | */
/*     |         IQ    --AN INTEGER WHICH INDICATES WHICH COL-  | */
/*     |                 UMNS OF Q TO COMPUTE (= 0 MEANS NONE,  | */
/*     |                 = 1 MEANS FIRST L, = 2 MEANS LAST M-L, | */
/*     |                 = 3 MEANS ALL M WHERE L = MIN(M,N))    | */
/*     |                                                        | */
/*     |         LP    --LEADING (ROW) DIMENSION OF ARRAY P     | */
/*     |                                                        | */
/*     |         IP    --AN INTEGER (LIKE IQ) WHICH INDICATES   | */
/*     |                 WHICH COLUMNS OF P TO COMPUTE          | */
/*     |                                                        | */
/*     |         A     --ARRAY CONTAINING COEFFICIENT MATRIX    | */
/*     |                                                        | */
/*     |         LA    --LEADING (ROW) DIMENSION OF ARRAY A     | */
/*     |                                                        | */
/*     |         M     --ROW DIMENSION OF MATRIX STORED IN A    | */
/*     |                                                        | */
/*     |         N     --COLUMN DIMENSION OF MATRIX STORED IN A | */
/*     |                                                        | */
/*     |         W     --WORK ARRAY(LENGTH AT LEAST MAX(M,3L-1))| */
/*     |                                                        | */
/*     |    OUTPUT:                                             | */
/*     |                                                        | */
/*     |         Q     --Q FACTOR IN THE SINGULAR VALUE DECOMP. | */
/*     |                                                        | */
/*     |         S     --SINGULAR VALUES IN DECREASING ORDER    | */
/*     |                                                        | */
/*     |         P     --P FACTOR IN THE SINGULAR VALUE DECOMP. | */
/*     |                                                        | */
/*     |         A     --THE HOUSEHOLDER VECTORS USED IN THE    | */
/*     |                 REDUCTION PROCESS                      | */
/*     |                                                        | */
/*     |    NOTE:                                               | */
/*     |                                                        | */
/*     |         EITHER P OR Q CAN BE IDENTIFIED WITH A BUT NOT | */
/*     |         BOTH. WHEN EITHER P OR Q ARE IDENTIFIED WITH A,| */
/*     |         THE HOUSEHOLDER VECTORS IN A ARE DESTROYED. IF | */
/*     |         IQ = 2, Q MUST HAVE M COLUMNS EVEN THOUGH THE  | */
/*     |         OUTPUT ARRAY HAS JUST M-L COLUMNS. SIMILARLY IF| */
/*     |         IP = 2, P MUST HAVE N COLUMNS EVEN THOUGH THE  | */
/*     |         OUTPUT ARRAY HAS JUST N-L COLUMNS.             | */
/*     |                                                        | */
/*     |    BUILTIN FUNCTIONS: MIN0                             | */
/*     |    PACKAGE SUBROUTINES: BIDAG2,EIG3,FGV,HSR1-HSR5,SCL, | */
/*     |                         SFT,SINGB,SNG0,SNG1,SORT2      | */
/*     |________________________________________________________| */

int
sing_(real *q, int *lq, int *iq, double *s, double *p,
      int *lp, int *ip, double *a, int *la, int *m, int *n, double *w)
{
    /* System generated locals */
    integer a_dim1, a_offset, q_dim1, q_offset, p_dim1, p_offset, i__1;

    /* Local variables */
    static integer i, l;
    static integer jl, iu, jr;

    /* Parameter adjustments */
    --w;
    a_dim1 = *la;
    a_offset = a_dim1 + 1;
    a -= a_offset;
    p_dim1 = *lp;
    p_offset = p_dim1 + 1;
    p -= p_offset;
    --s;
    q_dim1 = *lq;
    q_offset = q_dim1 + 1;
    q -= q_offset;

    /* Function Body */
    if (*iq >= 0) {
	goto L20;
    }
L10:
    fprintf(stderr,"ERROR: INPUT PARAMETER IQ FOR SUBROUTINE SING\n");
    fprintf(stderr,"EITHER LESS THAN 0 OR GREATER THAN 3\n");
    abort();
L20:
    if (*iq > 3) {
	goto L10;
    }
    jl = 0;
    if (*iq == 0) {
	goto L30;
    }
    if (*iq == 2) {
	goto L30;
    }
    jl = 1;
L30:
    if (*ip >= 0) {
	goto L50;
    }
L40:
    fprintf(stderr,"ERROR: INPUT PARAMETER IP FOR SUBROUTINE SING\n");
    fprintf(stderr,"EITHER LESS THAN 0 OR GREATER THAN 3\n");
    abort();
L50:
    if (*ip > 3) {
	goto L40;
    }
    jr = 0;
    if (*ip == 0) {
	goto L60;
    }
    if (*ip == 2) {
	goto L60;
    }
    jr = 1;
L60:
    bidag2_(&s[1], &w[1], &q[q_offset], lq, iq, &p[p_offset], lp, ip, &a[
	    a_offset], la, m, n);
    l = min(*m,*n);
    if (l > 1) {
	goto L80;
    }
    if (s[1] >= (double)0.) {
	return 0;
    }
    s[1] = -(double)s[1];
    if ((real) jl == (double)0.) {
	return 0;
    }
    i__1 = *m;
    for (i = 1; i <= i__1; ++i) {
/* L70: */
	q[i + q_dim1] = -(double)q[i + q_dim1];
    }
    return 0;
L80:
    iu = 0;
    if (*m >= *n) {
	goto L90;
    }
    iu = 1;
L90:
    singb_(&s[1], &l, &w[1], &iu, &q[q_offset], lq, m, &jl, &p[p_offset], lp, 
	    n, &jr, &w[l], &w[l + l]);
    return 0;
} /* sing_ */


/*      ________________________________________________________ */
/*     |                                                        | */
/*     |      REDUCE A GENERAL MATRIX A TO BIDIAGONAL FORM      | */
/*     |     A = Q TIMES BIDIAGONAL MATRIX TIMES P TRANSPOSE    | */
/*     |                                                        | */
/*     |    INPUT:                                              | */
/*     |                                                        | */
/*     |         D     --ARRAY WITH AT LEAST N ELEMENTS         | */
/*     |                                                        | */
/*     |         B     --ARRAY WITH AT LEAST M ELEMENTS         | */
/*     |                                                        | */
/*     |         LQ    --LEADING (ROW) DIMENSION OF ARRAY Q     | */
/*     |                                                        | */
/*     |         IQ    --AN INTEGER WHICH INDICATES WHICH COL-  | */
/*     |                 UMNS OF Q TO COMPUTE (= 0 MEANS NONE,  | */
/*     |                 = 1 MEANS FIRST L, = 2 MEANS LAST M-L, | */
/*     |                 = 3 MEANS ALL M WHERE L = MIN(M,N))    | */
/*     |                                                        | */
/*     |         LP    --LEADING (ROW) DIMENSION OF ARRAY P     | */
/*     |                                                        | */
/*     |         IP    --AN INTEGER (LIKE IQ) WHICH INDICATES   | */
/*     |                 WHICH COLUMNS OF P TO COMPUTE          | */
/*     |                                                        | */
/*     |         A     --ARRAY CONTAINING COEFFICIENT MATRIX    | */
/*     |                                                        | */
/*     |         LA    --LEADING (ROW) DIMENSION OF ARRAY A     | */
/*     |                                                        | */
/*     |         M     --ROW DIMENSION OF MATRIX STORED IN A    | */
/*     |                                                        | */
/*     |         N     --COLUMN DIMENSION OF MATRIX STORED IN A | */
/*     |                                                        | */
/*     |    OUTPUT:                                             | */
/*     |                                                        | */
/*     |         D     --DIAGONAL OF BIDIAGONAL FORM            | */
/*     |                                                        | */
/*     |         B     --SUPERDIAGONAL (IF M .GE. N) OR SUBDIAG-| */
/*     |                 ONAL (IF M .LT. N) OF BIDIAGONAL FORM  | */
/*     |                                                        | */
/*     |         A     --THE HOUSEHOLDER VECTORS USED IN THE    | */
/*     |                 REDUCTION PROCESS                      | */
/*     |                                                        | */
/*     |         Q     --THE Q FACTOR IN THE BIDIAGONALIZATION  | */
/*     |                                                        | */
/*     |         P     --THE P FACTOR IN THE BIDIAGONALIZATION  | */
/*     |                                                        | */
/*     |    NOTE:                                               | */
/*     |                                                        | */
/*     |         EITHER P OR Q CAN BE IDENTIFIED WITH A BUT NOT | */
/*     |         BOTH. WHEN EITHER P OR Q ARE IDENTIFIED WITH A,| */
/*     |         THEN THE HOUSEHOLDER VECTORS IN A ARE DESTROYED| */
/*     |                                                        | */
/*     |    BUILTIN FUNCTIONS: ABS,MIN0,SQRT                    | */
/*     |    PACKAGE SUBROUTINES: HSR1-HSR5                      | */
/*     |________________________________________________________| */

static int
bidag2_(double *d, double *b, double *q, int *lq, int *iq,
        double *p, int *lp, int *ip, double *a, int *la, int *m, int *n)
{
    /* System generated locals */
    integer a_dim1, a_offset, q_dim1, q_offset, p_dim1, p_offset, i__1, i__2;
    real r__1;

    /* Local variables */
    static integer h, i, j, k, l;
    static real r, s, t, u;
    static integer jp, jq;

    /* Parameter adjustments */
    a_dim1 = *la;
    a_offset = a_dim1 + 1;
    a -= a_offset;
    p_dim1 = *lp;
    p_offset = p_dim1 + 1;
    p -= p_offset;
    q_dim1 = *lq;
    q_offset = q_dim1 + 1;
    q -= q_offset;
    --b;
    --d;

    /* Function Body */
    l = min(*m,*n);
    if (*iq >= 0) {
	goto L20;
    }
L10:
    fprintf(stderr,"ERROR: INPUT PARAMETER IQ FOR SUBROUTINE BIDAG2\n");
    fprintf(stderr,"EITHER LESS THAN 0 OR GREATER THAN 3\n");
    abort();
L20:
    if (*iq > 3) {
	goto L10;
    }
    jq = *iq;
    if (*iq <= 1) {
	goto L30;
    }
    if (*iq == 3) {
	goto L30;
    }
    if (*m == l) {
	jq = 0;
    }
L30:
    if (*ip >= 0) {
	goto L50;
    }
L40:
    fprintf(stderr,"ERROR: INPUT PARAMETER IP FOR SUBROUTINE BIDAG2\n");
    fprintf(stderr,"EITHER LESS THAN 0 OR GREATER THAN 3\n");
    abort();
L50:
    if (*ip > 3) {
	goto L40;
    }
    jp = *ip;
    if (*ip <= 1) {
	goto L60;
    }
    if (*ip == 3) {
	goto L60;
    }
    if (*n == l) {
	jp = 0;
    }
L60:
    k = 1;
    h = 2;
    if (*m < *n) {
	goto L330;
    }
    if (*m > 1) {
	goto L70;
    }
    d[1] = a[a_dim1 + 1];
    if (*iq > 0) {
	q[q_dim1 + 1] = (double)1.;
    }
    if (*ip > 0) {
	p[p_dim1 + 1] = (double)1.;
    }
    return 0;
L70:
    j = k;
    k = h;
    i__1 = *m;
    for (i = k; i <= i__1; ++i) {
/* L80: */
	if (a[i + j * a_dim1] != (double)0.) {
	    goto L110;
	}
    }
    d[j] = a[j + j * a_dim1];
    a[j + j * a_dim1] = (double)0.;
    i__1 = *n;
    for (i = k; i <= i__1; ++i) {
/* L90: */
	d[i] = a[j + i * a_dim1];
    }
    if (jq == 0) {
	goto L200;
    }
    i__1 = *m;
    for (i = j; i <= i__1; ++i) {
/* L100: */
	q[i + j * q_dim1] = (double)0.;
    }
    goto L200;
L110:
    t = (r__1 = a[j + 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 L120;
	}
	u = (double)1. / s;
/* Computing 2nd power */
	r__1 = t * u;
	r = r * (r__1 * r__1) + (double)1.;
	t = s;
	goto L130;
L120:
/* Computing 2nd power */
	r__1 = s * u;
	r += r__1 * r__1;
L130:
	;
    }
    s = t * sqrt(r);
    r = a[j + j * a_dim1];
    u = (double)1. / sqrt(s * (s + fabs(r)));
    if (r < (double)0.) {
	s = -(double)s;
    }
    d[j] = -(double)s;
    a[j + j * a_dim1] = u * (r + s);
    i__1 = *m;
    for (i = k; i <= i__1; ++i) {
/* L140: */
	a[i + j * a_dim1] *= u;
    }
    if (jq == 0) {
	goto L160;
    }
    i__1 = *m;
    for (i = j; i <= i__1; ++i) {
/* L150: */
	q[i + j * q_dim1] = a[i + j * a_dim1];
    }
L160:
    if (k > *n) {
	goto L620;
    }
    i__1 = *n;
    for (l = k; l <= i__1; ++l) {
	t = (double)0.;
	i__2 = *m;
	for (i = j; i <= i__2; ++i) {
/* L170: */
	    t += a[i + j * a_dim1] * a[i + l * a_dim1];
	}
	a[j + l * a_dim1] -= t * a[j + j * a_dim1];
	d[l] = a[j + l * a_dim1];
	i__2 = *m;
	for (i = k; i <= i__2; ++i) {
/* L180: */
	    a[i + l * a_dim1] -= t * a[i + j * a_dim1];
	}
/* L190: */
    }
L200:
    h = k + 1;
    if (k < *n) {
	goto L210;
    }
    if (k > *n) {
	goto L620;
    }
    if (*m == *n) {
	goto L610;
    }
    b[j] = a[j + *n * a_dim1];
    goto L70;
L210:
    i__1 = *n;
    for (i = h; i <= i__1; ++i) {
/* L220: */