s_runlas1.cc

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

// File: s_runlas1.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_runlas1.h"

using namespace ssvd;

s_runlas1::s_runlas1(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;			// Ascii svecs desired?
  fscanf (options->in1,"%s %ld %ld %lf %lf %s %lf %d", options->name, &options->lanmax,
            &options->nums, &options->endl, &options->endr, options->vtf, &options->kappa,&asc);

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

  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_runlas1::s_runlas1(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);
  //}
}

/*
 * Sets up some default values for the parameters used by the las1
 * program. These are just sort of place holders for quick debugging, the
 * best way to use las1 is to actually make use of an options file or pass
 * a well formed OptionsStruct ...
 */
void s_runlas1::setupDefault()
{
  options = new OptionsStruct();
  /* Set up certain default values of parameters. Might not generally work
     though */
  strcpy(options->name, "Default");
  strcpy(options->out1, "las1.sval");
  strcpy(options->out2, "las1.vecs");
  
  options->lanmax = 10;
  options->endl    = -1e-30;
  options->endr    = 1e-30;
  options->nums  = 1;
  options->kappa   = 1e-6;
  options->vectors = false;
}


int s_runlas1::runIt()
{
  float t0, exetime;
  double tmp0, tmp1, tmp2, tmp3, xnorm;
  double *r, *ritz, *bnd, *d, *tptr1;
  long k, i, id, n, lanmax, maxprs, nnzero, size1, size2, memory;
  long *tptr3, hcount;
  long ida, count1, count2;
  long nn;

  n = nrow + ncol;

  /* write header of output file */
  write_data(stdout, options->lanmax, options->nums, options->endl,
	     options->endr, options->vectors, options->kappa,
	     "NONE", options->name, nrow, ncol, n);

  size1 = sizeof(double) * (n * (options->lanmax + 6));
  size2 = sizeof(long) * (ncol + m_nz + 1);

  if (!(tptr1    = (double*) malloc(size1))) {
    fprintf(stderr, "Malloc failed...so quitting\n");
    exit(1);
  }
  //tptr1 = value; // Not needed nemore..
  
  /* calculate memory allocated for problem (including work space in landr */
  memory = size1 + size2 + sizeof(double) * (5 * n + 4 * options->lanmax + 1);	
  // above value is wrong coz. i did some stuff....
  
  //rowind = pointr + (ncol + 1);
  //tptr1 += m_nz;
  r      = tptr1;
  tptr1 += n;
  ritz   = tptr1;
  tptr1 += n;
  bnd    = tptr1;
  tptr1 += n;
  d      = tptr1;
  tptr1 += n;
  a      = tptr1;

  sing   = new double[n];
  for (int zz = 0; zz < n; zz++)
    sing[zz] = 0.0;

  /* to get a random starting vector, the first n cells must be
   * initialize to zero */
  for (int i = 0; i < n; i++) r[i] = 0.;
  
  exetime = timer();
  
  /* make a lanczos run; exit upon error */
  if (landr(options->fpo2, options->vecs, options->ascii, n, options->lanmax,
	    options->nums, m_nz, options->endl, options->endr, options->vectors, 
	    options->kappa, ritz, bnd, r)) {
    
    delete[] value;
    delete[] pointr;
    
    if (options->in1)
      fclose(options->in1);
    fclose(options->of1);
    if (options->vectors) {
      close(options->fpo2);
      free (xv1);
      free (xv2);
    }
    exit(-1);
  }
  exetime = timer() - exetime;

  /* calculating singular vectors requires extra work space
   * (xv1, xv2, s) in landr() */
  if (options->vectors) memory += sizeof(double) * 
		    (n*(j+1) + n + (j+1)*(j+1));

  /* print error code if not zero */

  if (ierr)
    fprintf(stdout," ... RETURN FLAG = %9ld ...\n",ierr);

  /* ...Print ritz values and error bounds */
  fprintf(stdout,"\n");
  fprintf(stdout," ...... ALLOCATED MEMORY (BYTES)= %10.2E\n",(float)memory);
  fprintf(stdout," ...... LANSO EXECUTION TIME=%10.2E\n", exetime);
  fprintf(stdout," ...... \n");
  fprintf(stdout," ...... NUMBER OF LANCZOS STEPS = %3ld       NEIG = %3ld\n",j+1,neig);
  fprintf(stdout," ...... \n");
  fprintf(stdout," ......         COMPUTED RITZ VALUES  (ERROR BNDS)\n");
  fprintf(stdout," ...... \n");
  for (int i = 0; i <= j; i++)
    fprintf(stdout," ...... %3ld   %22.14E  (%11.2E)\n",i+1,ritz[i],bnd[i]);

  /* compute residual error when singular values and vectors are	
   * computed */
  if (options->vectors) {
    t0 = timer();
    id = 0;

    for (int i = 0; i < nsig; i++) {

      /* multiply by 2-cyclic indefinite matrix */
      opb(n, &xv1[id], xv2);
      tmp0 = ddot(n, &xv1[id], 1, xv2, 1);
      tmp1 = ddot(nrow, &xv1[id], 1, &xv1[id], 1);
      tmp2 = ddot(ncol, &xv1[id + nrow], 1, &xv1[id + nrow], 1);
      tmp3 = sqrt(tmp1 + tmp2);
      daxpy(n, -tmp0, &xv1[id], 1, xv2, 1);
      xnorm = ddot(n, xv2, 1, xv2, 1);
      xnorm = sqrt(xnorm);
      bnd[i] = xnorm / tmp3;
      d[i] = fabs(tmp0);
      id += n;
    }  
    exetime += (timer() - t0);
    hcount=mxvcount/2;
    fprintf(stdout," ...... \n");
    fprintf(stdout," ...... NO. MULTIPLICATIONS BY A  =%10ld\n", hcount);
    fprintf(stdout," ...... NO. MULT. BY TRANSPOSE(A) =%10ld\n", hcount);
    fprintf(stdout,"\n");
    fprintf(stdout," ...... LASVD EXECUTION TIME=%10.2E\n",exetime);
    fprintf(stdout," ...... \n");
    fprintf(stdout," ......         NSIG = %4ld\n", nsig);
    fprintf(stdout," ...... \n");
    fprintf(stdout," ......         COMPUTED S-VALUES     (RES. NORMS)\n");
    fprintf(stdout," ...... \n");
    for (int i = 0; i < nsig; i++) {
      sing[i] = d[i];
      fprintf(stdout," ...... %3ld   %22.14E  (%11.2E)\n", i + 1, d[i], bnd[i]);
    }
  } else {
    int zet;
    for (zet = j; zet >= 0; zet--)
      if (bnd[zet] > options->kappa * fabs(ritz[zet])) break;
    nsig = j - zet;
    hcount=mxvcount/2;
    fprintf(stdout," ...... \n");
    fprintf(stdout," ...... NO. MULTIPLICATIONS BY A  =%10ld\n", hcount);
    fprintf(stdout," ...... NO. MULT. BY TRANSPOSE(A) =%10ld\n", hcount);
    fprintf(stdout,"\n");
    fprintf(stdout," ...... LASVD EXECUTION TIME=%10.2E\n", exetime);
    fprintf(stdout," ...... \n");
    fprintf(stdout," ......         NSIG = %4ld\n", nsig);
    fprintf(stdout," ...... \n");
    fprintf(stdout," ......         COMPUTED S-VALUES   (ERROR BNDS)\n");
    fprintf(stdout," ...... \n");

    k = j + 1 - nsig;
    for (int i = 1 ; i <= nsig; i++)  {
      sing[i-1] = ritz[k];
      fprintf(stdout," ...... %3ld   %22.14E  (%11.2E)\n", i, ritz[k], bnd[k]);
      k++;
    }
  }
  
  // Write out the svals to file
  fprintf(options->of1, "%d\n", nsig);
  // Sort sing in place ...
  //std::sort(sing, sing + nsig);
  
  // Now make it descending
  //double* tmp = new double[nsig];

  //for (int zz = 0; zz < nsig; zz++) {
  // tmp[zz] = sing[nsig-zz-1];
  //}

  //delete[] sing; // damn a memory leak rite here!!
  //sing = tmp;

  for (int zz = 0; zz < nsig; zz++) 
    fprintf(options->of1, "%.14f\n", sing[zz]);

  if (options->in1)
    fclose(options->in1);
  fclose(options->of1);

  return 0;
}