svdpack.cc

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

// File: svdpack.cc -- Implements class svdpack
// Author: Suvrit Sra
// Date: 14th nov 2003


// (c) Suvrit Sra 2004 All Rights Reserved
/*************************************************************************
 THE ORIGINAL SVDPAKC COPYRIGHT
                           (c) Copyright 1993
                        University of Tennessee
                          All Rights Reserved                          
 *************************************************************************/

/* ****************************************************************
 * NOTES: Have to remove the LAPACK BLAS functions from here with calls
 * to just the system lapack library, because that might be ATLAS based
 * and the one here is just ordinary untuned blas
* ***************************************************************** */ 

#include <unistd.h>
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <iostream>
#include <fstream>


#include "svdpack.h"

using namespace ssvd;


void svdpack::write_header(long rfp, FILE* fp, bool ascii, size_t m, size_t n, int k, char* algo)
{
  if (ascii) {
    fprintf(fp, "%u %u %d %s\n", m, n, k, algo);
  } else {
    write(rfp, (char *)&m, sizeof(m));
    write(rfp, (char *)&n, sizeof(n));
    write(rfp, (char *)&k, sizeof(k));
    write(rfp, algo, 5);
    //    write(rfp, (char*)&nsig, sizeof(nsig));
  }
}

/*
 * This method loads the input sparse matrix...
 */
int svdpack::init(char* fname, char* txx)
{
  std::cerr << "Building sparsematrix from " << fname << " (" << txx << ") scaling\n";
  std::string scaling(txx);

  // Set up the file names that we are gonna open
  std::string dim(fname);
  std::string rows_file(dim + "_row_ccs");
  std::string cols_file(dim + "_col_ccs");
  std::string nz_file (dim + "_"+ scaling + "_nz");
  
  dim +=  "_dim";
  
  std::ifstream infile;
  infile.open(dim.data());
  if (infile.fail()) {
    std::cerr << "Error: could not open " << dim << std::endl;
    return -1;
  }
  
  infile >> nrow >> ncol >> m_nz;
  infile.close();
  infile.clear();
  
  //cout << "(m, n, nz)" << m_rows << " " << m_cols << " " << m_nz << endl;
  pointr  = new long[ncol+1];
  rowind  = new long[m_nz];
  value   = new double[m_nz];
  pointr[ncol] = m_nz; // Add the convenience value
  
  // Now read in all the columnpointers
  infile.open(cols_file.data());
  if (infile.fail()) {
    std::cerr << "Error: could not open " << cols_file << std::endl;
    return -1;
  }
  // Make sure the UB on i is correct.
  for (int i = 0; i < ncol+1; i++) 
    infile >> pointr[i];
  infile.close();
  infile.clear();
  
  infile.open(rows_file.data());
  if (infile.fail()) {
    std::cerr << "Error: could not open " << rows_file << std::endl;
    return -1;
  }

  for (int i = 0; i < m_nz; i++)
    infile >> rowind[i];
  infile.close();
  infile.clear();
  
    // Now read in the actual m_values
  infile.open(nz_file.data());
  if (infile.fail()) {
    std::cerr << "Error: could not open " << nz_file << std::endl;
    return -1;
  }
  for (int i = 0; i < m_nz; i++)
    infile >> value[i];
  infile.close();

//   for (int i = 0; i < ncol+1; i++) 
// 	std::cout << pointr[i] << " ";
//   std::cout << std::endl;

//   for (int i = 0; i < m_nz; i++) 
// 	std::cout << rowind[i] << " ";
//   std::cout << "\n And the nonzeros are\n";
//   for (int i = 0; i < m_nz; i++) 
// 	std::cout << value[i] << " ";

  return 0;
}

void svdpack::extern_init(SparseMatrix* s)
{
  pointr = s->getPointr();
  rowind = s->getIndx();
  value  = s->getVal();

  nrow   = s->numRows();
  ncol   = s->numCols();
  m_nz   = s->numNz();
}

void svdpack::extern_init(long* ptr, long* idx, double* val, long nr, long nc, long nz)
{
  pointr  = ptr;
  rowind  = idx;
  value   = val;

  nrow    = nr;
  ncol    = nc;
  m_nz    = nz;
} 


