tms2.c

求矩阵奇异分解svd算法

          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++) {
     opb(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);
}  
/* end of tsvd2 */
#include <stdio.h>
#include <math.h>
#include "tmsc.h"


/*#define ZERO 0.0
#define ONE  1.0*/

float epslon(float x)
{
  /* estimate unit roundoff inquantities of size x */

  float  a,b,c,eps;

  a = 4.0/3.0;
  do {
     b = a - ONE;
     c =b + b + b;
     eps = fabs(c -ONE);
  }
  while (eps == ZERO);  
  eps = eps*fabs(x);
  return(eps);
}
#include <stdio.h>
#include <math.h>
#include <errno.h>
#include <fcntl.h>

extern int ncol,nrow;
extern char *error[];
/***********************************************************************
 *                    check_parameters()                               *
 ***********************************************************************
   Description
   -----------
   Function validates input parameters and returns error code (int)

   Parameters
   ----------
  (input)
   p        desired number of eigenpairs of B.
   maxi     upper limit of desired number of trace steps or iterations.
   n        dimension of the eigenproblem for matrix B
   vec      1 indicates both singular values and singular vectors are wanted.
            0 indicates singular values only.
   nnzero   number of nonzero elements in input matrix (matrix A)

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

int check_parameter(FILE *fp_out1, int maxi, int nmax, int nzmax,
                     int p, int vec, int nnzero, int n)
{
   int error_index, ncells;

   error_index = 0;
   if(ncol >= nmax || nnzero > nzmax) error_index = 1;
   else if(p > maxi)  error_index = 3;
   else if(n <= 0)  error_index = 4;
   else if(maxi <= 0 )  error_index = 5;
   else if(p <= 0 || p > maxi)  error_index = 6;
   if(error_index)
     fprintf(fp_out1, "%s\n", error[error_index]);
   return(error_index);
}
/************************************************************** 
 * Function copies a vector x to a vector y	     	      *
 * Based on Fortran-77 routine from Linpack by J. Dongarra    *
 **************************************************************/ 

void dcopy(int n,float *dx,int incx,float *dy,int incy)

{
   int i;

   if (n <= 0 || incx == 0 || incy == 0) return;
   if (incx == 1 && incy == 1) 
      for (i=0; i < n; i++) *dy++ = *dx++;

   else {
      if (incx < 0) dx += (-n+1) * incx;
      if (incy < 0) dy += (-n+1) * incy;
      for (i=0; i < n; i++) {
         *dy = *dx;
         dx += incx;
         dy += incy;
      }
   }
   return;
}
/************************************************************** 
 * Function scales a vector by a constant.     		      *
 * Based on Fortran-77 routine from Linpack by J. Dongarra    *
 **************************************************************/ 

void dscal(int n,float da,float *dx,int incx)

{
   int i;

   if (n <= 0 || incx == 0) return;
   if (incx < 0) dx += (-n+1) * incx;
   for (i=0; i < n; i++) {
      *dx *= da;
      dx += incx;
   }
   return;
}
/************************************************************** 
 * Function interchanges two vectors		     	      *
 * Based on Fortran-77 routine from Linpack by J. Dongarra    *
 **************************************************************/ 

void dswap(int n,float *dx,int incx,float *dy,int incy)

