blkdirect.c

很经典的电磁计算（电容计算）软件 MIT90年代开发

/*
Copyright (c) 1990 Massachusetts Institute of Technology, Cambridge, MA.
All rights reserved.

This Agreement gives you, the LICENSEE, certain rights and obligations.
By using the software, you indicate that you have read, understood, and
will comply with the terms.

Permission to use, copy and modify for internal, noncommercial purposes
is hereby granted.  Any distribution of this program or any part thereof
is strictly prohibited without prior written consent of M.I.T.

Title to copyright to this software and to any associated documentation
shall at all times remain with M.I.T. and LICENSEE agrees to preserve
same.  LICENSEE agrees not to make any copies except for LICENSEE'S
internal noncommercial use, or to use separately any portion of this
software without prior written consent of M.I.T.  LICENSEE agrees to
place the appropriate copyright notice on any such copies.

Nothing in this Agreement shall be construed as conferring rights to use
in advertising, publicity or otherwise any trademark or the name of
"Massachusetts Institute of Technology" or "M.I.T."

M.I.T. MAKES NO REPRESENTATIONS OR WARRANTIES, EXPRESS OR IMPLIED.  By
way of example, but not limitation, M.I.T. MAKES NO REPRESENTATIONS OR
WARRANTIES OF MERCHANTABILITY OR FITNESS FOR ANY PARTICULAR PURPOSE OR
THAT THE USE OF THE LICENSED SOFTWARE COMPONENTS OR DOCUMENTATION WILL
NOT INFRINGE ANY PATENTS, COPYRIGHTS, TRADEMARKS OR OTHER RIGHTS.
M.I.T. shall not be held liable for any liability nor for any direct,
indirect or consequential damages with respect to any claim by LICENSEE
or any third party on account of or arising from this Agreement or use
of this software.
*/

#include "mulGlobal.h"

#define SQRMAT 0		/* for rdMat(), wrMat() */
#define TRIMAT 1
#define COLMAT 2
#define LOWMAT 0		/* for blkMatSolve() */
#define UPPMAT 1
#define LO2TR 3			/* for matXfer() */
#define UP2TR 4
#define L11 0			/*1/4 matrix file names, see getName() */
#define U11 1
#define U12 2
#define L21 3
#define LTIL 4
#define UTIL 5
#define PMODE 0644		/* protection code used with write() */

#define SQDEX(i, j, siz) ((i)*(siz) + (j))
#define LODEX(i, j, siz) (((i)*((i)+1))/2 + (j))
#define UPDEX(i, j, siz) (((i)*(2*(siz) - (i) - 1))/2 + (j))

double *trimat, *sqrmat;	/* flattened triangular and square matrices */
int *real_index;		/* real_index[i] = index in full array
				   (w/dummies) corresponding to entry i in
				   the condensed array used by blk routines */

/*
  index into a square siz x siz matrix stored linearly
*/
int sqrdex(i, j, siz)
int i, j, siz;			/* row, column and size */
{
  return(i*siz + j);
}

/*
  index into a lower triangular (w/diagonal) siz x siz matrix stored linearly
*/
int lowdex(i, j, siz)
int i, j, siz;			/* row, column and size */
{
  int ret;
  if(j > i || (ret = i*(i+1)/2 + j) > siz*(1+siz)/2) {
    fprintf(stderr, 
	    "lowdex: bad indices for lower triangular i=%d  j=%d siz =%d\n",
	    i, j, siz);
    exit(1);
  }
  else return(ret);
}

/*
  index into an upper triangular (w/diagonal) siz x siz matrix stored linearly
*/
int uppdex(i, j, siz)
int i, j, siz;			/* row, column and size */
{
  int ret;
  
  if(j < i || (ret = i*siz - i*(i-1)/2 + j - i) > siz*(1+siz)/2) {
    fprintf(stderr, 
	    "uppdex: bad indices for upper triangular i=%d  j=%d siz =%d\n",
	    i, j, siz);
    exit(1);
  }
  else return(ret);
}

/*
  for debug only - dumps upper left corner of a matrix
*/
void dumpMatCor(mat, vec, fsize)
int fsize;			/* full matrix size for flat storage */
double **mat, *vec;		/* does first one that's not null */
{
  int i, j, size = 5;

  if(mat != NULL) {		/* full matrix */
    fprintf(stdout, "\nUPPER LEFT CORNER - full %dx%d matrix\n", fsize, fsize);
    for(i = 0; i < MIN(size, fsize); i++) {
      for(j = 0; j < MIN(size, fsize); j++) {
	fprintf(stdout, "%.5e ", mat[i][j]);
      }
      fprintf(stdout, "\n");
    }

  }
  else if(vec != NULL){			/* flat matrix as in blkDirect.c */
    fprintf(stdout, "\nUPPER LEFT CORNER - flat %dx%d matrix\n", fsize, fsize);
    for(i = 0; i < MIN(size, fsize); i++) {
      for(j = 0; j < MIN(size, fsize); j++) {
	fprintf(stdout, "%.5e ", vec[SQDEX(i, j, fsize)]);
      }
      fprintf(stdout, "\n");
    }
  }
}

