tms1.c

求矩阵奇异分解svd算法

  /*for i = n step -1 until 2 do*/

  for (ii=3;ii<n+2-offset;ii++)
   {
     i=n+2-ii;
     l=i-1;
     h=ZERO; 
     scale=ZERO;

    /*scale row (algol tol then not needed)*/
     if (l>=1)
       for (k=offset;k<=l;k++)
        {
         scale+= fabs(d[k]);
        }
	
    if ((scale==ZERO)||(l<1))
     {
      e[i]=d[l];
      for (j=offset;j<=l;j++)
          {
            d[j]=z[l][j];
            z[i][j]=ZERO;
            z[j][i]=ZERO;
          }
     }
   else                   /*scale <> ZERO */
     {
       for (k=offset;k<=l;k++)
        {
         d[k]=d[k]/scale;
         h+=d[k]*d[k];
        }


       f=d[l];
       g=-fsign(sqrt(h), f);
       e[i]=scale * g;
       h-=f*g;
       d[l]=f-g;
   
       /* form A*u */
  
       for (j=offset; j<=l; j++)
          e[j]=ZERO;
          
          for (j=offset;j<=l;j++)
            {
             f=d[j];
             z[j][i]=f;
             g= e[j] + z[j][j] * f;
             
             jp1= j + 1;
   
             if (l >= jp1) 
                 {
                  for (k=jp1; k<=l; k++)
                   {
                     g+= z[k][j] * d[k];
                     e[k] += z[k][j] * f;
                   }
                 };
             e[j]=g;
           }

       /* form P */
 
       f= ZERO;
 
       for (j=offset; j<=l; j++)
        {
          e[j]=e[j]/h;
          f+= e[j] * d[j];
        }

       hh= f/ (h+h);
  
       /* form Q */
  
      for (j=offset; j<=l; j++)
       e[j] -= hh * d[j];

      /* form reduced A */

      for (j=offset; j<=l; j++)
       {
         f= d[j];
         g = e[j];

         for (k=j; k<=l; k++)
          z[k][j]= z[k][j] - f * e[k] - g * d[k];

         d[j]=z[l][j];
         z[i][j]=ZERO;
       }
    }  /* end scale <> zero */

    d[i]=h;
   }   /* end for ii */
   /*accumulation of transformation matrices */

   for (i=offset + 1;i<n;i++)
    {
     l=i-1;
     z[n-1][l] = z[l][l];
     z[l][l] = ONE;
     h=d[i];

     if (h != ZERO) 
       {
        for (k=offset; k<=l; k++)
          d[k]= z[k][i]/h;

        for (j=offset; j<=l; j++)
         {
           g= ZERO;
           
           for (k=offset;k<=l; k++)
            g+= z[k][i]*z[k][j];

	   for (k=offset;k<=l;k++)
            z[k][j] -= g * d[k];
         }
       }
       for (k=offset;k<=l;k++) z[k][i]=ZERO;
     }
  
     for (i=offset;i<n;i++)
       {
        d[i]=z[n-1][i];
        z[n-1][i]=ZERO;
       }
     z[n-1][n-1]=ONE;
     e[0]=ZERO;

/*preparation for tql2.c.. reorder e[]*/
for (i=1+offset;i<n;i++) e[i-1]=e[i]; 

/*preparation for tql2.c.. z has to be transposed for 
  tql2 to give correct eigenvectors */
for (ii=offset; ii<n; ii++)
 for (jj=ii; jj<n; jj++)
 {
   tmp=z[ii][jj];
  z[ii][jj]=z[jj][ii];
  z[jj][ii]=tmp;
 }

     return;
}
           
        
#include <stdio.h>
#include <math.h>

extern long ncol,nrow;
extern long *pointr,*rowind;
extern double *value,alpha;

void opb(long n, double *x, double *y, double shift)
{
/* multiplication of a shifted 2-cylic matrix c by a vector x , where

       [D    A]  
   c = [      ]
       [A'   D] , where A is nrow by ncol (nrow>>ncol),
                  and D = (alpha-shift)*I , alpha an upper bound
                  for the largest singular value of A, and 
                  shift is the approximate singular value of A.

   Hence, C is of order N=nrow+ncol (y stores product vector) */

  long i,j,upper;

  for(i=0; i<n; i++)
     y[i]=(alpha-shift)*x[i];

  for(i=0; i<ncol; i++) {
     upper = pointr[i+1];
     for(j=pointr[i]; j<upper; j++)
        y[rowind[j]] += value[j]*x[nrow+i];
  }

  for(i=nrow; i<n; i++) {
     upper = pointr[i-nrow+1];
     for(j=pointr[i-nrow]; j<upper; j++)
        y[i] += value[j]*x[rowind[j]];
  }

  return;
}
/************************************************************** 
 * Constant times a vector plus a vector     		      *
 * Based on Fortran-77 routine from Linpack by J. Dongarra    *
 **************************************************************/ 

void daxpy (long n,double da,double *dx,long incx,double *dy,long incy)

{
   long i;

   if (n <= 0 || incx == 0 || incy == 0 || da == 0.0) return;
   if (incx == 1 && incy == 1) 
      for (i=0; i < n; i++) {
	 *dy += da * (*dx++);
	 dy++;
      }
   else {
      if (incx < 0) dx += (-n+1) * incx;
      if (incy < 0) dy += (-n+1) * incy;
      for (i=0; i < n; i++) {
         *dy += da * (*dx);
         dx += incx;
         dy += incy;
      }
   }
   return;
}
/************************************************************** 
 * Function forms the dot product of two vectors.      	      *
 * Based on Fortran-77 routine from Linpack by J. Dongarra    *
 **************************************************************/ 

double ddot(long n,double *dx,long incx,double *dy,long incy)

{
   long i;
   double dot_product;

   if (n <= 0 || incx == 0 || incy == 0) return(0.0);
   dot_product = 0.0;
   if (incx == 1 && incy == 1) 
      for (i=0; i < n; i++) dot_product += (*dx++) * (*dy++);
   else {
      if (incx < 0) dx += (-n+1) * incx;
      if (incy < 0) dy += (-n+1) * incy;
      for (i=0; i < n; i++) {
         dot_product += (*dx) * (*dy);
         dx += incx;
         dy += incy;
      }
   }
   return(dot_product);
}
#include <stdio.h>

void 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;
}
#include <stdio.h>
#include <math.h>
#include "tmsc.h"

void 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;
}
#include <stdio.h>
#include <math.h>
#include "tmsc.h"

double dasum(long, double *, long);

void 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;
}
#include <stdio.h>
#include <math.h>
#include "tmsc.h"

extern double *z,*v2,*r,*r2,*pp,*z2;

void daxpy(long, double, double *,long, double *, long);
double ddot(long, double *,long, double *, long);
void dscal(long, double, double *,long);
void dgemv2(long, long, long, double, double **, long, long,
        double *, double, double *);
void opb(long, double *, double *, double);

void 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, temp;
  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 */

  opb(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++) {

     opb(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;
}
#include <stdio.h>
#include <math.h>
#include "tmsc.h"

extern double *z,*v2,*r,*r2,*pp,*z2;

void daxpy