tms.cc

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

// File: tms.cc -- Implements class tms
// Author: Suvrit Sra
// Date: 15, nov 2003

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

#include <cmath>
#include "tms.h"

using namespace ssvd;
double tms::epslon(double x)
{
  /* Estimate unit roundoff in quantities of size x 
     
     This function should work properly on all systems
     satisfying the following two assumptions,
     1. The base used in representiong floating point
        numbers is not a power of three.
     2. The quantity a in statement 10 is represented to 
        the accuracy used in floating point variables
        that are stored in memory.
     The statement number 10 and the while are intended 
     to force optimizing compilers to generate code
     satisfying assumption 2.
     Under these assumptions, it should be true that,
         A is not exactly equal to four-thirds,
         B has a zero for its last bit or digit,
         C is not exactly equal to one,
         EPS measures the separation of 1.0 from
             the next larger floating point number.
     This routine is based on Eispack. The developers of
     Eispack would appreciate being informed about any
     systems where these assumptions do not hold. */

  double  a,b,c,eps;

  a = 4.0/3.0;
  do {
     b = a - ONE; /* statement 10 */
     c = b + b + b;
     eps = fabs(c -ONE);
  }
  while (eps == ZERO);  
  eps = eps*fabs(x);
  return(eps);
}


void tms::disk(long n, double *sig, double *rad, long *csize, long *clus)
{
  /* monitor separation of gershgorin disks */

  double radi,radipi,tmp;
  long   cptr,i,k,flag,upper;

  /* assume ordering: sig[0] <= sig[1] <= ...<= sig[n-1] for sig 
     array in tsvd1.
     rad := radii for approximate e-values  

     csize array indicates location and size of clusters:
     csize[i]=k  (k != 0) gives cluster of size k with first disk
     centered at sig[i]. */

  upper = n-1;
  for(i=0; i<upper; i++) {
     radi = rad[i];
     radipi = rad[i+1];
     tmp = (radi+radipi) - (sig[i+1] - sig[i]);
     if(tmp<=0.0) clus[i] = 1;
     else  clus[i] = 0;
  }
  /* clustering vector filled, now locate clusters and their size */

  for(i=0; i<n; i++) csize[i] = 1;

  for(i=0; i<upper; i++) {
     cptr = i;
     flag = 1;
     k = i-1;
     while((flag) && (k>=0)) {
        if(csize[k])  flag = 0;
        else  cptr=k;
        k -= 1;
     }
   if(!clus[i]) {
     if(csize[i]) csize[i] += 1;
     else csize[cptr-1] += 1;
     csize[i+1] -= 1;
   }
  }
return;
}


void tms::isol(long i, double resid, double tol, long *init, double *ireso,
          double *creso, long s)
{
  /* monitor residual reduction for isolated singuar value
     approximations.
     Assume ordering : sig[0]<=sig[1]<=...<=sig[s-1] for sig 
     array in tsvd1. 
     Resid : 2-norm of B-eigenpair residual   */

  if((ireso[i]<ZERO) || (creso[i]>=ZERO)) {

      if(creso[i]<ZERO) ireso[i]=resid;
      else {
           ireso[i]=resid;
           creso[i]=-ONE;
      }
    }
  else 
       if(resid<=tol*ireso[i]) init[i] = TRUE;

  return;
}

void tms::clus(long i, long size, double *resid, double tol, long *init,
          double *ireso, double *creso, double *tmp, long s)
{
  /* Monitor residual reduction for clustered singular value 
     approxmations.
     Assume ordering : sig[0] <= sig[1] <=...<= sig[n-1] for 
     sig array in tsvd1.

     resid = 2-norm of B-eigenpair residuals. 
     i     = first disk in cluster.
     size  = number of disks in cluster. */


  double error;
  long   k,upper;

  for(k=0; k<size; k++) tmp[k] = resid[k]*resid[k];
  error = sqrt(dasum(size, tmp, 1));
  if(creso[i] < ZERO) { 
    if(ireso[i] < ZERO) creso[i] = error;
    else {
         creso[i] = error;
         ireso[i] = -ONE;
    }
  }
  else
     if(error <= (tol*creso[i])) {
       upper = i+size;
       for(k=i; k<upper; k++) init[k] = TRUE;
     }
  return;
}