/*
  gets the file name from the flag
*/
char *getName(file, name)
int file;
char *name;
{
  if(file == L11 || file == L21 || file == LTIL) name[0] = 'L';
  else name[0] = 'U';

  if(file == L11 || file == U11 || file == U12) name[1] = '1';
  else if(file == L21) name[1] = '2';
  else name[1] = 'T';

  if(file == L11 || file == U11 || file == L21) name[2] = '1';
  else if(file == U12) name[2] = '2';
  else name[2] = 'I';

  name[3] = '\0';
  return(name);
}

/*
  converts a square matrix to its transpose (used to write columnwise)
*/
void transpose(mat, siz)
double *mat;
int siz;
{
  int i, j, sqrdex();
  double temp;

  for(i = 0; i < siz; i++) {
    for(j = 0; j < i; j++) {
      temp = mat[SQDEX(i, j, siz)];
      mat[SQDEX(i, j, siz)] = mat[SQDEX(j, i, siz)];
      mat[SQDEX(j, i, siz)] = temp;
    }
  }
}


/*
  writes full or triangular matrices 
*/
void wrMat(mat, siz, file, type)
double *mat;
int siz, file, type;		/* siz is #rows and cols */
{
  int i, j, sqrdex(), ds = sizeof(double), fdis;
  int realsiz, actsiz;			/* size in chars */
  char name[3], *getName();			/* name of file */

  /* figure the real size */
  if(type == TRIMAT) realsiz = ds*siz*(siz+1)/2;
  else if(type == SQRMAT || type == COLMAT) realsiz = sizeof(double)*siz*siz;
  else {
    fprintf(stderr, "wrMat: bad type flag %d\n", type);
    exit(1);
  }

  /* figure name of file and create, open to write */
  if((fdis = creat(getName(file, name), PMODE)) == -1) {
    fprintf(stderr, "wrMat: can't creat '%s'\n", name);
    exit(1);
  }

  fprintf(stderr,"Writing %s...", name);

  /* write the data and close */
  if(type == COLMAT) transpose(mat, siz);	/* store columnwise */
  if((actsiz = write(fdis, (char *)mat, realsiz)) != realsiz) {
    fprintf(stderr, 
	    "wrMat: buffer write error to '%s,' wrote %d of %d dbls\n", 
	    name, actsiz/ds, realsiz/ds);
    exit(1);
  }
  close(fdis);

  fprintf(stderr, "done.\n");
}

/*
  reads full or triangular matrices 
*/
void rdMat(mat, siz, file, type)
double *mat;
int siz, file, type;		/* siz is #rows and cols */
{
  int i, j, sqrdex(), fdis;
  int realsiz;			/* size in chars */
  char name[3], *getName();   	/* name of file */

  /* figure the real size */
  if(type == TRIMAT) realsiz = sizeof(double)*siz*(siz+1)/2;
  else if(type == SQRMAT) realsiz = sizeof(double)*siz*siz;
  else {
    fprintf(stderr, "rdMat: bad type flag %d\n", type);
    exit(1);
  }

  /* figure name of file and open to read */
  if((fdis = open(getName(file, name), 0)) == -1) {
    fprintf(stderr, "rdMat: can't open '%s'\n", name);
    exit(1);
  }

  fprintf(stderr,"Reading %s...", name);

  /* read the data and close */
  if(realsiz != read(fdis, (char *)mat, realsiz)) {
    fprintf(stderr, "rdMat: read error to '%s'\n", name);
    exit(1);
  }
  close(fdis);

  fprintf(stderr, "done.\n");
}

/*
  transfer part of the square matrix to the triangular matrix
*/
void matXfer(matsq, matri, siz, type)
int siz, type;
double *matsq, *matri;
{
  int i, j, uppdex(), sqrdex(), temp;

  if(type == UP2TR) {		/* mv upper triangular part */
    for(i = 0; i < siz; i++) {	/* from row zero up */
      for(j = i; j < siz; j++) {
	temp = UPDEX(i, j, siz);
	matri[temp] = matsq[SQDEX(i, j, siz)];
      }
    }
  }
  else if(type == LO2TR) {	/* mv lower triangular part, put 1's on diag */
    for(i = 0; i < siz; i++) {	/* from row zero up */
      for(j = 0; j <= i; j++) {
	if(i == j) matri[LODEX(i, j, siz)] = 1.0;
	else matri[LODEX(i, j, siz)] = matsq[SQDEX(i, j, siz)];
      }
    }
  }
  else {
    fprintf(stderr, "matXfer: bad type %d\n", type);
    exit(1);
  }
}  

