sis2.c

svd 算法代码 This directory contains instrumented SVDPACKC Version 1.0 (ANSI-C) programs for compiling

     m=1;
     k=1;
   }
   else {
     e2=2.0e0/e;
     e1=5.1e-1 * e2;
     k= 2 * imax((long)(4.0e0/cx[0]), 1);
     if (m>k) m=k;
   }
   /*reduce kem if convergence would be too slow */
   xkm=(double)km;
   if ((f[kem-1] != ZERO) &&
       (xks < 9.0e-1 *xkm)){
         xk=(double)k;
         s[0]=xk * cx[kem-1];
         if (s[0] < 5.0e-2) t=5.0e-1 * s[0] *cx[kem-1];
         else t=cx[kem-1] + log(5.0e-1 * (ONE + exp(-2.0e0 * s[0])))/xk;
         s[0]=log(d[kem-1] * f[kem-1]/(eps * (d[kem-1] -e)))/t;
         if ((xkm - xks) * xkm < s[0] * xks) kem--;
   }
   inf(ks,kg,kh,f,m);

  if ((kg >= kem-1) ||
      (ks >= km))      break;
 
  for (k=ig; k<kp; k++)  rq[k] = d[k];

  do {
      /*statements 680-700 */
      if (ks + m > km) {
           kz2=-1;
           if (m >1) m = 2* ((km -ks +1)/2);
      }
      else
          m1=m;

      /*shortcut last intermediate block if all error quantities f are
        sufficiently small */
      if (l >= kem){
        s[0]= d[kem-1] * f[kem-1]/(eps *(d[kem-1] -e));
        t= s[0] * s[0] - ONE;
        if (t <= ZERO) break;
        s[0] = log(s[0] + sqrt(t))/(cx[kem-1] -cx[kh+1]);
        m1=2 * (long)(5.0e-1 * s[0] + 1.01e0);
        if (m1<=m) kz2=-1;
        else m1=m;
     }
     xm1=(double) m1;
     
                                   /*chebyshev iteration */
     if (m==1)
       for (k=ig; k<kp; k++){
       opb(n, &x[k][0], &w[0][0]);
       for (i=0; i<n; i++) 
           x[k][i] = w[0][i];
       }
     else                                /*degree != ONE */
       for (k=ig; k<kp; k++){
          opb(n, &x[k][0], &w[0][0]);

          for (i=0; i<n; i++) u[i]=e1 * w[0][i];

          opb(n, u, &w[0][0]);

          for (i=0; i<n; i++) x[k][i]= e2 *w[0][i] - x[k][i];

          if (m1>=4)
               for (j=3; j<m1; j+=2){
                 opb(n, &x[k][0], &w[0][0]);
 
                 for (i=0; i<n; i++)
                    u[i]=e2 * w[0][i] - u[i];

                 opb(n, u, &w[0][0]);

                 for (i=0; i<n; i++)
                     x[k][i] = e2 * w[0][i] - x[k][i];
                }
       }

       
       l1=ig+1;

    /* extend orthonormalization to all kp rows of x

     Variables used here:
 
     NCol                     : number of columns of x
     NRow                     : number of rows of x
     ii, ik, jp, k, l1        : indexing integers
     t, TmpRes                : double precision storage.


     at the end of this section of code, 

     x   : transpose(Q)
     b   : R

     where X = QR is the orthogonal factorization of the matrix X
     stored in array x

     invokes: ddot, daxpy, dscal (BLAS)
              sqrt               (math.h)
  */

    ik = l1 -1; 
  for (i=0; i< jp; i++){
    for (k=l1-1; k<jp; k++) b[i][k]= ZERO;
   
    if (i >= l1-1){
      ik=i+1;
      b[i][i]=sqrt(ddot(NCol, &x[i][0], 1 ,  &x[i][0], 1));
      t=ZERO;
      if (b[i][i] != ZERO) t=ONE / b[i][i];
      dscal(NCol, t, &(x[i][0]),  1);
     }

   for (ii=ik; ii< NRow; ii++){
     TmpRes=ZERO;
     for (jj=0; jj<NCol; jj++)
       TmpRes+= x[ii][jj] * x[i][jj];
     b[i][ii]=TmpRes;
  }
  
    for (k=ik; k<jp;k++) 
        daxpy(NCol, -b[i][k], &x[i][0], 1, &x[k][0], 1);

  } /* end for i */

      /*discounting the error quantities F */
      if (kg!=kh){
         if (m ==1) 
              for (k=ig; k<=kh; k++) f[k]= f[k] * (d[kh+1]/d[k]);
         else {
               t=exp(-xm1 * cx[kh+1]);
               for (k=ig; k<=kh; k++){
                   s[0]=exp(-xm1 * (cx[k]-cx[kh+1]));
                   f[k] = s[0] * f[k] * (ONE + t*t)/(ONE + (s[0] *t)*(s[0]*t));
               }
         } 
      }
      ks=ks+m1;
      kz2=kz2 - m1;
     
   } /*possible repetition of intermediate steps*/
   while (kz2>=0);

   kz1++;
   kz2 = 2 * kz1;
   m = 2 * m; 

}  /* end while kz2<=0 ... go to 70*/

