s_runbls2.cc

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

// File: s_runbls2.cc
// Author: Suvrit Sra
// Date: 23 Nov, 2003

/*************************************************************************
 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_runbls2.h"


using namespace ssvd;
s_runbls2::s_runbls2(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%ld%lf%s%d", options->name, &options->maxit, 
	  &options->nc, &options->nb, &options->nums, &options->tol, options->vtf,&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_runbls2::s_runbls2(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 bls2
 * program. These are just sort of place holders for quick debugging, the
 * best way to use bls2 is to actually make use of an options file or pass
 * a well formed OptionsStruct ...
 */
void s_runbls2::setupDefault()
{
  options = new OptionsStruct();
  /* Set up certain default values of parameters. Might not generally work
     though */
  strcpy(options->name, "Default");
  strcpy(options->out1, "bls2.sval");
  strcpy(options->out2, "bls2.vecs");
  
  options->maxit = 10;
  options->nc    = 100;		// upper bound krylov
  options->nb    = 4;
  options->nums  = 1;
  options->tol   = 1e-6;
  options->vectors = false;
}

int s_runbls2::runIt()
{
  double t0, exetime;
  long i, k;
  long ncold, nbold, nk, size1, size2, indexu, indexv;
  long memory, vectors;
  double tol,  *v, *u, *res, *tmpv, *tptr1;
  double tmp1, tmp2, xnorm;
  nn = imin (nrow, ncol);

  /* even though the validity of the parameters will be checked in the
   * SVD code, some parameter checking should also be done before
   * allocating memory to ensure that they are nonnegative */
  if (validate(stderr, m_nz, nrow, ncol, options->nc, options->nb, options->nums, options->tol)) {
    if (options->in1)
      fclose(options->in1);
    fclose(options->of1);
    if (options->vectors) close(options->fpo2);
    return -1;
  }

  size1 = sizeof(double) * ( options->nums * (ncol+2) + nrow * options->nc);
  if (!(tptr1  = (double *) malloc(size1))   ||
      !(ztempp = (double **) malloc(sizeof(double *) * options->nc)))    {
    perror("MALLOC FAILED in MAIN()");
    return -1;
  }

  /* calculate memory allocated for problem */
  memory = size1 + ncol + m_nz + 1 + sizeof(double *) * options->nc;

  v      = tptr1;
  tptr1 += ncol * options->nums;
  sing   = tptr1;
  tptr1 += options->nums;
  res    = tptr1;
  tptr1 += options->nums;
  ztemp  = tptr1;
  j = 0;
  k = options->nc * nrow;
  for (i = 0; i < k; i += nrow) ztempp[j++] = &ztemp[i];

  ncold = options->nc;
  nbold = options->nb;
  nk    = options->nums;
  mxvcount = 0;
  mtxvcount = 0;

  exetime = timer();
  if (blklan2(stderr, m_nz, nrow, ncol, options->nums, v, sing, ncold, 
	      nbold, options->tol, res, options->maxit, &nk, &options->nc, &options->nb, &memory)) {
    
    delete[] value;
    delete[] rowind;
    if (options->in1)
      fclose(options->in1);
    fclose(options->of1);
    if (options->vectors) close(options->fpo2);
    return -1;
  }

  exetime = timer() - exetime;
  if (options->vectors) {
    /*******************************************************************
     * allocate memory						       *
     * tmpv   -                                               (ncol)   *
     * u      -                                          (nrow * nk)   *
     *******************************************************************/
    size1 = sizeof(double) * (ncol + nk * nrow);
    memory += size1;
    if (!(tmpv = (double *) malloc(size1))) {
      perror("SECOND MALLOC FAILED in MAIN()");
      return -2;
    }
    u = tmpv + ncol;
    size1 = sizeof(double) * nrow;
    size2 = sizeof(double) * ncol;
    indexu = 0;
    indexv = 0;

    t0 = timer();
    if (options->vectors)
      write_header(options->fpo2, options->vecs, options->ascii, nrow, ncol, nk, "bls2");
    for (i = 0; i < nk; i++) {

      /* multiply by A'A */
      opb(nrow, ncol, &v[indexv], tmpv);
      tmp1 = ddot(ncol, &v[indexv], 1, tmpv, 1);
      daxpy(ncol, -tmp1, &v[indexv], 1, tmpv, 1);
      tmp1 = sqrt(tmp1);
      xnorm = sqrt(ddot(ncol, tmpv, 1, tmpv, 1));

      /* multiply by A to get scaled left singular vectors */
      opa(nrow, ncol, &v[indexv], &u[indexu]);
      tmp2 = ONE / tmp1;
      dscal(nrow, tmp2, &u[indexu], 1);
      res[i] = xnorm * tmp2;
      sing[i] = tmp1;

      if (options->ascii) {
	for (int zz = 0; zz < nrow; zz++)
	  fprintf(options->vecs, "%f ", u[indexu+zz]);
	fprintf(options->vecs, "\n");
	for (int zz = 0; zz < ncol; zz++)
	  fprintf(options->vecs, "%f ", v[indexv + zz]);
	fprintf(options->vecs, "\n");
	
      } else {
	write (options->fpo2, (char *)&u[indexu], size1);
	write (options->fpo2, (char *)&v[indexv], size2);
      }
      indexu += nrow;
      indexv += ncol;
    }
    exetime += timer() - t0;
  }
  /* write results to output file */
  fprintf(stdout,"\n");
  fprintf(stdout, " ... \n");
  fprintf(stdout, " ... HYBRID BLOCK LANCZOS [A^TA]\n");
  fprintf(stdout, " ... NO. OF EQUATIONS          =%10ld\n", nn);
  fprintf(stdout, " ... MAX. NO. OF ITERATIONS    =%10ld\n", options->maxit);
  fprintf(stdout, " ... NO. OF ITERATIONS TAKEN   =%10ld\n", iter);
  fprintf(stdout, " ... NO. OF TRIPLETS SOUGHT    =%10ld\n", options->nums);
  fprintf(stdout, " ... NO. OF TRIPLETS FOUND     =%10ld\n", nk);
  fprintf(stdout, " ... INITIAL BLOCKSIZE         =%10ld\n", nbold);
  fprintf(stdout, " ... FINAL   BLOCKSIZE         =%10ld\n", options->nb);
  fprintf(stdout, " ... MAXIMUM SUBSPACE BOUND    =%10ld\n", ncold);
  fprintf(stdout, " ... FINAL   SUBSPACE BOUND    =%10ld\n", options->nc);
  fprintf(stdout, " ... NO. MULTIPLICATIONS BY A  =%10ld\n", mxvcount);
  fprintf(stdout, " ... NO. MULT. BY TRANSPOSE(A) =%10ld\n", mtxvcount);
  fprintf(stdout, " ... TOTAL SPARSE MAT-VEC MULT.=%10ld\n", 
	  mxvcount + mtxvcount);
  fprintf(stdout, " ... MEMORY NEEDED IN BYTES    =%10ld\n", memory);
  if (options->vectors) 
    fprintf(stdout, " ... WANT S-VECTORS ON OUTPUT?           T\n");
  else
    fprintf(stdout, " ... WANT S-VECTORS ON OUTPUT?           F\n");
  fprintf(stdout, " ... TOLERANCE                 =%10.2E\n\n", options->tol);
  //fprintf(stdout, " %s\n", title);
  fprintf(stdout, "           %s\n", options->name);
  fprintf(stdout, " ... NO. OF TERMS (ROWS OF A)   = %8ld\n", nrow);
  fprintf(stdout, " ... NO. OF DOCS  (COLS OF A)   = %8ld\n", ncol);
  fprintf(stdout, " ... \n\n");
  fprintf(stdout, " ...... BLSVD EXECUTION TIME=%10.2E\n",exetime);
  fprintf(stdout, " ...... \n");
  fprintf(stdout, " ......         COMPUTED S-VALUES     (RES. NORMS)\n");
  fprintf(stdout, " ...... \n");
  fprintf(options->of1, "%d\n", nk);
  for (i = 0; i < nk; i++) {
    fprintf(stdout," ...... %3ld   %22.14E  (%11.2E)\n",i+1,sing[i],res[i]);
    fprintf(options->of1, "%.14f\n", sing[i]);
  }

  delete[] value;
  delete[] rowind;
  fclose(options->in1);

  fclose(options->of1);
  if (options->vectors) {
    free(tmpv);
    close(options->fpo2);
  }
  return 0;
}