void tms::cgts(long n, long left, double *w, double **v, long *cgiter,
          double sig, double sigold, double sigmax, double shift,
          long *kount, double eps)
{
  /* cg for independent systems in trace min. optimization step.
   (shift incorporated)
   v stores current orthog basis for r[y]. */

  double denom, bound, a, a0, b, error, rnorm, rnorm0, vnorm;
  long   maxi, i, j, k, ii;

  maxi = n;
  *kount = 0;
  if(sig!=shift) {
    bound = (sig-shift)/(sigmax-shift);
    bound=bound*bound;
  }
  else {
       bound = (sigold-shift)/(sigmax-shift);
       bound=bound*bound;
  }
  /* w0=w by definition , get first redidual residual via orthog.
     projector */

  myopb(n, w, z, shift);
  *kount += 1;
  dgemv2(TRANSP, n, left, ONE, v, 0, 0, z, ZERO, r);
  dgemv2(NTRANSP, n, left, ONE, v, 0, 0, r, ZERO, v2);
  for(i=0; i<n; i++) r[i] = z[i]-v2[i];
  rnorm0 = sqrt(ddot(n,r,1,r,1));
  for(i=0; i<n; i++) pp[i] = r[i];

  /* main iteration loop 30*/

  for( ii=0; ii<maxi; ii++) {

     myopb(n, pp, v2, shift);
     *kount += 1;
     denom = ddot(n,pp,1,v2,1);
     if(denom<=ZERO){
       return;
     }
     a = ddot(n, r, 1, r, 1) / denom;
     if(ii==0) a0 = a;
     daxpy(n, -a, pp, 1, w, 1);
     for(j=0; j<n; j++) r2[j] = r[j];
     dgemv2(TRANSP, n, left, ONE, v, 0, 0, v2, ZERO, z);
     dgemv2(NTRANSP, n, left, ONE, v, 0, 0, z, ZERO, z2);
     for(i=0; i<n; i++) v2[i] -= z2[i];
     daxpy(n, -a, v2, 1, r2, 1);
     rnorm = sqrt(ddot(n,r2,1,r2,1));
     for(k=0; k<n; k++) v2[k] = pp[k];
     dscal(n, a, v2, 1);
     vnorm = sqrt(ddot(n,v2,1,v2,1));

  /* early termination code: */

     error = fabs(a*rnorm*rnorm/(a0*rnorm0*rnorm0));
     if((error<=bound) || (vnorm <= eps))  {
       *cgiter += ii+1;
       return;
       }
     else if (ii==maxi-1) {
          *cgiter += maxi;
           printf("cgts failed to converge in  %ld  iterations\n",maxi);
           return;
          }

     for(j=0; j<n; j++) v2[j] = r2[j];
     b=ddot(n, r2, 1, r2, 1) / ddot(n, r, 1, r, 1);
     daxpy(n, b, pp, 1, v2, 1);
     for(j=0; j<n; j++) pp[j] = v2[j];
     for(j=0; j<n; j++) r[j] = r2[j];
  }
return;
}



