s_sis1.cc

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

  } /* 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;}
  
         myopb(n, &x[k][0], &w[0][0], alpha);
         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]) {
     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--;
   }
   intros(fptr, 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++){
       myopb(n, &x[k][0], &w[0][0], alpha);
       for (i=0; i<n; i++) 
           x[k][i] = w[0][i];
       }
     else                                /*degree != ONE */
       for (k=ig; k<kp; k++){
          myopb(n, &x[k][0], &w[0][0], alpha);

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

          myopb(n, u, &w[0][0], alpha);

          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){
                 myopb(n, &x[k][0], &w[0][0], alpha);
 
                 for (i=0; i<n; i++)
                    u[i]=e2 * w[0][i] - u[i];

                 myopb(n, u, &w[0][0], alpha);

                 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++) {
   myopb(n, &x[k][0], &w[0][0], alpha);

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



void s_sis1::intros(FILE* fp_out2, long ks, long kg, long kh, double *f, long m)

{
  long l, i, j;

  l=kh+1;

  ksteps=ks;
  maxm=imax(maxm,m);

  fprintf(fp_out2,"\n%8.ld%8.ld%8.ld%8.ld%15.6e\n",m,ks,kg+1,kh+1,f[0]);
 
 if (l>=1)
  for (i=1;i<=l;i++){
   for (j=0;j<32;j++) fprintf(fp_out2," ");
   fprintf(fp_out2,"%15.6e\n",f[i]);
  }
}

/************************************************************** 
 * multiplication of 2-cyclic matrix B by a vector x, where   *
 *							      *
 * B = [D  A]						      *
 *     [A' D] , where A is nrow by ncol (nrow >> ncol). Hence,*
 * B is of order n = nrow+ncol (y stores product vector).     *
 **************************************************************/ 
void s_sis1::myopb(long n,double *x, double *y , double Alpha)

{
   long i, j, end;
   double *tmp;
 
   mxvcount += 2;
   for (i = 0; i < n; i++) y[i] = Alpha*x[i];

   tmp = &x[nrow]; 
   for (i = 0; i < ncol; i++) {
      end = pointr[i+1];
      for (j = pointr[i]; j < end; j++)

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

   return;
}