lsq.sav.txt

射线追踪程序

float cvar(float *x,int n)
{ int i;
  double M=0.;
  double V=0.;
  for (i=0; i<n; i++) M+=x[i];
  M/=(double)n;
  for (i=0; i<n; i++) V+=(x[i]-M)*(x[i]-M);
  V/=(double)n;
  return((float)V);
}

int main(int argc,char **argv)
{ struct matrix_t *G;
  float **R=NULL;
  float **Ginv=NULL;
  float *d=NULL,*d2=NULL,*sd=NULL,*res=NULL;
  float *m=NULL,*sm=NULL;
  float *dampfact=NULL;
  float damp0=-1.;
  int N,M,i,j,k,compinf=0,normdamp=0;
  int WS=0;
  FILE *f=NULL;
  char sysname[255],paraname[255],resolname[255],inmname[255],
       infoname[255],logname[255],varname[255],outname[255];
  char *basename=NULL;
  float trR,trI;
  float var,x;
  int VAR=0;

  for (i=1; i<argc; i++) {
      if (argv[i][0]=='-') {
         switch (argv[i][1]) {
                case 'v' : VAR=1;
                           break;
                case 's' : WS=1;
                           break;
                case 'i' : compinf=1;
                           break;
                case 'n' : if (argc<=i+1) goto error;
                           basename=argv[++i];
                           break;
                case 'D' : normdamp=1;
                case 'd' : if (argc<=i+1) goto usage;
                           if (!sscanf(argv[++i],"%f",&damp0)) goto error;
                           if (damp0<0.) goto error;
                           break;
                default  : goto error;
         }
      } else goto error;
  }
  if (!basename) goto usage;
  sprintf(sysname,"%s.sys",basename);
  sprintf(paraname,"%s.par",basename);
  sprintf(inmname,"%s.inm",basename);
  sprintf(resolname,"%s.res",basename);
  sprintf(infoname,"%s.inf",basename);
  sprintf(logname,"%s.log",basename);
  sprintf(varname,"%s.var",basename);
  sprintf(outname,"%s.out",basename);

  flog=fopen(logname,"wt");
  if (!flog) {
     perror(logname);
     exit(1);
  }

  if (damp0<0.) {
     fprintf(stderr,"Warning: damp=%f\n",damp0=0.);
     fprintf(flog,"Warning: damp=%f\n",damp0=0.);
  }

  fprintf(stderr,"Reading %s...",sysname);
  fprintf(flog,"Reading %s...",sysname);
  if (NULL==(f=fopen(sysname,"rb"))) {
     perror(sysname);
     fclose(flog);
     exit(1);
  }
  fprintf(stderr,"[G]... ");
  fprintf(flog,"[G]... ");
  G=readmatrix(f);
  if (!G) {
     fprintf(flog,"Error reading %s\n",sysname);
     fclose(flog);
     exit(1);
  }

  N=G->nl;
  M=G->nc;

  fprintf(stderr,"%d:%d...",N,M);
  fprintf(flog,"%d:%d...",N,M);
  fprintf(flog,"%d:%d...",N,M);
  d=(float*)calloc(N,sizeof(float));
  sd=(float*)calloc(N,sizeof(float));
  d2=(float*)calloc(N,sizeof(float));
  m=(float*)calloc(M,sizeof(float));
  dampfact=(float*)calloc(M,sizeof(float));
  res=(float*)calloc(M,sizeof(float));
  sm=(float*)calloc(M,sizeof(float));
  if ((!d)||(!m)||(!sd)||(!sm)) {
     fprintf(stderr,"malloc error (x)\n");
     fprintf(flog,"malloc error (x)\n");
     fclose(flog);
     exit(1);
  }
  fprintf(stderr,"[d]... ");
  fprintf(flog,"[d]... ");
  if (fread(d,sizeof(float),N,f)!=N) {
     fprintf(stderr,"Error (2) reading %s\n",sysname);
     fprintf(flog,"Error (2) reading %s\n",sysname);
     fclose(flog);
     exit(1);
  }
  fprintf(stderr,"[sd]... ");
  fprintf(flog,"[sd]... ");
  if (fread(sd,sizeof(float),N,f)!=N) {
     fprintf(stderr,"Error (3) reading %s\n",sysname);
     fprintf(flog,"Error (3) reading %s\n",sysname);
     fclose(flog);
     exit(1);
  }
  fprintf(stderr,"[damp]... ");
  fprintf(flog,"[damp]... ");
  if (fread(dampfact,sizeof(float),M,f)!=M) {
     fprintf(stderr,"Error (4) reading %s\n",sysname);
     fprintf(flog,"Error (4) reading %s\n",sysname);
     fclose(flog);
     exit(1);
  }
  fprintf(stderr,"OK\n");
  fprintf(flog,"OK\n");
  fclose(f);

  fprintf(stderr,"Init... inv... ");
  fprintf(flog,"Init... inv... ");
  Ginv=(float**)calloc(M,sizeof(float*));
  for (i=0; i<M; i++)
      if ((Ginv[i]=(float*)calloc(N,sizeof(float)))==NULL) {
         fprintf(stderr,"malloc error (3:%d)\n",i);
         fprintf(flog,"malloc error (3:%d)\n",i);
         fclose(flog);
         exit(1);
      }
  fprintf(stderr,"res... ");
  fprintf(flog,"res... ");
  R=(float**)calloc(M,sizeof(float*));
  for (i=0; i<M; i++)
      if ((R[i]=(float*)calloc(M,sizeof(float)))==NULL) {
         fprintf(stderr,"malloc error (4:%d)\n",i);
         fprintf(flog,"malloc error (4:%d)\n",i);
         fclose(flog);
         exit(1);
      }
  fprintf(stderr,"\n\nvar(d)=%10.3e",var=cvar(d,N));
  fprintf(flog,"\n\nvar(d)=%10.3e",var);

  if (normdamp) for (j=0; j<M; j++) dampfact[j]*=damp0;
  else for (j=0; j<M; j++) dampfact[j]=damp0;

  inverse(G,Ginv,R,sm,dampfact,N,M);

  fprintf(stderr,"\nSolving... ");
  fprintf(flog,"\nSolving... ");
  memset(m,0,M*sizeof(float));
  for (i=0; i<M; i++)
      for (j=0; j<N; j++)
          m[i]+=Ginv[i][j]*d[j];

  if (WS) {
     fprintf(stderr,"\nWriting inverse matrix (%s)... ",inmname);
     fprintf(flog,"\nWriting inverse matrix (%s)... ",inmname);
     if (NULL==(f=fopen(inmname,"wb"))) {
        perror(inmname);
        fclose(flog);
        exit(1);
     }
     for (i=0; i<M; i++)
         for (j=0; j<N; j++)
             fwrite(&(Ginv[i][j]),sizeof(float),1,f);
     for (i=0; i<M; i++) fwrite(&(sm[j]),sizeof(float),1,f);
  }

  fprintf(stderr,"\nWriting model (%s)... ",paraname);
  fprintf(flog,"\nWriting model (%s)... ",paraname);
  if (NULL==(f=fopen(paraname,"wt"))) {
     perror(paraname);
     fclose(flog);
     exit(1);
  }
  for (j=0; j<M; j++)
      fprintf(f,"%d %e %e %e\n",j,m[j],sd[i]*sqrt(sm[j]),R[j][j]);
  fclose(f);

  fprintf(stderr," resolution (%s)... ",resolname);
  fprintf(flog," resolution (%s)... ",resolname);

  dump(resolname,R,M,M);

  memset(d2,0,N*sizeof(float));
  for (i=0; i<N; i++)
      for (j=0; j<M; j++)
          if (sparsegetval(G,i,j,&x))
             d2[i]+=x*m[j];

  fprintf(stderr,"\n\nvar(d)=");
  fprintf(flog,"\n\nvar(d)=");
  for (i=0; i<N; i++) d2[i]-=d[i];

  if (NULL==(f=fopen(outname,"wt"))) {
     perror(outname);
     fclose(flog);
     exit(1);
  }
  for (j=0; j<N; j++)
      fprintf(f,"%d %e %e\n",j,d[j],d2[j]);
  fclose(f);

  fprintf(stderr,"%10.3e\n",var=cvar(d2,N));
  fprintf(flog,"%10.3e\n",var);

  fprintf(stderr,"var(m)=");
  fprintf(flog,"var(m)=");
  fprintf(stderr,"%10.3e\n",var=cvar(m,M));
  fprintf(flog,"%10.3e\n",var);

  trR=0.;
  for (i=0; i<M; i++) trR+=R[i][i];

  fprintf(stderr,"tr(R)= %10.3e\n",trR);
  fprintf(flog,"tr(R)=%.2e\n",trR);

  if (VAR) {
     FILE *fvar=fopen(varname,"at");
     float vi=cvar(d,N);
     float vd=cvar(d2,N);
     float vm=cvar(m,M);
     if (!fvar) {
        perror(varname);
        exit(0);
     }
     fprintf(fvar,"%f %f %f %f %f %d %d\n",vi,vd,vm,trR,damp0,N,M);
     fclose(fvar);
  }

  if (compinf) {
     float *v;
     fprintf(stderr,"\nInformation matrix...");
     fprintf(flog,"\nInformation matrix...");
     if (NULL==(f=fopen(infoname,"wb"))) {
        perror(paraname);
        fclose(flog);
        exit(1);
     }
     trI=0.;
     for (j=0; j<N; j++)
         for (i=j; i<N; i++) {
             float inf=0.,fi=(float)i,fj=(float)j;
             for (k=0; k<M; k++)
                 if (sparsegetval(G,i,k,&x))
                    inf+=x*Ginv[k][j];
             fwrite(&fi,sizeof(float),1,f);
             fwrite(&fj,sizeof(float),1,f);
             fwrite(&inf,sizeof(float),1,f);
             if (i==j) { trI+=inf; continue; }
             fwrite(&fj,sizeof(float),1,f);
             fwrite(&fi,sizeof(float),1,f);
             fwrite(&inf,sizeof(float),1,f);
         }
     fclose(f);
     fprintf(stderr,"\ntr(I)=%.2e\n",trI);
     fprintf(flog,"\ntr(I)=%.2e\n",trI);
  }
  fprintf(stderr,"\n-----------\nOutput:\n");
  fprintf(flog,"\n-----------\nOutput:\n");
  fprintf(stderr,"   param:\t%s\n",paraname);
  fprintf(flog,"   param:\t%s\n",paraname);
  fprintf(stderr,"residual:\t%s\n",outname);
  fprintf(flog,"residual:\t%s\n",outname);
  fprintf(stderr,"   resol:\t%s\n",resolname);
  fprintf(flog,"   resol:\t%s\n",resolname);
  if (VAR) fprintf(stderr,"     var:\t%s\n",varname);
  if (VAR) fprintf(flog,"     var:\t%s\n",varname);
  if (compinf) {
     fprintf(stderr,"   info:\t%s\n",infoname);
     fprintf(flog,"   info:\t%s\n",infoname);
  }

  return(0);
error:
  fprintf(stderr,"%s: bad argument: '%s'\n",argv[0],argv[i]);
usage:
  fprintf(stderr,"%s -n name [-d|D damp] [-r] [-i]\n",
          argv[0]);
  return(1);
}