cmatrix.c

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

        c[j+1][i]=t01;
        }
      else {
        c[i][j]=t00;
        c[i][j+1]=t01;
        }
      }
    }
  if(odd_nc) {
    for(i=0; i < nr-1 ; i+=2) {
      att=a[i]; bt=b[j];
      at1=a[i+1];
      if(add) {
        if(tc) {
          t00 = c[j][i];
          t10 = c[j][i+1];
          }
        else {
          t00 = c[i][j];
          t10 = c[i+1][j];
          }
        }
      else t00=t10=0.0;
      for(k= nl; k ; k--,att++,bt++,at1++) {
        t00 += *att * *bt;
        t10 += *at1 * *bt;
        }
      if(tc) {
        c[j][i]=t00;
        c[j][i+1]=t10;
        }
      else {
        c[i][j]=t00;
        c[i+1][j]=t10;
        }
      }
    if(odd_nr) {
      att=a[i]; bt=b[j];
      if(add) t00 = (tc) ? c[j][i] : c[i][j];
      else t00=0.0;
      for(k=nl; k ; k--,att++,bt++) t00 += *att * *bt;
      if(tc) c[j][i]=t00;
      else c[i][j]=t00;
      }
    }

  if(ta) {
      cmat_transpose_matrix(a,nr,nl);
      if (old_a) {
          free(a);
          a = old_a;
        }
      cmat_matrix_pointers(a,a[0],nr,nl);
    }
  if(!tb) {
      cmat_transpose_matrix(b,nc,nl);
      if (old_b) {
          free(b);
          b = old_b;
        }
      cmat_matrix_pointers(b,b[0],nl,nc);
    }
  }

/*
 * a is symmetric (na,na) in a triangular storage format
 * b is rectangular (na,nb)
 * a (+)= b * transpose(b) (+= if add)
 */
void
cmat_symmetric_mxm(double**a,int na, /* a is (na,na) */
                   double**b,int nb, /* b is (na,nb) */
                   int add)
{
  int i,j,k;
  for (i=0; i<na; i++) {
      double*ai=a[i];
      for (j=0; j<=i; j++) {
          double*bi=b[i];
          double*bj=b[j];
          double tmp;
          if (add) tmp = 0.0;
          else tmp = ai[j];
          for (k=nb; k; k--,bi++,bj++) {
              tmp += *bi * *bj;
            }
          ai[j] = tmp;
        }
    }
}

/*
 * a is symmetric (na,na) in a triangular storage format
 * b is symmetric (nb,nb) in a triangular storage format
 * a (+)= c * b * transpose(c) (+= if add)
 */
void
cmat_transform_symmetric_matrix(double**a,int na, /* a is (na,na) */
                                double**b,int nb, /* b is (nb,nb) */
                                double**c,        /* c is (na,nb) */
                                int add)
{
  int i,j,k;
  double**t;
  double* brow;

  /* create a temporary matrix, t */
  t = cmat_new_rect_matrix(na,nb);

  /* t = transpose(b * transpose(c)) */
  brow = (double*) malloc(sizeof(double)*nb);
  if (!brow) {
    fprintf(stderr,"cmat_transform_symmetric_matrix: malloc brow failed\n");
    abort();
  }
  for (i=0; i<nb; i++) {
      for (k=0; k<=i; k++) brow[k] = b[i][k];
      for (   ; k<nb; k++) brow[k] = b[k][i];
      for (j=0; j<na; j++) {
          double*bi = brow;
          double*cj = c[j];
          double tmp = 0.0;
          for (k=nb; k; k--,bi++,cj++) tmp += *bi * *cj;
          t[j][i] = tmp;
        }
    }
  free(brow);

  /* a = c * transpose(t) */
  for (i=0; i<na; i++) {
      for (j=0; j<=i; j++) {
          double*ci = c[i];
          double*tj = t[j];
          double tmp;
          if (add) tmp = a[i][j];
          else tmp = 0.0;
          for (k=nb; k; k--,ci++,tj++) tmp += *ci * *tj;
          a[i][j] = tmp;
        }
    }

  /* delete the temporary */
  cmat_delete_matrix(t);
}

/*
 * a is symmetric (na,na) in a triangular storage format
 * b is diagonal (nb,nb) in a vector storage format
 * a (+)= c * b * transpose(c) (+= if add)
 */
void
cmat_transform_diagonal_matrix(double**a,int na, /* a is (na,na) */
                               double*b,int nb,  /* b is (nb,nb) */
                               double**c,        /* c is (na,nb) */
                               int add)
{
  int i,j,k;
  double t;

  for (i=0; i < na; i++) {
    for (j=0; j <= i; j++) {
      t=0;
      for (k=0; k < nb; k++)
        t += c[i][k] * c[j][k] * b[k];
      if (add)
        a[i][j] += t;
      else
        a[i][j] = t;
    }
  }
}

