s_runsis2.cc

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

// File: s_runsis2.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_runsis2.h"

using namespace ssvd;

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

  /* open files for input/output */
  if (!(options->in1 = fopen(optionsFile, "r"))) { 
    fprintf(stderr, "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, "could not open output file %s \n", options->out1);
    exit (-1);
  }


  options->of1 = fopen(options->out1, "w");
  if (!options->of1) {
    fprintf(stderr, "Error opening %s. Quitting!\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);
  }

  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, "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);
      }
    }
  }
  else
    options->vectors = false;
}


s_runsis2::s_runsis2(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->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_runsis2::setupDefault()
{
  options = new OptionsStruct();
  /* Set up certain default values of parameters. Might not generally work
     though */
  
}


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

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


    n = ncol; 
 
    if (m_nz > NZMAX) {
      fprintf(stderr,"sorry, your matrix is too big\n");
      return -1;
    }
    
    p=  options->maxsubsp + options->numextra;
    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)       *
                         (squares of the singular values of A)         *
     * f      - temporary storage array                      (p)       *
     * cx     - temporary storage array (control quantities) (n)       *
     * u      - work array                                   (nrow)    *
     *******************************************************************/
    size1 = sizeof(double) * (p* 2 + n  + nrow + NSIG * n);
    //size2 = sizeof(long) * (ncol + nnzero + 1);
    //if (!(pointr = (long *)   malloc(size2))   ||
    if ( !(tptr  = (double *) malloc(size1))) {
      perror("MALLOC FAILED in MAIN()");
      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;
    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, "sis2");
    }
    for (j=0;j<em2;j++) {
      mopb(n, &x[j][0], cx);
      tmp1=ddot(n,&(x[j][0]), 1, cx, 1);
      daxpy(n, -tmp1, &(x[j][0]), 1, cx, 1);
      tmp1=sqrt(tmp1);
      xnorm=sqrt(ddot(n, cx, 1, cx, 1));
      /* multiply by matrix A to get (scaled) left S-vector */
      opa(n, &x[j][0], u);
      tmp2=ONE/tmp1;
      dscal(nrow, tmp2, u, 1);
      xnorm=xnorm*tmp2;
      f[j]=xnorm;
      d[j]=tmp1;
      if (options->vectors) {
	if (options->ascii) {
	  for (int zz = 0; zz < nrow; zz++)
	    fprintf(options->vecs, "%f ", u[zz]);
	  fprintf(options->vecs, "\n");
	  for (int zz = 0; zz < n; zz++)
	    fprintf(options->vecs, "%f ", x[j][zz]);
	  fprintf(options->vecs, "\n");
	} else {
	  write(options->fpo2, (char *)u, nrow * sizeof(double));
	  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 ......    COMPUTED SINGULAR VALUES  (RESIDUAL NORMS)\n ......", mxvcount);
    
    for (ii=0;ii<em2;ii++)
      fprintf(stdout,"\n ...... %3.ld   %22.14e  (%11.2e)",
	      ii+1,     d[ii],    f[ii]);
    
    
    fprintf(stdout,"\n");
    
    fprintf(options->of1, "%d\n", em2);
    for (ii=0; ii < em2; ii++)
      fprintf(options->of1, "%.14f\n", d[ii]);
    if (options->in1)
      fclose(options->in1);
    if (options->vectors) close(options->fpo2);
    fclose(options->of2);
    fclose(options->of1);
    return 0;
}


/*********************************************************************/
void s_runsis2::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 [A^TA] 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," ... 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);
}