svdpack::svdpack(long nm, long nz, long sd)
{
  NZMAX   = nz;
  NMAX    = nm;
  SEED    = sd;
}

/**************************************************************
 *                                                            *
 * multiplication of matrix A by vector x, where A is m by n  *
 * (m >> n) and is stored using the Harwell-Boeing compressed *
 * column sparse matrix format.  y stores product vector.     *
 *                                                            *
 **************************************************************/

void svdpack::opa(long m, long n, double *x, double *y)

{
   long end,i,j;

   mxvcount += 1;
   for (i = 0; i < m; i++) y[i] = ZERO;

   for (i = 0; i < n; i++) {
      end = pointr[i+1];
      for (j = pointr[i]; j < end; j++)
	 y[rowind[j]] += value[j] * x[i]; 
   }
   return;
}


/**************************************************************
 *                                                            *
 * multiplication of an n by m matrix A' by a vector X, store *
 * in Y.                                                      *
 *                                                            *
 **************************************************************/

void svdpack::opat(long n, double *x, double *y)
{
   long end,i,j;
   
   mtxvcount += 1;
   for (i = 0; i < n; i++) y[i] = ZERO;

   for (i = 0; i < n; i++) {
      end = pointr[i+1];
      for (j = pointr[i]; j < end; j++)
	 y[i] += value[j] * x[rowind[j]]; 
   }
   return;
}


double svdpack::fsign(double a,double b)
  /************************************************************** 
   * returns |a| if b is positive; else fsign returns -|a|      *
   **************************************************************/ 
{
  
  if ((a>=0.0 && b>=0.0) || (a<0.0 && b<0.0))return(a);
  if  ((a<0.0 && b>=0.0) || (a>=0.0 && b<0.0))return(-a);
}

double svdpack::dmax(double a, double b)
  /************************************************************** 
   * returns the larger of two double precision numbers         *
   **************************************************************/ 
{
  
  if (a > b) return(a);
  else return(b);
}

double svdpack::dmin(double a, double b)
  /************************************************************** 
   * returns the smaller of two double precision numbers        *
   **************************************************************/ 
{
  
  if (a < b) return(a);
  else return(b);
}

long svdpack::imin(long a, long b)
  /************************************************************** 
   * returns the smaller of two integers                        *
   **************************************************************/ 
{
  
  if (a < b) return(a);
  else return(b);
}

long svdpack::imax(long a,long b)
  /************************************************************** 
   * returns the larger of two integers                         *
   **************************************************************/ 
{
  
  if (a > b) return(a);
  else return(b);
}
 
/***********************************************************************
 *                                                                     *
 *                        orthg()                                      *
 *         Gram-Schmidt orthogonalization procedure                    *
 *                                                                     *
 ***********************************************************************/
/***********************************************************************

   Description
   -----------

   The p by n matrix Z stored row-wise in rows f to (f+p-1) of
   array X is reorthogonalized w.r.t. the first f rows of array X.
   The resulting matrix Z is then factored into the product of a
   p by n orthonormal matrix (stored over matrix Z) and a p by p
   upper-triangular matrix (stored in the first p rows and columns 
   of array B).  (Based on orthog from Rutishauser) 


   Parameters
   ----------

   (input)
   p           number of consecutive vectors of array x (stored row-wise)
	       to be orthogonalized
   f           number of rows against which the next p rows are to be
	       orthogonalized
   n           column dimension of x
   x           2-dimensional array whose p rows are to be orthogonalized
	       against its first f rows
   temp        work array


   (output)
   x           output matrix whose f+p rows are orthonormalized
   b           p by p upper-triangular matrix


   Functions called
   --------------

   BLAS         dgemv, ddot, dscal, daxpy, dcopy

 ***********************************************************************/

void svdpack::orthg(long p, long f, long n, double **b, double **x, double *temp)