/*
 * Argument a contains pointers to the rows of a symmetrix matrix.  The
 * in each row is the row number + 1.  These rows are stored in
 * contiguous memory starting with 0.  Evecs also contains pointers to
 * contiguous memory.  N is the dimension.
 */
void
cmat_diag(double**a, double*evals, double**evecs, int n,
              int matz, double tol)
{
  int i,j;
  int diagonal=1;
  double*fv1;

 /* I'm having problems with diagonalizing matrices which are already
  * diagonal.  So let's first check to see if _a_ is diagonal, and if it
  * is, then just return the diagonal elements in evals and a unit matrix
  * in evecs
  */

  for (i=1; i < n; i++) {
    for (j=0; j < i; j++) {
      if (fabs(a[i][j]) > tol) diagonal=0;
      }
    }

  if (diagonal) {
    for(i=0; i < n; i++) {
      evals[i] = a[i][i];
      evecs[i][i] = 1.0;

      for(j=0; j < i; j++) {
        evecs[i][j] = evecs[j][i] = 0.0;
        }
      }
    eigsort(n,evals,evecs);
    return;
    }

  fv1 = (double*) malloc(sizeof(double)*n);
  if (!fv1) {
    fprintf(stderr,"cmat_diag: malloc fv1 failed\n");
    abort();
  }

  for(i=0; i < n; i++) {
      for(j=0; j <= i; j++) {
          evecs[i][j] = evecs[j][i] = a[i][j];
        }
    }

  tred2(n,evecs,evals,fv1,1);

  cmat_transpose_square_matrix(evecs,n);
  tqli(n,evals,evecs,fv1,1,tol);
  cmat_transpose_square_matrix(evecs,n);

  eigsort(n,evals,evecs);

  free(fv1);
  }

#define dsign(a,b) (((b) >= 0.0) ? fabs(a) : -fabs(a))

static void
tred2(int n,double** a,double* d,double* e,int matz)
{
  int i,j,k,l;
  double f,g,h,hh,scale,scale_inv,h_inv;
  if (n == 1) return;

  for(i=n-1; i > 0; i--) {
    l = i-1;
    h = 0.0;
    scale = 0.0;
    if(l) {
      for(k=0; k <= l; k++) scale += fabs(a[i][k]);
      if (scale == 0.0) e[i] = a[i][l];
      else {
        scale_inv=1.0/scale;
        for (k=0; k <= l; k++) {
          a[i][k] *= scale_inv;
          h += a[i][k]*a[i][k];
          }
        f=a[i][l];
        g= -(dsign(sqrt(h),f));
        e[i] = scale*g;
        h -= f*g;
        a[i][l] = f-g;
        f = 0.0;
        h_inv=1.0/h;
        for (j=0; j <= l; j++) {
          if (matz) a[j][i] = a[i][j]*h_inv;
          g = 0.0;
          for (k=0; k <= j; k++) g += a[j][k]*a[i][k];
          if (l > j) for (k=j+1; k <= l; k++) g += a[k][j]*a[i][k];
          e[j] = g*h_inv;
          f += e[j]*a[i][j];
          }
        hh = f/(h+h);
        for (j=0; j <= l; j++) {
          f = a[i][j];
          g = e[j] - hh*f;
          e[j] = g;
          for (k=0; k <= j; k++) a[j][k] -= (f*e[k] + g*a[i][k]);
          }
        }
      }
    else {
      e[i] = a[i][l];
      }
    d[i] = h;
    }
  if(matz) d[0] = 0.0;
  e[0] = 0.0;

  for(i=0; i < n; i++) {
    l = i-1;
    if (matz) {
      if(d[i]) {
        for(j=0; j <= l; j++) {
          g = 0.0;
          for(k=0; k <= l; k++) g += a[i][k]*a[k][j];
          for(k=0; k <= l; k++) a[k][j] -= g*a[k][i];
          }
        }
      }
    d[i] = a[i][i];
    if(matz) {
      a[i][i] = 1.0;
      if(l >= 0) for (j=0; j<= l; j++) a[i][j] = a[j][i] = 0.0;
      }
    }
  }