{
   int i;
   float dtemp;

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

void   dgemv(int, int, int, float, float **,
             float *, float, float *);
float ddot(int, float *,int, float *, int);
void   dscal(int, float, float *,int);
void   daxpy(int, float, float *,int, float *, int);
void   dcopy(int, float *, int, float *, int);

/***********************************************************************
 *                                                                     *
 *                        orthg()                                      *
 *         Gram-Schmidt orthogonalization procedure                    *
 *                                                                     *
 ***********************************************************************/
/***********************************************************************

   Description
   -----------

   The p by n matrix Z stored row-wise in rows f to (f+p-1) of
   array X is reorthogonalized w.r.t. the first f rows of array X.
   The resulting matrix Z is then factored into the product of a
   p by n orthonormal matrix (stored over matrix Z) and a p by p
   upper-triangular matrix (stored in the first p rows and columns 
   of array B).  (Based on orthog from Rutishauser) 


   Parameters
   ----------

   (input)
   p           number of consecutive vectors of array x (stored row-wise)
	       to be orthogonalized
   f           number of rows against which the next p rows are to be
	       orthogonalized
   n           column dimension of x
   x           2-dimensional array whose p rows are to be orthogonalized
	       against its first f rows
   temp        work array


   (output)
   x           output matrix whose f+p rows are orthonormalized
   b           p by p upper-triangular matrix


   Functions called
   --------------

   BLAS         dgemv, ddot, dscal, daxpy, dcopy

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

void orthg(int p, int f, int n, float **b, float **x, float *temp)

{
   int fp, k, km1;
   int orig, small;
   float t, s;

   if (!p) return;
   if (f == 0 && p > n) {
      fprintf(stderr,"%s\n",
         "*** ON ENTRY TO ORTHG, MATRIX TO BE ORTHONORMALIZED IS SINGULAR");
      exit(-1);
   }
   fp = f + p;

   for (k = f; k < fp; k++) {
      km1 = k - 1;
      orig = TRUE;

      while(TRUE) {
         t = ZERO;

	 if (km1 >= 0) {
	    if (km1 > 0) {
	       dgemv(NTRANSP, k, n, ONE, x, x[k], ZERO, temp);
	       t += ddot(k, temp, 1, temp, 1);
	    }

	    else {
	       temp[0] = ddot(n, x[0], 1, x[k], 1);
	       t += temp[0] * temp[0];
	    }

	    if (orig && km1 >= f) 
               dcopy(k - f, &temp[f], 1, &b[k - f][0], 1); 

            if (km1 > 0) 
	       dgemv(TRANSP, k, n, -ONE, x, temp, ONE, &x[k][0]);
            else
	       daxpy(n, -temp[0], x[0], 1, x[k], 1);
         }

	 if (km1 < 0 || p != 1) {
	    s = ddot(n, x[k], 1, x[k], 1);
	    t += s;
	    if (s > t/CONST) {
	       small = FALSE;
	       s = sqrt(s);
               b[k - f][k - f] = s;
	       if (s != ZERO) s = ONE/s;
	       dscal(n, s, x[k], 1);
	    }
	    else {
	       small = TRUE;
	       orig  = FALSE;
	    }
	 }
	 if (small == FALSE || p == 1) break;
      }
   }
}
#include <stdio.h>
#include <math.h>

void dscal(int, float, float *, int);
void opb(int, float *, float *, float);
void daxpy(int, float, float *, int, float *, int);

void pmul(int n, float *y, float *z, float *z2, float *z3,
          float *z4, float e, int degree, float alpha)

/* Compute the multiplication z=P(B)*y (and leave y untouched).
   P(B) is constructed from chebyshev polynomials. The input
   parameter degree specifies the polynomial degree to be used.
   Polynomials are constructed on interval (-e,e) as in Rutishauser's
   ritzit code                                                    */
{

float tmp2,alpha2,eps,tmp;
 
int i,kk;

  alpha2 = alpha*alpha;
  tmp2 = 1.0;
  eps = tmp2 + 1.0e-6;
  if(!degree) for(i=0; i<n; i++) z[i] = y[i];
  else if(degree==1) opb(n,y,z,alpha2);
     else if(degree>1) {
          tmp2 = 2.0/e;
          tmp = 5.1e-1*tmp2;
          for(i=0; i<n; i++)
             z[i] = y[i];
          dscal(n,tmp,z,1);
          opb(n,z,z3,alpha2);
          for(kk=1; kk<degree; kk++) {
             for(i=0; i<n; i++) {
                z4[i] = z3[i];
                z[i] = -z[i];
             }
             opb(n,z3,z2,alpha2);
             daxpy(n,tmp2,z2,1,z,1);
             if(kk != degree-1) {
               for(i=0; i<n; i++) z3[i] = z[i];