{
   long fp, k, km1;
   long orig, small;
   double t, s;

   if (!p) return;
   if (f == 0 && p > n) {
      fprintf(stderr,"%s\n",
         "*** ON ENTRY TO ORTHG, MATRIX TO BE ORTHONORMALIZED IS SINGULAR");
      return;
   }
   fp = f + p;

   for (k = f; k < fp; k++) {
      km1 = k - 1;
      orig = TRUE;

      while(TRUE) {
         t = ZERO;

	 if (km1 >= 0) {
	    if (km1 > 0) {
	       dgemv(NTRANSP, k, n, ONE, x, x[k], ZERO, temp);
	       t += ddot(k, temp, 1, temp, 1);
	    }

	    else {
	       temp[0] = ddot(n, x[0], 1, x[k], 1);
	       t += temp[0] * temp[0];
	    }

	    if (orig && km1 >= f) 
               dcopy(k - f, &temp[f], 1, &b[k - f][0], 1); 

            if (km1 > 0) 
	       dgemv(TRANSP, k, n, -ONE, x, temp, ONE, &x[k][0]);
            else
	       daxpy(n, -temp[0], x[0], 1, x[k], 1);
         }

	 if (km1 < 0 || p != 1) {
	    s = ddot(n, x[k], 1, x[k], 1);
	    t += s;
	    if (s > t/CONST) {
	       small = FALSE;
	       s = sqrt(s);
               b[k - f][k - f] = s;
	       if (s != ZERO) s = ONE/s;
	       dscal(n, s, x[k], 1);
	    }
	    else {
	       small = TRUE;
	       orig  = FALSE;
	    }
	 }
	 if (small == FALSE || p == 1) break;
      }
   }
}


/***********************************************************************
 *                                                                     *
 *                         dgemv()                                     *
 * A C-translation of the level 2 BLAS routine DGEMV by Dongarra,      *
 * du Croz, and Hammarling, and Hanson (see LAPACK Users' Guide).      *
 *                                                                     *
 ***********************************************************************/
/***********************************************************************

   Description
   -----------

   dgemv() performs one of the matrix-vector operations

   y := alpha * A * x + beta * y  or  y := alpha * A' * x + beta * y

   where alpha and beta are scalars, X, Y are vectors and A is an
   m by n matrix.

void dgemv(long transa, long m, long n, 
           double alpha, double **a, double *x, double beta, double *y)

   Parameters
   ----------

   (input)
   transa   TRANSP indicates op(A) = A' is to be used in the multiplication
	    NTRANSP indicates op(A) = A is to be used in the multiplication

   m        on entry, m specifies the number of rows of the matrix A.
	    m must be at least zero.  Unchanged upon exit.

   n        on entry, n specifies the number of columns of the matrix A.
	    n must be at least zero.  Unchanged upon exit.

   alpha    a scalar multiplier.  Unchanged upon exit.

   a        matrix A as a 2-dimensional array.  Before entry, the leading
	    m by n part of the array a must contain the matrix A.

   x        linear array of dimension of at least n if transa = NTRANSP
	    and at least m otherwise.

   beta     a scalar multiplier.  When beta is supplied as zero then y
	    need not be set on input.  Unchanged upon exit.

   y        linear array of dimension of at least m if transa = NTRANSP
	    and at leat n otherwise.  Before entry with beta nonzero,
	    the array y must contain the vector y.  On exit, y is 
	    overwritten by the updated vector y.


 ***********************************************************************/

void svdpack::dgemv(long transa, long m, long n, 
           double alpha, double **a, double *x, double beta, double *y)

{
   long info, leny, i, j;
   double temp, *ptrtemp;

   info = 0;
   if      ( transa != TRANSP && transa != NTRANSP ) info = 1;
   else if ( m < 0 ) 				     info = 2;
   else if ( n < 0 )				     info = 3;

   if (info) {
      fprintf(stderr, "%s %1ld %s\n",
      "*** ON ENTRY TO DGEMV, PARAMETER NUMBER",info,"HAD AN ILLEGAL VALUE");
      return;
      //throw new svdpack_error(info);
   }

   if (transa) leny = n;
   else        leny = m;

   if (!m || !n || (alpha == ZERO && beta == ONE))
      return;

   ptrtemp = y; 

   /* form Y := beta * Y */
   if (beta == ZERO) 
      for (i = 0; i < leny; i++) *ptrtemp++ = ZERO;
   else if (beta != ONE) 
      for (i = 0; i < leny; i++) *ptrtemp++ *= beta;

   if (alpha == ZERO) return;

   switch(transa) {