void tms::cgt(long n, long left, double *w, double **v, long *cgiter, 
         double sig, double sigmax, long *kount, double eps)
{
  /* cg for independent systems in trace min. optimization step.
     v stores current orthog basis for r[y]. */

  double denom, bound, a, a0, b, error, rnorm, rnorm0, vnorm;
  long   maxi, i, j, k, ii;

  maxi = n;
  *kount = 0;
  bound = sig/sigmax ;
  bound = bound*bound;

  /* w0=w by definition , get first residual via orthog, projector */

  myopb(n, w, z, ZERO);
  *kount += 1;
  dgemv2(TRANSP, n, left, ONE, v, 0, 0, z, ZERO, r);
  dgemv2(NTRANSP, n, left, ONE, v, 0, 0, r, ZERO, v2);
  for(i=0; i<n; i++) r[i] = z[i]-v2[i];
  rnorm0 = sqrt(ddot(n,r,1,r,1));
  for(i=0; i<n; i++) pp[i] = r[i];

  /* main iteration loop */

  for( ii=0; ii<maxi; ii++) {

     myopb(n, pp, v2, ZERO);
     *kount += 1;
     denom = ddot(n, pp, 1, v2, 1);
     if(denom <= ZERO) {
       return;
     }
     a = ddot(n, r, 1, r, 1)/denom;
     if(ii==0) a0 = a;
     daxpy(n, -a, pp, 1, w, 1);
     for(j=0; j<n; j++) r2[j] = r[j];
     dgemv2(TRANSP, n, left, ONE, v, 0, 0, v2, ZERO, z);
     dgemv2(NTRANSP, n, left, ONE, v, 0, 0, z, ZERO, z2);
     for(i=0; i<n; i++) v2[i] -= z2[i];
     daxpy(n, -a, v2, 1, r2, 1);
     rnorm = sqrt(ddot(n, r2, 1,r2,1));
     for(k=0; k<n; k++) v2[k] = pp[k];
     dscal(n, a, v2, 1);
     vnorm = sqrt(ddot(n,v2,1,v2,1));

    /* early termination code: */

     error = fabs(a*rnorm*rnorm/(a0*rnorm0*rnorm0));
     if((error <= bound) || (vnorm <= eps))  {
       *cgiter += ii+1;
       return;
     }
     else if(ii==maxi-1) {
           *cgiter += maxi;
           printf("cgt failed to converge in  %ld   iterations\n",maxi);
           return;
          }
   
     for(j=0; j<n; j++) v2[j] = r2[j];
     b=ddot(n, r2, 1, r2, 1)/ddot(n, r, 1, r, 1);
     daxpy(n, b, pp, 1, v2, 1);
     for(j=0; j<n; j++) pp[j]=v2[j];
     for(j=0; j<n; j++) r[j]=r2[j];
  }
  return;
}


void tms::porth(long p, long f, long n, double **x, double *tmp,
           double *tmp1, double *tmp2, double *tmp3,
           double *tmp4, double e, long degree, double alpha)
/* 
     C translation of blas2  Gram-Schmidt orthogonalization procedure 
     for polynomial accelerated trace minimization (job=2)
     the n by p matrix z stored in columns f+1 to f+p of
     array x is reorthogonalized w.r.t. the first f columns of
     array x.  the resulting matrix z contains the desired P-orthogonal
     columns, where P is a polynomial in the matrix b from tsvd2.
     (based on orthog from Rutishauser)  */

{
long fpp,k,km1,j;

double s,t;

  if(!p) return;
  fpp = f+p;

  for(k=f; k<fpp; k++) {
     km1 = k-1;
L10:
    t = 0.0;
    if(km1<0) goto L25;
    pmul(n,x[k],tmp1,tmp2,tmp3,tmp4,e,degree,alpha);

  if(km1>0) {
    dgemv2(TRANSP,n,km1+1,1.0,x,0,0,tmp1,0.0,tmp);
    t += ddot(km1+1,tmp,1,tmp,1);
  }
  else {
     tmp[0] = ddot(n,x[0],1,tmp1,1);
     t += tmp[0]*tmp[0];
  }

  if(km1>0) dgemv2(NTRANSP,n,km1+1,-1.0,x,0,0,tmp,1.0,x[k]);
  else daxpy(n,-tmp[0],x[0],1,x[k],1);
  if(p==0) goto L50;
L25:
  pmul(n,x[k],tmp1,tmp2,tmp3,tmp4,e,degree,alpha);
  s = ddot(n,x[k],1,tmp1,1);