s_runsis1.cc

矩阵奇异分解(svd)最新c++版本

// File: s_runsis1.cc
// Author: Suvrit Sra
// Date: 23 Nov 03

/*************************************************************************
 THE ORIGINAL SVDPAKC COPYRIGHT
                           (c) Copyright 1993
                        University of Tennessee
                          All Rights Reserved                          
 *************************************************************************/


#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <cerrno>
#include <cstring>
#include <unistd.h>
#include <fcntl.h>

#include "s_runsis1.h"

using namespace ssvd;

s_runsis1::s_runsis1(char* optionsFile)
{
  options = new OptionsStruct();

  /* open files for input/output */
  if (!(options->in1 = fopen(optionsFile, "r"))) { 
    fprintf(stderr, "s_runsis1()::could not open file %s for reading\n", optionsFile);
    exit (-1);
  }
    
  fscanf(options->in1, "%s %s %s", options->out1, options->out2, options->out3);

  if (!(options->of1 = fopen(options->out1, "w"))) { 
    fprintf(stderr, "s_runsis1()::could not open output file %s \n", options->out1);
    exit (-1);
  }


  options->of1 = fopen(options->out1, "w");
  if (!options->of1) {
    fprintf(stderr, "s_runsis1()::Error opening %s. Quitting!\n", options->out1);
    exit(-1);
  }

  options->of2 = fopen(options->out2, "w");
  if (!options->of2) {
    fprintf(stderr, "s_runsis1()::Error opening %s. Quitting!\n", options->out2);
    exit(-1);
  }

  int asc;
  fscanf(options->in1,"%s%ld%ld%ld%lf%s%d", options->name, &options->maxsubsp,
	 &options->numextra, &options->km, &options->eps, options->vtf, &asc);

  if (!(strcmp(options->vtf, "TRUE"))) {
    options->vectors = true;
    if (asc > 0) {
      options->ascii = true;
      options->vecs = fopen(options->out3, "w");
      if (!options->vecs) {
	fprintf(stderr, "s_runsis1()::cannot open output file %s \n", options->out3);
	exit (-1);
      }
    } else {
      if ((options->fpo2 = open(options->out3, O_CREAT | O_RDWR)) == -1) {
	fprintf(stderr, "s_runsis1()::cannot open output file %s \n", options->out3);
	exit (-1);
      }
    }
  }
  else
    options->vectors = false;
}


s_runsis1::s_runsis1(OptionsStruct* o)
{ 
  if (!o)
    setupDefault();
  else 
    options = o;

  /* Now open the files etc. specified by the options*/
  if (!(options->of1 = fopen(options->out1, "w"))) { 
    fprintf(stderr, "could not open output file %s \n", options->out1);
    exit (-1);
  }
  options->of2 = fopen(options->out2, "w");
  if (!options->of2) {
    fprintf(stderr, "Error opening %s. Quitting!\n", options->out2);
    exit(-1);
  }

  if (options->vectors) {
    if (options->ascii) {
      options->vecs = fopen(options->out3, "w");
      if (!options->vecs) {
	fprintf(stderr, "cannot open output file %s \n", options->out3);
	exit (-1);
      }
    } else {
      if ((options->fpo2 = open(options->out3, O_CREAT | O_RDWR)) == -1) {
	fprintf(stderr, "cannot open output file %s \n", options->out3);
	exit (-1);
      }
    }
  }
  //if (init(options->prefix, options->txx) <  0) {
  //  fprintf(stderr, "Could not build sparse matrix %s\n", options->prefix);
  //  exit (-1);
  //}
}

void s_runsis1::setupDefault()
{
  options = new OptionsStruct();
  /* Set up certain default values of parameters. Might not generally work
     though */
  
}