/*
  does many rhs, triangular problem to get L21 and U12 for blk factorization
*/
void blkMatsolve(matsq, matri, siz, type)
int siz, type;
double *matsq, *matri;
{
  int i, j, k, sqrdex(), lowdex(), uppdex();
  extern int fulldirops;

  if(type == LOWMAT) {
    /* a row in the result matrix is a linear comb of the previous rows */
    for(i = 1; i < siz; i++) {	/* loop on rows */
      fprintf(stderr, "%d ", i);
      for(j = 0; j < siz; j++) { /* loop on columns */
	for(k = 0; k < i; k++) {	/* loop on previous rows */
	  matsq[SQDEX(i, j, siz)]
	      -= (matri[LODEX(i, k, siz)]*matsq[SQDEX(k, j, siz)]);
	  fulldirops++;
	}
      }
    }
    fprintf(stderr, "\n");
  }
  else if(type == UPPMAT) {
    /* a column in the result matrix is a linear comb of the prev columns */
    for(j = 0; j < siz; j++) {	/* loop on columns */
      fprintf(stderr, "%d ", j);
      for(i = 0; i < siz; i++) { /* loop on rows */
	for(k = 0; k < j; k++) {	/* loop on previous columns */
	  matsq[SQDEX(i, j, siz)]
	      -= (matri[UPDEX(k, j, siz)]*matsq[SQDEX(i, k, siz)]);
	  fulldirops++;
	}
	matsq[SQDEX(i, j, siz)] /= matri[UPDEX(j, j, siz)];
	fulldirops++;
      }
    }
    fprintf(stderr, "\n");
  }
  else {
    fprintf(stderr, "blkMatsolve: bad type %d\n", type);
    exit(1);
  }

}

/*
  figures the difference A11-(L21)(U12) - part of block factorization
*/
void subInnerProd(matsq, matri, siz, matl, matu)
double *matsq, *matri;
int siz, matl, matu;		/* size in doubles; matrices to multiply */
{
  int i, j, k, matrisiz, rowlim, rowliml, colimu, fdl, fdu, lowdex(), sqrdex();
  int froml, fromu, ds = sizeof(double), readl, readu;
  char *getName(), name[3];
  double *matriu, temp;
  extern int fulldirops;
  extern double lutime;

  /* figure how many doubles will fit in matrix
  matrisiz = siz*(siz + 1)/2; */

  /* figure max number of square matrix rows that can fit (truncate) */
  rowlim = (siz+1)/2;	/* total #rows */
  rowliml = rowlim/2;		/* number of rows in first chunk */
  colimu = rowlim - rowliml;	/* number of rows in second chunk */
  matriu = matri + rowliml*siz; /* pointer to second block of rows */

  /* matri always holds rows from matl; matri1 holds columns from matu */
  /* open the relvant files */
  if((fdl = open(getName(matl, name), 0)) == -1) {
    fprintf(stderr, "subInnerProd: can't open '%s'\n", name);
    exit(1);
  }

  /* for each matl chunk, read in all chunks of matu and do inner products */
  for(froml = 0; froml < siz; froml += rowliml) { /* loop on rows in l part */
    readl = read(fdl, (char *)matri, rowliml*siz*ds);
    if(readl % (siz*ds) != 0) {	/* must read in row size chunks */
      fprintf(stderr, "subInnerProd: read error from '%s'\n", 
	      getName(matl, name));
      exit(1);
    }
    readl /= (siz*ds);
    if((fdu = open(getName(matu, name), 0)) == -1) { /* (re)open u part */
      fprintf(stderr, "subInnerProd: can't open '%s'\n", name);
      exit(1);
    }
    for(fromu = 0; fromu < siz; fromu += colimu) { /* loop on cols in u part */
      fprintf(stderr, "%d-%d ", froml, fromu);
      readu = read(fdu, (char *)matriu, colimu*siz*ds);
      if(readu % (siz*ds) != 0) {	/* must read in col size chunks */
	fprintf(stderr, "subInnerProd: read error from '%s'\n", 
		getName(matu, name));
	exit(1);
      }
      readu /= (siz*ds);
      /* do the inner product/subtractions possible with these chunks */
      starttimer;
      for(i = 0; i < readl; i++) { /* loop on rows */
	for(j = 0; j < readu; j++) {	/* loop on columns */
	  for(k = 0, temp = 0.0; k < siz; k++) { /* inner product loop */
	    /* matriu indices are flipped since it was stored as columns */
	    temp += matri[SQDEX(i, k, siz)] * matriu[SQDEX(j, k, siz)];
	    fulldirops++;
	  }
	  /* indices must be offset to get right square matrix entry */
	  matsq[SQDEX(i + froml, j + fromu, siz)] -= temp;
	  fulldirops++;
	}
      }
      stoptimer;
      lutime += dtime;
    }

    /* close the u file so it can be reread for next set of l rows */
    close(fdu);
  }
  fprintf(stderr, "\n");
  close(fdl);
}


/*