s_tms2.cc

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

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

     if((job==2) && (!shift)) { 
       if(ndegre) {
         degree = ndegre;
         *maxd = imax(*maxd,degree);
         polyac = TRUE;
       }
       else polyac = FALSE;
     }
/* section formation */

     t1 = timer();

     for(j=0; j<s; j++) {
        for(i=0; i<n; i++) work2[j][i] = y[j][i];
     }
     if(polyac) {
       porth(s,0,n,work2,work1[0],work1[1],work4[0],work4[1],
             work4[2],e,degree,alpha);
       mxv[0] += degree*s;
       mxv[1] += degree*s;
     }
     else {
       orthg(s,0,n,work1,work2,&work4[0][0]);
     }

   t1 = timer() - t1;
   sec1 = sec1 + t1;
  
     /* form work1(1:n,iptr:s) = b(1:n,1:n)*work2(1:n,iptr:s) */
  
     t1 = timer();
  
  if(polyac) {
    for(j=0; j<left; j++) {
       for(i=0; i<n; i++) work1[iptr+j-1][i] = work2[iptr+j-1][i];
    }
  }
  else {
     for(j=0; j<left; j++)
        myopb(n,work2[iptr+j-1],work1[iptr+j-1],ZERO);
  }

 
/* form work3(1:left,1:left) = work2(1:n,iptr:s)'*work1(1:n,iptr:s) */

  dgemm3(TRANSP, NTRANSP, left, left, n, ONE, work2, 0,iptr-1,
         work1, 0, iptr-1, ZERO, work3, 0, 0);
     
  if(job >= 1) {
    for(j=0; j<left; j++)
       work4[0][j] = dasum(left, &work3[0][j], s) - fabs(work3[j][j]);
  }

     t1 = timer() - t1;
     sec21 = sec21 + t1;

     /*    compute the eigenvalues and eigenvectors of work3:  */
     
     t1 = timer();

     /*    eigenvectors overwrite array work3(:,:)
         store current e-values in work5(:,2)  */
  if(polyac) {
    tred2(0,left,work3,work5[0],work4[1],work3);
    ierr = tql2(0,left,work5[0],work4[1],work3);
    if(ierr) printf("IERR FROM TQL2 = %ld\n",ierr);
  }
  else {
     
     for(j=iptr-1; j<s; j++) work5[1][j] = sig[j];
     
     tred2(0,left, work3, &sig[iptr-1], &work4[1][0], work3);
     ierr = tql2(0,left, &sig[iptr-1], &work4[1][0], work3);
    if(ierr) printf("IERR FROM TQL2 = %ld\n",ierr);
  } 
  t1 = timer() - t1;
  sec22 = sec22 + t1;
     
     /* form y(1:n,iptr:s) = work2(1:n,iptr:s)*work3(1:left,1:left) */

     t1 = timer();

  dgemm3(NTRANSP, NTRANSP, n, left, left, ONE, work2, 0, iptr-1,
           work3, 0, 0, ZERO, y, 0, iptr-1);

     t1 = timer() - t1;
     sec23 = sec23 + t1;
     
     /*    test for convergence here */
     
     t1 = timer();

     for(j=iptr-1; j<s; j++) {
        myopb(n, y[j], work2[j], ZERO);
        if(polyac) {
          work5[1][j]=sig[j];
          sig[j] = ddot(n,y[j],1,work2[j],1)/ddot(n,y[j],1,y[j],1);
        }
        daxpy(n, -sig[j], y[j], 1, work2[j],1);
        work4[2][j-iptr+1] = sqrt(ddot(n,work2[j],1,work2[j],1));
        if(j < p) titer[j] += 1;
     }
     mxv[0] += s-iptr+1;
     mxv[1] += s-iptr+1;
     t1 = timer() - t1;
     sec3 = sec3 + t1;
    temp_left=left; 
    /* array work4(:,3) stores the vector 2-norms of residual vectors */
     for(k=0; k<temp_left; k++) {
         if(fabs(work4[2][k]) <= tol) {
            iptr = iptr + 1;
            if(iptr > p)  {
                          goto L10;
                          }
            left = s-iptr+1;
         }
     }

     if((!shift) && (job>=1)) {
     
       /* work4(:,1) stores gershgorin radii  *
        * iwork(:,1) is the clustering vector */ 

       disk(left, &sig[iptr-1], work4[0], iwork[0], iwork[1]);
       
       for(k=iptr-1; k<s; k++) {
          if(iwork[0][k-iptr+1] == 1)
            isol(k, work4[2][k-iptr+1], red, lwork, work5[2],
                 work5[3],s);
          else
             if(iwork[0][k-iptr+1] > 1)
               clus(k,iwork[0][k-iptr+1],&work4[2][k-iptr+1],red,lwork,
                    work5[2], work5[3], work4[1],s);
       }
     }             
       
     /* continue algorithm...  */
       
     /* shift start */
     
     if((iter>0) && (!shift) && (job>=1)) {
        if(lwork[iptr-1]) {
           shift = TRUE;
           polyac = FALSE;
           orthg(s, 0, n, work1, y, &work4[0][0]);
        }
     }

     if(shift) {
     
     /* compute shifts in work5(:,1) here: */
  

      for(k=iptr-1; k<s; k++) {
         if((k>iptr-1) && (k<=p-1)) {
           work5[0][k] = ZERO;
           for(j=iptr-2; j<k; j++) {
              if((j!=-1) && (sig[j] <= (sig[k] - work4[2][k-iptr+1]))) 
                work5[0][k] = sig[j];
              }
 
           if((work5[0][k] == sig[k-1]) && (work5[0][k-1]==sig[k-1]))
             work5[0][k] = sig[k];
           if(work5[0][k] == ZERO) 
             work5[0][k] = work5[0][k-1];
         }
         else if(k>p-1) work5[0][k] = work5[0][p-1];
              else if(k==iptr-1) work5[0][k] = sig[k];
      }
     }

     t1 = timer();

