s_runtms1.cc

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

// File: s_runtms1.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_runtms1.h"

using namespace ssvd;

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

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

  long    *workptr3, **workptr4, *tempptr3, **tempptr4,length,memsize;
  long    k, i, id, size1, size2,  mem;
  long    n,iwall, mxv[3], itime[4], *itcgt, *titer;


  n = nrow + ncol;

   
   /*******************************************************************
    * 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)  *
    *******************************************************************/
   //size1 = sizeof(double) * (nnzero);
   //size2 = sizeof(long) * (ncol + nnzero + 1);
   /* calculate memory allocated */
   mem = m_nz + ncol + 1;

   //rowind = pointr + (ncol + 1);

  /* allocate memory for arrays */

   long int s = options->maxsubsp;

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

   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  ( n x  s)                          *
    *    sig                    s                                     *
    *    res                    s                                     *
    *                                                                 *
    *    (memory allocated separately)                                *
    *                                                                 *
    *    iwork(iw   )           s  ( s x  2)                          *
    *    lwork                  s                                     *
    *    titer                  s                                     *
    *    itcgt                  s                                     *
    *                                                                 *
    *******************************************************************/
   if(!(workptr1 = (double *)malloc(memsize))) {
      perror("SECOND MALLOC FAILED IN TMS1()");
      return -2;
   }
  
   memsize = sizeof(double *) * (s + s + s + 3 + 4 + s);
   mem   += memsize;

   if(!(workptr2 = (double **)malloc(memsize))) {
     perror("THIRD MALLOC FAILED IN TMS1()");
     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    = n * s;
   yy        = tempptr1;
   tempptr1 += length;
   y         = tempptr2;
   tempptr2 += s;
   j=0;
   for(i=0; i<length; i+=n) y[j++] = &yy[i];
 
   sig       = tempptr1;
   tempptr1 += s;

   sing     = sig;		// HOLD SINGULAR VALUES
   res       = tempptr1;
   tempptr1 += s;
 
   /* Allocate memory for logical array lwork (long 0= false and 1=true)
      and integer array iwork                                    */
 
   memsize = sizeof(long) * ( s * (2 + 1) + s + s );
   mem   += memsize;
 
   if(!(workptr3 = (long *)malloc(memsize)) ) {
      perror("FOURTH MALLOC FAILED IN TMS1()");
      return -2;
   }
 
   memsize = sizeof(long *) * 2;
   mem   += memsize;
 
   if(!(workptr4 = (long **)malloc(memsize))) {
      perror("FIFTH MALLOC FAILED IN TMS1()");
      return -2;
   }
 
   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;
 
   titer     = tempptr3;
   tempptr3 += s;
 
   itcgt     = tempptr3;
   tempptr3 += s;
 
    /* estimate alpha (via matrix 1-norm)
    * used for upper bound on max singular value of matrix A */

   alpha = norm_1();
   exetime = timer();

   /* make a trace min. run; exit upon error */ 
   fprintf(stdout, "s_tms1::runIt() about to start computing svd\n");
   ierr = tsvd1(stdout, n, options->nums, s, options->accel,
		options->tol, red, sig, options->maxit, 
                &mem, itcgt, titer, time, res,  mxv, work1,
                work2, work3, work4, work5, y, iwork, lwork);
   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])); 

   /* write output */

   fprintf(stdout, " ... \n");
   fprintf(stdout, " ... TRACE MINIMIZATION ON THE  \n");
   fprintf(stdout, " ... EQUIVALENT CYCLIC   E_PROBLEM\n");
   fprintf(stdout, " ... NO. OF EQUATIONS            =%10ld\n", n);
   fprintf(stdout, " ... MAX. NO. OF ITERATIONS      =%10ld\n", options->maxit);
   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",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. TSVD1       =%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)       %ld           %ld \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, "tms1");
       for (i = 0;  i < options->nums;  i++)
          write(options->fpo2, (void *) &y[i][0], sizeof(double) * n);
   }

   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;
}