int s_runsis1::runIt()
{

  float t0, exetime;
  long k,  i, ii, n, size1, size2, vectors;
  double *x[NSIG], *cx, *f, *d, *u, *tptr ;
  long em, p, em2, imem;
  double dsum, tmp1, tmp2, tmp3, tmp4, xnorm;
  
  /*write header of output file*/   

  fprintf(options->of2, "\n\n ----------------------------------------\n INTERMEDIATE OUTPUT PARMS:\n\n M:=CURRENT DEGREE OF CHEBYSHEV POLYNOMIAL\n S:=NEXT ITERATION STEP\n G:=NUMBER OF ACCEPTED EIGENVALUES OF B\n H:=NUMBER OF ACCEPTED EIGEN VECTORS OF B\n F:=VECTOR OF ERROR FOR EIGENPAIRS OF B\n ---------------------------------------------\n     M         S       G       H         F     \n ---------------------------------------------");

  n = nrow + ncol;
  


  if (n >= NMAX || m_nz > NZMAX) {
    fprintf(stderr,"runIt()::sorry, your matrix is too big\n");
    return -1;
  }

  em = options->maxsubsp;
  p = em + options->numextra;
  //  fprintf(stderr, "stuff  == %u %u %u\n", em, n, p);
  em2 = options->maxsubsp;
  /*******************************************************************
   * allocate memory						       *
   * pointr - column start array of harwell-boeing sparse matrix     *
   *          format                                       (ncol+1)  *
   * rowind - row indices array of harwell-boeing format   (nnzero)  *
   * value  - nonzero values array of harwell-boeing sparse matrix   *
   *          format                                       (nnzero)  *
   * x      - 2 dimensional array of  iteration vectors    (NSIG*n)  *
   * d      - array of eigenvalues of B                    (p)       *
  (singular values of A)                        *
  * f      - temporary storage array                      (p)       *
  * cx     - temporary storage array (control quantities) (n)       *
  * u    - work array                                     (n)       *
  *******************************************************************/
  size1 = sizeof(double) * (p* 2 + 2 * n + NSIG * n);
  //  size2 = sizeof(long) * (ncol + nnzero + 1);
  //if (!(pointr = (long *)   malloc(size2))   ||
  if (!(tptr  = (double *) malloc(size1))) {
    perror("MALLOC FAILED in MAIN()");
    fprintf(stderr,"runIt()::Tried to allocate %ul bytes\n", size1);
    return -2;
  }

  imem=size1; // - (nnzero * sizeof(double));     

  for (ii=0; ii < NSIG; ii++){
    x[ii]= tptr;
    tptr += n;
  }
  d=tptr;
  tptr+=p;
  f=tptr;
  tptr+=p;
  cx=tptr;
  tptr+=n;
  u=tptr;
	 

  sing = d;
  alpha = ZERO;

  /*compute alpha as 1-nrm of matrix A */
  for (k=0;k<ncol;k++){
    dsum=ZERO;
    for (i=pointr[k];i<=pointr[k+1]-1;i++)
      dsum+=value[i];
    alpha=dmax(alpha,dsum);
  }
  exetime = timer();
    
  /*  call ritzit */

  ritzit(options->of2, n,p,options->km,
		 options->eps,options->maxsubsp,x,d,f,cx,u,&imem);

  exetime = timer() - exetime;
  write_data(n,options->km,em2,options->maxsubsp,p,
			 imem, options->vectors,
			 options->eps, "NONE", options->name);
  t0=timer();

  if (options->vectors)
    write_header(options->fpo2,  options->vecs, options->ascii,nrow, ncol, em2, "sis1");
  
  for (j=0;j< em2; j++) {
    myopb(n, &x[j][0], cx, ZERO);
    tmp1=ddot(n,&(x[j][0]), 1, cx, 1);
    tmp2=ddot(nrow,&(x[j][0]), 1, &(x[j][0]), 1);
    tmp3=ddot(ncol, &(x[j][nrow]),1, &(x[j][nrow]), 1);
    tmp4=sqrt(tmp2 + tmp3);
    daxpy(n, -tmp1, &(x[j][0]), 1, cx, 1);
    xnorm= ddot(n, cx,1, cx,1);
    xnorm=sqrt(xnorm);
    f[j]=xnorm/tmp4;
    d[j]=fabs(tmp1);
    if (options->vectors) {
      if (options->ascii) {
	for (int zz = 0; zz < n; zz++) {
	  fprintf(options->vecs, "%f ", x[j][zz]);
	  if (zz == nrow-1 ||  zz == n-1) 
	    fprintf(options->vecs, "\n");
	}
      } else {
	write(options->fpo2, (char *)&x[j][0], n * sizeof(double));
      }
    }
  }
  exetime=exetime + timer() -t0;
  
  fprintf(stdout,"\n ...... SISVD EXECUTION TIME= %10.2e\n",exetime);
  fprintf(stdout," ......\n ......");
  fprintf(stdout," MXV   = %12.ld\n", mxvcount);
  fprintf(stdout," ......    COMPUTED SINGULAR VALUES");
  fprintf(stdout,"  (RESIDUAL NORMS)\n ......");

  for (ii=0;ii<em2;ii++)
    fprintf(stdout,"\n ...... %3.ld   %22.14e  (%11.2e)",
	    ii+1,     d[ii],    f[ii]);

  
  fprintf(options->of1, "%d\n", em2);
  for (ii=0; ii < em2; ii++)
    fprintf(options->of1, "%.14f\n", d[ii]);

  fprintf(stdout,"\n");
  if (options->in1)
    fclose(options->in1);
  if (options->vectors) close(options->fpo2);
  fclose(options->of2);
  fclose(options->of1);
  return 0;
}

void s_runsis1::write_data(long n, long km, long em2, long em, long p, 
						   long imem, long vectors, double eps, 
						   char *title, char *name)
{
  char svectors;
  if (vectors) svectors='T'; else svectors='F';
  fprintf( stdout,"\n ...\n ... SOLVE THE CYCLIC EIGENPROBLEM\n");
  fprintf(stdout," ... NO. OF EQUATIONS           = %10.ld\n",n);
  fprintf( stdout," ... MAX. NO. OF ITERATIONS     = %10.ld\n ",km);
  fprintf(stdout,"... NO. OF DESIRED EIGENPAIRS  = %10.ld\n ",em2);
  fprintf(stdout,"... NO. OF APPROX. EIGENPAIRS  = %10.ld\n ",em);
  fprintf(stdout,"... INITIAL SUBSPACE DIM.      = %10.ld\n ",p);
  fprintf(stdout,"... FINAL   SUBSPACE DIM.      = %10.ld\n ",p-em);
  fprintf(stdout,"... MEM REQUIRED (BYTES)       = %10.ld\n",imem);
  fprintf(stdout," ... WANT S-VECTORS? [T/F]      = %10.c\n ",svectors);
  fprintf(stdout,"... ALPHA                      = %10.2e\n",alpha);
  fprintf(stdout," ... TOLERANCE                  = %10.2e\n ",eps);
  fprintf(stdout,"... NO. OF ITERATIONS TAKEN    = %10.ld\n ",ksteps);
  fprintf(stdout,"... MAX. DEG. CHEBYSHEV POLY.  = %10.ld\n ...\n ", maxm);
  fprintf(stdout,"%s\n %s\n ... NO. OF TERMS     (ROWS)    = %10.d\n",title, name, nrow);
  fprintf(stdout," ... NO. OF DOCUMENTS (COLS)    = %10.ld\n ",ncol);
  fprintf(stdout,"... ORDER OF MATRIX B          = %10.ld\n ...\n", n);
  fflush(stdout);
}