/* monitor polynomial degree (if job=2) */
   if(job==2) {
     enew = sig[s-1];
     e = fabs(enew-alpha*alpha);
     if((sig[iptr-1]>pparm1*sig[s-1]) && (enew<alpha*alpha) &&
        (!shift)) {
       tmp = acosh(sig[s-1]/sig[iptr-1]);
       itmp = 2*imax(((int) (pparm2/tmp)),1);
       if(degree) ndegre = imin(2*degree,itmp);
       else ndegre = 2;
     }
   }
     
  if(polyac) {
    for(j=iptr-1; j<s; j++) {
       pmul(n,y[j],work2[j],work4[0],work4[1],work4[2],
            e,degree,alpha);
       }
    orthg(left,0,n,work1,&work2[iptr-1],work4[0]);

    for(j=iptr-1; j<s; j++) {
       for(i=0; i<n; i++)
          y[j][i]=work2[j][i];
    }

    dtrsm('R','U',TRANSP,'N',n,left,ONE,work1,&y[iptr-1]);

    mxv[0] += degree*left;
    mxv[1] += degree*left;
 }
 else {
  
     /* do cg iterations */

     /*load y into work1 array for orthog. projector in subroutine cgt*/

     for(j=0; j<left; j++) {
        for(i=0; i<n; i++)
            work1[j][i] = y[iptr+j-1][i];
     }
           
    /* cg loop for independent systems (no shifting used)  */

     if(!shift) 
       for(i=0; i<left; i++)
             cgt(n, left, y[iptr+i-1], work1, &itcgt[iptr+i-1],
                 sig[iptr+i-1], sig[s-1], &iwork[0][i], meps);
     else {
     
         /*shift with work5(:,1) in cg iterations */
       for(i=0; i<left; i++) {
                cgts(n, left, y[iptr+i-1], work1, &itcgt[iptr+i-1],
                     sig[iptr+i-1], work5[1][iptr+i-1], sig[s-1],
                     work5[0][iptr+i-1], &iwork[0][i], meps);
                }
     }

     for(i=0; i<left; i++) {
        mxv[0] += iwork[0][i];
        mxv[1] += iwork[0][i];
     }
  }
     t1 = timer() - t1;
     sec4 += t1;
}
L10:     
   /* compute final 2-norms of residual vectors w.r.t. B */

   for(j=0; j<p; j++) {
     myopb(n, y[j], work2[j], alpha*alpha);
     tmp = ddot(n,y[j],1,work2[j],1)/ddot(n,y[j],1,y[j],1);
     daxpy(n,-tmp,y[j],1,work2[j],1);
     tmp = sqrt(fabs(tmp));
     work4[2][j] = sqrt(ddot(n,work2[j],1,work2[j],1));
     opa(y[j],&y[j][n]);
     tmpi = 1.0/fabs(tmp);
     dscal(nrow,tmpi,&y[j][n],1);
     work4[2][j] = work4[2][j]*tmpi;
     sig[j]=tmp;
   }
   mxv[0] += 2*p;
   mxv[1] += p;
   
   sec2   = sec21 + sec22 + sec23;
   total += sec1 + sec2 + sec3 + sec4;
   
   /* load output time and mxv arrays */
   time[0]=sec1;
   time[1]=sec2;
   time[2]=sec3;
   time[3]=sec4;
   time[4]=total;
   mxv[2]=mxv[0]+mxv[1];
  
  /* load residual vector */
  for(i=0; i<p; i++)
      res[i] = work4[2][i];
  return(ierr);
}