s_runlas2.cc

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

// File: s_runlas2.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_runlas2.h"

#include <algorithm>

using namespace ssvd;

s_runlas2::s_runlas2(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 %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_runlas2::s_runlas2(OptionsStruct* o)
{
  DBG ("s_runlas2::s_runlas2(OptionsStruct*)\n");
  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 las2
 * program. These are just sort of place holders for quick debugging, the
 * best way to use las2 is to actually make use of an options file or pass
 * a well formed OptionsStruct ...
 */
void s_runlas2::setupDefault()
{
  options = new OptionsStruct();
  /* Set up certain default values of parameters. Might not generally work
     though */
  strcpy(options->name, "Default");
  strcpy(options->out1, "las2.sval");
  strcpy(options->out2, "las2.vecs");
  
  options->lanmax = 10;
  options->endl    = -1e-30;	
  options->endr    = 1e-30;
  options->nums  = 1;
  options->kappa   = 1e-6;
  options->vectors = false;
}


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

  DBG("In s_rulas2::runIt()\n");

  n = ncol;
  nn   = ncol + nrow;

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

  ierr = 0;
  size1 = sizeof(double) * (6 * n + nrow  + (n * options->lanmax));
  //size2 = sizeof(long) * (ncol + 1 + m_nz);
  //memory = size1 + size2 + sizeof(double) * (5 * n + lanmax * 4 + 1);
  //long * rind = (long*)malloc(size2);

  tptr1 = (double*) malloc(size1);
  
  if (!tptr1) {
    fprintf(stderr, "Error allocating memory in s_las2::runIt()\n");
    return -1;
  }

  r      = tptr1;
  tptr1 += n;
  ritz   = tptr1;
  tptr1 += n;
  bnd    = tptr1;
  tptr1 += n;
  d      = tptr1;
  tptr1 += n;
  ztemp     = tptr1;
  tptr1 += nrow;
  a      = tptr1;

  sing   = new double[options->nums * options->nums]; 	// keep a copy of
							// singular
							// values...overallocating for vagaries of las2
  for (int zz = 0; zz < options->nums * options->nums; zz++)
    sing[zz] = 0;

  exetime = timer();

  /* 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.;

  /* make a lanczos run; see landr for meaning of parameters */
  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) != 0) {
    DBG("s_las2::landr() failed\n");
    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;

  /* memory allocated for xv1, xv2 and s in landr() */
  if (options->vectors) memory += sizeof(double) * (nn * (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 (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.  This is done only for singular values that are
   * within the relative accuracy (kappa) */
  if (options->vectors) {
    DBG("s_runlas2::options->vectors\n");
    size1 = sizeof(double) * nrow;
    t0 = timer();
    id = 0;

    DBG("s_runlas2::writing out: " << nsig << " vectors\n");
    //    std::cerr << "Num of sig vectors == " << nsig << std::endl;
    for (i = 0; i < nsig; i++) {
      /* multiply by matrix B first */
      mopb(n, &xv1[id], xv2);
      tmp0 = ddot(n, &xv1[id], 1, xv2, 1);
      daxpy(n, -tmp0, &xv1[id], 1, xv2, 1);
      tmp0 = sqrt(tmp0);
      xnorm = sqrt(ddot(n, xv2, 1, xv2, 1));
      ida = id + ncol;

      /* multiply by matrix A to get (scaled) left s-vector */
      mopa(&xv1[id], &xv1[ida]);
      tmp1 = 1.0 / tmp0;
      dscal(nrow, tmp1, &xv1[ida], 1);
      xnorm *= tmp1;
      bnd[i] = xnorm;
      d[i] = tmp0;

      /* write left s-vector to output file */
      if (options->ascii) {
	for (int zz = 0; zz < nrow; zz++)
	  fprintf(options->vecs, "%f ", xv1[ida+ zz]);
	fprintf(options->vecs, "\n");
      } 
      else
	write(options->fpo2, (char *)&xv1[ida], size1);
      id += nn;
    }
    exetime += (timer() - t0);
    count1=(mxvcount-nsig)/2 + nsig;
    count2=(mxvcount-nsig)/2;
    fprintf(stdout, " ...... \n");
    fprintf(stdout, " ...... NO. MULTIPLICATIONS BY A  =%10ld\n", count1);
    fprintf(stdout, " ...... NO. MULT. BY TRANSPOSE(A) =%10ld\n", count2);
    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 (i = 0; i < nsig; i++) {
      fprintf(stdout, " ...... %3ld   %22.14E  (%11.2E)\n", i + 1, d[i], bnd[i]);
      sing[i] = d[i];
    }
  }
  else {
    DBG("s_runlas2 spitting out singular values\n");
    for (i = j; i >= 0; i--)
      if (bnd[i] > options->kappa * fabs(ritz[i])) break;
    nsig = j - i;
    
    count1=(mxvcount-nsig)/2 + nsig;
    count2=(mxvcount-nsig)/2;
    fprintf(stdout," ...... \n");
    fprintf(stdout," ...... NO. MULTIPLICATIONS BY A  =%10ld\n", count1);
    fprintf(stdout," ...... NO. MULT. BY TRANSPOSE(A) =%10ld\n", count2);
    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 (i = 1 ; i <= nsig; i++) {
      sing[i-1] = sqrt(ritz[k]);
      fprintf(stdout," ...... %3ld   %22.14E  (%11.2E)\n", i, sing[i-1], 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; // Giving segfaults so have commented it out...
  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;
}