/*statement 900 of original RITZIT program begins here*/

kem = kg;
l1=1; 
jp=kp-1;

  /* extend orthonormalization to all kp rows of x 

     Variables used here:
 
     NCol                     : number of columns of x
     NRow                     : number of rows of x
     ii, ik, jp, k, l1        : indexing integers
     t, TmpRes                : double precision storage.


     at the end of this section of code, 

     x   : transpose(Q)
     b   : R

     where X = QR is the orthogonal factorization of the matrix X
     stored in array x

     invokes: ddot, daxpy, dscal (BLAS)
              sqrt               (math.h)
  */

    ik = l1 -1; 
  for (i=0; i< jp; i++){
    for (k=l1-1; k<jp; k++) b[i][k]= ZERO;
   
    if (i >= l1-1){
      ik=i+1;
      b[i][i]=sqrt(ddot(NCol, &x[i][0], 1 ,  &x[i][0], 1));
      t=ZERO;
      if (b[i][i] != ZERO) t=ONE / b[i][i];
      dscal(NCol, t, &(x[i][0]),  1);
     }

   for (ii=ik; ii< NRow; ii++){
     TmpRes=ZERO;
     for (jj=0; jj<NCol; jj++)
       TmpRes+= x[ii][jj] * x[i][jj];
     b[i][ii]=TmpRes;
  }
  
    for (k=ik; k<jp;k++) 
        daxpy(NCol, -b[i][k], &x[i][0], 1, &x[k][0], 1);

  } /* end for i */

/*statements 920 up to last card of the original RITZIT program */
 
for (k=0; k<ip; k++)
  for (i=k; i<ip; i++) b[k][i]=ZERO;

for (k=0; k<ip; k++) {
   opb(n, &x[k][0], &w[0][0]);

   for (i=0; i<=k; i++)
        for (l=0; l<n; l++)
          b[i][k]=b[i][k] - x[i][l] * w[0][l];
   }

tred2(0, ip, b, d, u, b);
tql2(0, ip, d, u, b);

/*reordering of eigenvalues and eigenvectors according to the magnitudes 
  of the former */

for (i=0; i< ip; i++)
  if (i!=ip-1){ 
       k=i;
       t=d[i];
       ii=i+1;
       for (j=ii; j<ip; j++)
         if (fabs(d[j]) > fabs(t)){
             k=j; 
             t=d[j];
         }
       if (k!=i) {
          d[k] = d[i];
          d[i] = t;

          for (j=0; j<ip; j++){
            s[j]= b[i][j];
            b[i][j] = b[k][j];
            b[k][j] = s[j];
          }
      }
   d[i]=-d[i];
   }

  for (i=0; i<ip; i++)
   for(j=0; j<n; j++)
    w[i][j]=ZERO;

  for (i=0; i<ip; i++)
   for (k=0; k<ip; k++){
    TmpRes=b[i][k];
    for (j=0; j<n; j++)
      w[i][j]+=TmpRes * x[k][j];
   }

     for (i=0; i<ip; i++)
        for (j=0; j<n; j++)
           x[i][j] = w[i][j];

     d[kp-1]=e;

  /* free memory at the end of ritzit*/
  free(rq);
  return;

}  /*end ritzit*/

