s_runtms2.cc

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

// File: s_runtms2.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_runtms2.h"

using namespace ssvd;

s_runtms2::s_runtms2(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", options->out1, options->out2);

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

  fscanf (options->in1,"%s %ld %ld %ld %lf %lf %s %ld %d", options->name,
	  &options->nums, &options->maxsubsp, &options->accel, &options->tol,
	  &options->red, options->vtf, &options->maxit, &asc);

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



s_runtms2::s_runtms2(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);
  }

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

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

int s_runtms2::runIt()
{
  double time[5],exetime,*sig,*res;
  double  *tempptr1,*workptr1;
  double  **tempptr2,**workptr2;

  long  *workptr3,**workptr4,*tempptr3,**tempptr4,memsize,length;

  long n,i,j, size1, size2, s,mem;
  long iwall,mxv[3], itime[4],*itcgt,*titer,maxd;

  n = imin (nrow, ncol);

  mem = m_nz + ncol + 1;

  /* allocate memory */
  s = options->maxsubsp;

  size1 = sizeof(double) * s;
  size2 = sizeof(long) * s;
  mem = 2 * (size1 + size2);

  if(!(sig = (double *) malloc(size1)) || 
     !(res = (double *) malloc(size1)) ||
     !(itcgt = (long *) malloc(size2)) || 
     !(titer = (long *) malloc(size2))) {
    perror("SECOND MALLOC FAILED IN TMS2()");
    return -2;
  }
  /* allocate memory for arrays */
  memsize = sizeof(double) *
    ((n * (s + s + 3)) + ( s*(4 + s) ) + s*(nrow+ncol) );

  mem += memsize;

  /************************************************************************
   *        Allocate work area and initialize pointers                     *
   *                                                                       *
   *        pointer                size                                    *
   *                                                                       *
   *        w1   (work1)           n  ( n x  s)                            *
   *        w2   (work2)           n  ( n x  s)                            *
   *        w3   (work3)           s  ( s x  s)                            *
   *        w4   (work4)           n  ( n x  3)                            *
   *        w5   (work5)           s  ( s x  4)                            *
   *        yy   (y     )          n  ( nrow+ncol x  s)
   *
   *                                                                       *
   *        (memory allocated separately)                                  *
   *                                                                       *
   *        iwork(iw   )           s  ( s x  2)                            *
   *        lwork                  s                                       *
   *                                                                       *
   ***********************************************************************/

  if(!(workptr1 = (double *)malloc(memsize))) {
    perror("THIRD MALLOC FAILED IN TMS2()");
    return -2;
  }

  memsize = sizeof(double *) * (s + s + s + 3 + 4 + s);
  mem += memsize;

  if(!(workptr2 = (double **)malloc(memsize))) {
    perror("FOURTH MALLOC FAILED IN TMS2()");
    return -2;
  }
  tempptr1 = workptr1;
  tempptr2 = workptr2;

  length    = n * s;
  w1        = tempptr1;
  tempptr1 += length;
  work1     = tempptr2;
  tempptr2 += s;
  j=0;
  for(i=0; i<length; i+=n) work1[j++] = &w1[i];

  length    = n * s;
  w2        = tempptr1;
  tempptr1 += length;
  work2     = tempptr2;
  tempptr2 += s;
  j=0;
  for(i=0; i<length; i+=n) work2[j++] = &w2[i];

  length    = s * s;
  w3        = tempptr1;

  tempptr1 += length;
  work3     = tempptr2;
  tempptr2 += s;
  j=0;
  for(i=0; i<length; i+=s) work3[j++] = &w3[i];

  length    = n * 3;
  w4        = tempptr1;
  tempptr1 += length;
  work4     = tempptr2;
  tempptr2 += 3;
  j=0;
  for(i=0; i<length; i+=n) work4[j++] = &w4[i];

  length    = s * 4;
  w5        = tempptr1;
  tempptr1 += length;
  work5     = tempptr2;
  tempptr2 += 4;
  j=0;
  for(i=0; i<length; i+=s) work5[j++] = &w5[i];

  length    = (nrow+ncol) * s;
  yy        = tempptr1;
  tempptr1 += length;
  y         = tempptr2;
  tempptr2 += s;
  j=0;
  for(i=0; i<length; i+=(nrow+ncol)) y[j++] = &yy[i];


  /* Allocate memory for logical array lwork ( long 0= false and 1=true) *
     and integer array iwork                                            */

  memsize = sizeof(long) * ( s * (2 + 1));
  mem   += memsize;

  if(!(workptr3 = (long *)malloc(memsize)) ) {
    perror("FIFTH MALLOC FAILED IN TMS2()");
    return -1;
  }

  memsize = sizeof(long *) * 2;
  mem   += memsize;

  if(!(workptr4 = (long **)malloc(memsize))) {
    perror("SIXTH MALLOC FAILED IN TMS2()");
    return -1;
  }

  tempptr3 = workptr3;
  tempptr4 = workptr4;

  length   = s * 2;
  iw       = tempptr3;
  tempptr3+= length;
  iwork    = tempptr4;
  tempptr4 += 2;
  j=0;

  for(i=0; i<length; i+=s) iwork[j++] = &iw[i];

  lwork   = tempptr3;
  tempptr3 += s;

  sing   = sig;

  /* estimate alpha (via matrix 1-norm)
   * used for uppper bound on max singular value of matrix A */
  alpha = norm_1();
  exetime = timer();

  /* make a trace min. run; exit upon error */ 
  ierr = tsvd2(stdout, n, options->nums, options->maxsubsp,options->accel, 
	       options->tol, options->red, sig, options->maxit, &mem, itcgt,
	       titer, time, res, mxv, work1, work2, work3, work4, work5,
	       y, iwork, lwork,&maxd);
  if(ierr) {
    if (options->in1)
      fclose(options->in1);
    fclose(options->of1);
    if(options->vectors) close(options->fpo2);
    free(workptr1);
    free(workptr2);
    free(workptr3);
    free(workptr4);
    return -1;
  } 
  exetime  = timer() - exetime;
  iwall    = (int) (fmod(exetime+1000000.0,1000000.0));
  itime[0] = (int) (100*(time[0]/time[4]));
  itime[1] = (int) (100*(time[1]/time[4]));
  itime[2] = (int) (100*(time[2]/time[4]));
  itime[3] = (int) (100*(time[3]/time[4])); 

  if(options->accel != 2) maxd = 0;

  /* write output */

  fprintf(stdout, " ... \n");
  fprintf(stdout, " ... TRACE MINIMIZATION ON THE  \n");
  fprintf(stdout, " ... EQUIVALENT A^TA  E_PROBLEM\n");
  fprintf(stdout, " ... NO. OF EQUATIONS            =%10ld\n", n);
  fprintf(stdout, " ... MAX. NO. OF ITERATIONS      =%10ld\n", options->maxit);
  fprintf(stdout, " ... MAX. DEG. CHEBYSHEV POLY.   =%10ld\n",maxd);
  fprintf(stdout, " ... NO. OF ITERATIONS TAKEN     =%10ld\n",titer[options->nums-1]);
  fprintf(stdout, " ... NO. OF DESIRED EIGENPAIRS   =%10ld\n",options->nums);
  fprintf(stdout, " ... INITIAL SUBSPACE DIM.       =%10ld\n",s);
  fprintf(stdout, " ... FINAL   SUBSPACE DIM.       =%10ld\n",s-options->nums);
  fprintf(stdout, " ... MEMORY REQUIRED (BYTES)     =%10ld\n",mem);
  fprintf(stdout, " ... JOB: 0=NO SHIFT, 1=SHIFT    =%10ld\n",options->accel);
  if (options->vectors)
    fprintf(stdout, " ... WANT S-VECTORS?   [T/F]   =           T\n");
  else
    fprintf(stdout, " ... WANT S-VECTORS?   [T/F]   =           F\n");
  fprintf(stdout, " ... ALPHA                       =%10.2E\n",alpha);
  fprintf(stdout, " ... FINAL RESIDUAL TOLERANCE    =%10.2E\n" ,options->tol);
  fprintf(stdout, " ... RESIDUAL REDUCTION TOL.     =%10.2E\n",options->red);
  fprintf(stdout, " ... MULTIPLICATIONS BY A        =%10ld\n",mxv[0]);
  fprintf(stdout, " ... MULTIPLICATIONS BY A^T      =%10ld\n",mxv[1]);
  fprintf(stdout, " ... TOTAL NUMBER OF MULT.       =%10ld\n",mxv[2]);
  fprintf(stdout, " ... IERR FROM SUBR. TSVD2       =%10ld\n\n",ierr);

  fprintf(stdout,"... CPU TIME BREAKDOWN:\n");
  fprintf(stdout,"... GRAM-SCHMIDT ORTHOG.          =%10.2E (%3ld%%)\n",
	  time[0],itime[0]);
  fprintf(stdout,"... SPECTRAL DECOMPOSITION        =%10.2E (%3ld%%)\n",
	  time[1],itime[1]);
  fprintf(stdout,"... CONVERGENCE CRITERIA          =%10.2E (%3ld%%)\n",
	  time[2],itime[2]);
  fprintf(stdout,"... CONJUGATE GRADIENT            =%10.2E (%3ld%%)\n",
	  time[3], itime[3]);
  fprintf(stdout,"... TOTAL CPU  TIME (SEC)         =%10.2E\n", time[4]);
  fprintf(stdout,"... WALL-CLOCK TIME (SEC)         =%10ld\n\n",iwall);
  //fprintf(stdout, " %s\n", title);
  fprintf(stdout, "           %s\n", options->name);
  fprintf(stdout, " ... NO. OF TERMS     (ROWS)   = %10ld\n", nrow);
  fprintf(stdout, " ... NO. OF DOCUMENTS (COLS)   = %10ld\n", ncol);
  fprintf(stdout, " ... ORDER OF MATRIX A         = %10ld\n", n);

  fprintf(stdout,"......\n");
  fprintf(stdout,"...... COMPUTED SINGULAR VALUES  (RESIDUAL NORMS)  T-MIN STEPS CG STEPS\n"); 
  fprintf(stdout,"......\n");
  fprintf(options->of1, "%d\n", options->nums);
  for(i=0; i< options->nums; i++) {
    fprintf(stdout, "... %2ld  %24.14E ( %10.2E)      %3ld           %3ld \n", i+1, sig[i], res[i],titer[i],itcgt[i]); 
    fprintf(options->of1, "%.14f\n", sig[i]);
  }
  if (options->vectors) {
    write_header(options->fpo2, options->vecs, options->ascii, nrow, ncol, options->nums, "tms2");
    for (i = 0;  i < options->nums;  i++)
      write(options->fpo2, (void *) &y[i][0], sizeof(double) * (nrow+ncol)); }

  if (options->in1)
    fclose(options->in1);
  fclose(options->of1);
  if (options->vectors) close(options->fpo2);
  free(workptr1);
  free(workptr2);
  free(workptr3);
  free(workptr4);
  return 0;
}