sis2.c

求矩阵奇异分解svd算法

 BLAS:  ddot, dscal, daxpy
 EISP:  tred2, tql2, pythag
 MISC:  opb, inf, dmax, imax
 
*****************************************************************/
void ritzit( long n, long kp, long km, double eps, void (*opb) (long , double *, double * ), void (*inf)(long , long , long , double *, long ),long kem, double **x, double *d, double *f, double *cx, double *u, long *imem)

{
  long i, j, l, i1, l1, size, kg, kh, kz, kz1, kz2, ks, m, m1;
  long  jj, ii, ig, ip, ik, jp, k, flag1 ;  
  double  *w[NMAX], *rq, *b[NSIG], *s, *tptr, TmpRes;
  double xks, ee2, e, e2, e1, xkm, xk, xm1, t;

  /*******************************************************************
   * allocate memory for temporary storage arrays used in ritzit only*
   *                                                                 *
   *            rq - (kp)                                            *
   *             w - (NSIG * n)                                      *
   *             b - (NSIG * kp)                                     *
   *             s - (kp)                                            *
   *******************************************************************/

   size =    kp * (NSIG +2) + NSIG * n;
   if (!(rq=(double *) malloc(size * sizeof(double)))){
      perror("MALLOC FAILED IN RITZIT()");
      exit(errno);
   }
   tptr= rq + kp;
   for (i=0; i< NSIG; i++){
    w[i]= tptr;
    tptr+=n;
   }

   /*w=tptr; tptr+=(n * kp);*/
   for (i=0; i<NSIG; i++) {
    b[i]=tptr;
    tptr+=kp;
   };
   s=tptr;

   /*finished allocating memory*/


  *imem += size;
  *imem = sizeof(double) * (*imem);

  ee2 = ONE + 1.0e-1 * eps;
  e   = ZERO;
  kg  = -1;
  kh  = -1;
  kz  = 1367;
  kz1 = 0;
  kz2 = 0;
  ks  = 0;
  m   = 1;

  for (l=0; l<kp; l++){
    f[l] = 4.0e+0;
   cx[l] = ZERO;
   rq[l] = ZERO;
 } 
 
  if (km >= 0)    /*generate random initial iteration vectors*/
      for (j=0; j<kp; j++)
       for (l=0; l<n; l++){
         kz = (3125 * kz) % 65536;  
        x[j][l] = (double) (kz - 32768);
      }
  km =abs(km);
  l1 = 1;
  i1 = 1;
  jp=kp;
flag1=0;
#define NRow kp
#define NCol n

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

  ig = 0;
  ip = kp - 1;
  
 /*statement 70  from original RITZIT begins here*/
while (1){

  flag1=1;  /* so that jacobi step is done at least once */

  while (flag1){
        /* jacobi step modified */
        for(k=ig; k< kp; k++){
          opb(n, &x[k][0], &w[0][0]);
          for (j=0; j<n; j++) x[k][j]=w[0][j];
         }
        l1=ig + 1;
        jp=kp;
        flag1=0; /*flag1 is set only if re-orthog needs to be done*/

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

       if ( ks<=0 ) {
         /* measures against unhappy choice of initial vectors*/
         for (k=0;k<kp; k++)
           if (b[k][k] ==ZERO)
              for (j=0;j<n;j++){
                 kz= (3125 * kz) % 65536;
                 x[k][j] = (double)(kz-32768);
                 flag1=1;
               }
       }
     if (flag1){
       l1=1;
       ks=1; /*we dont want to re-initialize x[][] again*/

   /* 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 */
      } 
    } /*end while flag1 */

  for (k= ig; k<kp; k++)
    for (l=k; l<kp; l++){
      t = ZERO;

      for (i=l; i<kp; i++)
        t+= b[k][i] * b[l][i];

      /* negate matrix to reverse eigenvalue ordering */
      b[k][l] = -t;
    }
  j=kp - kg - 1;

  tred2(ig, kp, b, d, u, b);
  ii=tql2(ig, kp, d, u, b);


  for (k=ig; k< kp; k++)
   d[k]=sqrt(dmax(-d[k], ZERO));
  

 for (j=ig; j<kp; j++)
  for (k=0; k<n; k++) 
   w[j][k]=ZERO;

   for (j=ig; j<kp; j++)
     for (l=ig; l<kp; l++){
       TmpRes=b[j][l];
       for (k=0; k<n; k++)
         w[j][k] += TmpRes  * x[l][k];}   /* w is going to be transposed as 
                                         compared to the fortran version */
     

  for (j=ig; j<kp; j++)
   for (k=0; k<n; k++)
      x[j][k]=w[j][k]; 

  xks=(double)(++ks);
  if (d[kp-1] > e) e=d[kp-1 ];

  /* randomization */
  if (kz1<3) {
    for (j=0; j<n; j++) {
       kz= (3125 * kz) % 65536;
       x[kp-1][j] = (double) (kz -32768);
    }
    l1=kp;
    jp=kp;

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

  }
  
  /*compute control quantities cx */
  for (k=ig; k<ip;k++){
    t=(d[k] - e)*(d[k] + e);
    if (t <= ZERO) cx[k]=ZERO; 
    else 
           if (e ==ZERO) cx[k] = 1.0e+3 + log(d[k]);
           else  cx[k]= log( (d[k]+sqrt(t)) / e);
  }
    
    /*acceptance test for eigenvalues including adjustment of 
      kem and kh such that 
           d[kem] > e, 
           d[kh]  > e, and,
           d[kem] does not oscillate strongly */

    for (k=ig; k<kem; k++)
       if ((d[k] <= e) ||
           ((kz1>1) && (d[k] <= 9.99e-1 * rq[k]))){
          kem=k-1;
          break;
       }

   if (!kem) break; 

   k=kh+1;

   while ((d[k] != ZERO) && ( d[k] <= ee2 * rq[k])){
      kh=k++;
   }

   do {
    if (d[k] <= e) kh= k- 1;
    --k;
   }
   while (k > kem-1);

  /*acceptance test for eigenvectors */

  l= kg;
  e2= ZERO;
 
  for (k=ig; k<ip; k++){
     if (k == l+1){
       /*check for nested eigenvalues */
       l=k;
       l1=k;
       if (k != ip){
          ik= k + 1;
          s[0] = 5.0e-1 / xks;
          t= ONE / (double)(ks * m);
          for (j=ik; j<ip; j++){
            if ((cx[j] * (cx[j] + s[0]) + t) <= (cx[j-1] * cx[j-1]))
                break;
            l=j;
          }
       }
     }
     if (l > kh) {l = l1; break;}
  
         opb(n, &x[k][0], &w[0][0]);
         s[0]= ZERO;
      
         for (j=0; j<=l; j++)
            if (fabs(d[j] - d[k]) < 1.0e-2 * d[k]){
                t=ZERO;
 
                for (i=0; i<n; i++) 
                    t+= w[0][i] * x[j][i];

                for (i=0; i<n; i++) {
                    w[0][i] =w[0][i] - t * x[j][i];
                }

                s[0]+= t*t;
            }
  
         t=ZERO;

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

         if (s[0] != ZERO) t=sqrt(t/(s[0] + t)); 
         else t=ONE;
   
         if (t > e2) e2=t;

         if (k==l){
             /* test for acceptance of group of eigenvectors*/
             if ((l>=kem-1) &&
                 (d[kem] * f[kem] < eps * (d[kem] -e)))
               kg=kem;

             if (e2 < f[l])
               for (j=l1; j<=l; j++) f[j]=e2;

             if ((l <= kem) &&
                 (d[l] * f[l] < eps * (d[l]-e)))
               kg=l;

             ig=kg+1;
         }
}

  /* statements 570 to 660  from original RITZIT*/

   if (e<= 4.0e-2 * d[0]) {