#include "sisc.h"

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

/**************************************************************
 * multiplication of matrix A by a vector x, where            *
 *                                                            *
 * A is nrow by ncol (nrow >> ncol)                           *
 * y stores product vector                                    *
 **************************************************************/
void opa(long n,double *x, double *y)

{
   long i, j, end;

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

/* multiply by sparse C */
   for (i = 0; i < ncol; i++) {
      end = pointr[i+1];
      for (j = pointr[i]; j < end; j++)
	y[rowind[j]] += value[j] * x[i];
   }

   return;
}
#include "sisc.h"

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

/************************************************************** 
 * multiplication of matrix B by a vector x, where            *
 *							      *
 * B =  A'A, where A is nrow by ncol (nrow >> ncol)           *
 * Hence, B is of order n:=ncol                               *
 * y stores product vector                                    *
 **************************************************************/ 
void opb(long n,double *x, double *y)

{
   long i, j, end;
   double *ztemp;
 

  ztemp=(double *) malloc(nrow * sizeof(double));
   mxvcount += 2;

   for (i = 0; i < ncol; i++) y[i] = ZERO;

   for (i=0; i<nrow; i++) ztemp[i] = ZERO;

/* multiply by sparse C */
   for (i = 0; i < ncol; i++) {
      end = pointr[i+1];
      for (j = pointr[i]; j < end; j++)
	ztemp[rowind[j]] += value[j] * x[i];
   }

/*multiply by sparse C' */
   for (i=0;i<ncol; i++){
      end=pointr[i+1];
      for (j=pointr[i]; j< end; j++)
         y[i] += value[j] * ztemp[rowind[j]];
   }

   free(ztemp);
   return;
}
#include <math.h>
double fsign(double, double);
/***********************************************************************
 *                                                                     *
 *                              tred2()                                *
 *                                                                     *
 ***********************************************************************/
/***********************************************************************

  Description
  -----------
  
  tred2() is a translation of the algol procedure TRED2, Num. Math. 11, 
  181-195 (1968) by Martin, Reinsch, and Wikinson.  Handbook for Auto.
  Comp., Vol. II- Linear Algebra, 212-226 (1971)

  This subroutine reduces a real symmetric matrix to a symmetric
  tridiagonal matrix using and accumulating orthogonal similarity
  transformations.

  Arguments
  ---------

  (input)
  offset index of the leading element of the matrix to be
         tridiagonalized. The matrix tridiagonalized should be 
         stored in a[offset:n-1, offset:n-1]

  n	 order of the matrix

  a	 contains the real symmetric input matrix. Only the upper
	 triangle of the matrix need be supplied

  (output)
  d	 contains the diagonal elements of the tridiagonal matrix.
  
  e	 contains the subdiagonal elements of the tridiagonal matrix
	 in its first n-1 positions.

  z	 contains the orthogonal transformation matrix produced in the
	 reduction.

  a and z may coincide. If distinct, a is unaltered.

  Functions used:
  UTILITY: fsign

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

void tred2(long offset, long n, double **a, double *d, double *e, double **z)
{
 long jj,ii,i,j,k,l, jp1;
 double h, scale, f, g,  hh, tmp;

 

 for (i=offset;i<n;i++) 
  { 
   for (j=i;j<n;j++)
     {
      z[j][i]=a[i][j];   /*fix this later?.. the rest of the routine 
                           assumes that z has the lower triangular part
                           of the symmetric matrix */
     }
   d[i]=a[i][n-1];
  }


  if (n==1) 
   {
    for (i=offset;i<n;i++)
     {
       d[i]=z[n-1][i];
       z[n-1][i]=ZERO;
     }
    z[n-1][n-1]=ONE;
    e[1]=ZERO;
    return;
   }

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