static void
tqli(int n, double* d, double** z, double* e, int matz, double toler)
{
  register int k;
  int i,l,m,iter;
  double g,r,s,c,p,f,b;
  double azi;

  f=0.0;
  if (n == 1) {
    d[0]=z[0][0];
    z[0][0] = 1.0;
    return;
    }

  for (i=1; i < n ; i++) e[i-1] = e[i];
  e[n-1] = 0.0;
  for (l=0; l < n; l++) {
    iter = 0;
L1:
    for (m=l; m < n-1;m++) if (fabs(e[m]) < toler) goto L2;
    m=n-1;
L2:
    if (m != l) {
      if (iter++ == 30) {
        fprintf (stderr,"tqli not converging %d %g\n",l,e[l]);
        continue;
        }

      g = (d[l+1]-d[l])/(2.0*e[l]);
      r = sqrt(g*g + 1.0);
      g = d[m] - d[l] + e[l]/((g + dsign(r,g)));
      s=1.0;
      c=1.0;
      p=0.0;
      for (i=m-1; i >= l; i--) {
        f = s*e[i];
        b = c*e[i];
        if (fabs(f) >= fabs(g)) {
          c = g/f;
          r = sqrt(c*c + 1.0);
          e[i+1] = f*r;
          s=1.0/r;
          c *= s;
          }
        else {
          s = f/g;
          r = sqrt(s*s + 1.0);
          e[i+1] = g*r;
          c = 1.0/r;
          s *= c;
          }
        g = d[i+1] - p;
        r = (d[i]-g)*s + 2.0*c*b;
        p = s*r;
        d[i+1] = g+p;
        g = c*r-b;

        if (matz) {
          double *zi = z[i];
          double *zi1 = z[i+1];
          for (k=n; k ; k--,zi++,zi1++) {
            azi = *zi;
            f = *zi1;
            *zi1 = azi*s + c*f;
            *zi = azi*c - s*f;
            }
          }
        }

      d[l] -= p;
      e[l] = g;
      e[m] = 0.0;
      goto L1;
      }
    }
  }

static void
eigsort(int n, double* d, double** v)
{
  int i,j,k;
  double p;

  for(i=0; i < n-1 ; i++) {
    k=i;
    p=d[i];
    for(j=i+1; j < n; j++) {
      if(d[j] < p) {
        k=j;
        p=d[j];
        }
      }
    if(k != i) {
      d[k]=d[i];
      d[i]=p;
      for(j=0; j < n; j++) {
        p=v[j][i];
        v[j][i]=v[j][k];
        v[j][k]=p;
        }
      }
    }
  }

void
cmat_schmidt(double **C, double *S, int nrow, int nc)
{
  int i,j,ij;
  int m;
  double vtmp;
  double *v = (double*) malloc(sizeof(double)*nrow);

  if (!v) {
    fprintf(stderr,"cmat_schmidt: could not malloc v(%d)\n",nrow);
    abort();
  }
  
  for (m=0; m < nc; m++) {
    v[0] = C[0][m] * S[0];

    for (i=ij=1; i < nrow; i++) {
      for (j=0,vtmp=0.0; j < i; j++,ij++) {
        vtmp += C[j][m]*S[ij];
        v[j] += C[i][m]*S[ij];
      }
      v[i] = vtmp + C[i][m]*S[ij];
      ij++;
    }

    for (i=0,vtmp=0.0; i < nrow; i++)
      vtmp += v[i]*C[i][m];

    if (!vtmp) {
      fprintf(stderr,"cmat_schmidt: bogus\n");
      abort();
    }

    if (vtmp < 1.0e-15)
      vtmp = 1.0e-15;

    vtmp = 1.0/sqrt(vtmp);
    
    for (i=0; i < nrow; i++) {
      v[i] *= vtmp;
      C[i][m] *= vtmp;
    }

    if (m < nc-1) {
      for (i=m+1,vtmp=0.0; i < nc; i++) {
        for (j=0,vtmp=0.0; j < nrow; j++)
          vtmp += v[j] * C[j][i];
        for (j=0; j < nrow; j++)
          C[j][i] -= vtmp * C[j][m];
      }
    }
  }
}

/* Returns the number of linearly independent vectors
   orthogonal wrt S. */
int
cmat_schmidt_tol(double **C, double *S, int nrow, int ncol,
                 double tolerance, double *res)
{
  int i,j,ij;
  int m;
  double vtmp;
  int northog = 0;
  double *v = (double*) malloc(sizeof(double)*nrow);

  if (res) *res = 1.0;

  if (!v) {
    fprintf(stderr,"cmat_schmidt_tol: could not malloc v(%d)\n",nrow);
    abort();
  }

  /* Orthonormalize the columns of C wrt S. */
  for (m=0; m < ncol; m++) {
    v[0] = C[0][m] * S[0];

    for (i=ij=1; i < nrow; i++) {
      for (j=0,vtmp=0.0; j < i; j++,ij++) {
        vtmp += C[j][m]*S[ij];
        v[j] += C[i][m]*S[ij];
      }
      v[i] = vtmp + C[i][m]*S[ij];
      ij++;
    }

    for (i=0,vtmp=0.0; i < nrow; i++)
      vtmp += v[i]*C[i][m];

    if (vtmp < tolerance) continue;

    if (res && (m == 0 || vtmp < *res)) *res = vtmp;

    vtmp = 1.0/sqrt(vtmp);
    
    for (i=0; i < nrow; i++) {
      v[i] *= vtmp;
      C[i][northog] = C[i][m] * vtmp;
    }

    for (i=m+1,vtmp=0.0; i < ncol; i++) {
        for (j=0,vtmp=0.0; j < nrow; j++)
            vtmp += v[j] * C[j][i];
        for (j=0; j < nrow; j++)
            C[j][i] -= vtmp * C[j][northog];
      }
    northog++;
  }
  return northog;
}