sis2.c

求矩阵奇异分解svd算法

/*************************************************************************
                           (c) Copyright 1993
                        University of Tennessee
                          All Rights Reserved                          
 *************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <errno.h>
#include <math.h>
#include <fcntl.h>
#include "sisg.h"
#include "sisc.h"

/* #define  UNIX_CREAT */

#ifdef UNIX_CREAT
#define PERMS 0664
#endif

void   dscal(long, double, double *, long);
void   daxpy(long, double, double *,long, double *, long);
double ddot(long, double *,long, double *, long);
void   opa(long, double *, double *);
void   opb(long, double *, double * );
void   intros(long , long , long , double * , long);
void   write_data(long , long , long , long , long ,  long , long , double , char *, char *);
float  timer(void);
void   ritzit( long , long , long , double , void (*) (long , double *, double *), void (*)(long , long , long , double *, long ),long , double **, double *, double *, double *, double *, long *);
double fsign(double, double);
double dmax(double, double);
double dmin(double, double);
long   imin(long , long );
long   imax(long, long);
#define ABS(x) ((x) > (ZERO) ? (x) : (-(x)))



/***************************************************************
 *                                                             *
 *                    main()                                   *
 *    Sparse SVD via Eigensystem of them matrix A'A            *
 *              (double precision)                             *
 *                                                             *
 ***************************************************************/
/***************************************************************

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

 this sample program uses ritzit to compute singular triplets of
 the matrix A via the equivalent symmetric eigenvalue problem

 B=A'A, A is m (nrow) by n (ncol)    (nrow >> ncol), and X'=[u', v'],

 so that {u, sqrt(lambda, v} is a singular triplet of A.
 (A' = transpose of A)

 user supplied routines: opb, opa

 opa( n, x, y) takes an n-vector x and should return A*x
                      in y
 opb( n, x, y) takes an n-vector x and should return B*x
                      in y.

 User should edit timer() with an appropriate call to an
 intrinsic timing routine that returns elapsed user cpu time.


 External parameters (constants)
 -------------------------------
 Defined and documented in sisg.h (sisc.h)
 All constants are in upper-case.

 Local parameters
 ----------------
 (input)

 em2		no. of  eigenpairs of B approximated
 numextra	number of extra vectors to carry
 km		max. number of iterations
 p		initial subspace dimension (= em2 + numextra)
 em2		no. of desired eigenpairs
 eps		tolerance for eigenpair residuals
 x		two dimensional array of iteration vectors
 cx,u   	work arrays

 (output)

 d		array of eigenvalues of B (squares of the singular
                values of A)
 f		array of residual norms
 u		left-singular vectors
 x		right-singular vectors
 
 Functions used
 --------------
 BLAS         daxpy, ddot
 USER         opb, timer, intros
 MISC         write_data 
 SIS2         ritzit

 Precision
 ---------
 All floating-point calculations use double precision;
 Variables are declared as long and double

 SIS2 development
 ----------------
 SIS2 is a C translation of the Fortran-77 SIS2 from the SVDPACK
 library written by Michael W. Berry, University of Tennessee,
 Dept of Computer Science, 107 Ayres hall, Knoxville, 
 TN 37996-1301

 Date Written: 21 Jun 1992

 Sowmini Varadhan
 University of Tennessee
 Dept. of Computer Science
 107 Ayres Hall
 Knoxville, TN 37996-1301
 
 internet: [e-mail]

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


main() 

{
   float t0, exetime;
   long k, j, i, ii, n,  nnzero, size1, size2,  vectors;
   double *x[NSIG], *cx, *f, *d, *u  , *tptr ;
   long em, numextra, km, p, em2, imem;
   double eps, dsum, tmp1, tmp2, tmp3, tmp4, xnorm;
   char title[73], name[41], v[6];
   char *in1, *in2, *out1, *out2, *out3;
   FILE *fp_in1, *fp_in2 ;
   long fp_out3;
 
   in1 = "sip2"     /* input parameters    */;
   in2 = "matrix"   /* matrix              */;
   out1 = "sio2"    /* results             */;
   out2 = "sio5"    /* info                */;
   out3 = "siv2"    /* singular vectors    */;

   /* open files for input/output */
    if (!(fp_in1 = fopen(in1, "r"))) { 
       printf("cannot open file %s for reading\n", in1);
       exit(-1);
    }
    if (!(fp_in2 = fopen(in2, "r"))) {
       printf("cannot open file %s for reading\n", in2);
       exit(-1);
    }
    if (!(fp_out1 = fopen(out1, "w"))) { 
       printf("cannot open output file %s \n", out1);
       exit(-1);
    }
    if (!(fp_out2 = fopen(out2, "w"))) { 
       printf("cannot open output file %s \n", out2);
       exit(-1);
    }
	
 
   /*write header of output file*/   

   fprintf(fp_out2,"\n\n ----------------------------------------\n");
   fprintf(fp_out2," INTERMEDIATE OUTPUT PARMS:\n\n");
   fprintf(fp_out2," M:=CURRENT DEGREE OF CHEBYSHEV POLYNOMIAL\n");
   fprintf(fp_out2," S:=NEXT ITERATION STEP\n");
   fprintf(fp_out2," G:=NUMBER OF ACCEPTED EIGENVALUES OF B\n");
   fprintf(fp_out2," H:=NUMBER OF ACCEPTED EIGEN VECTORS OF B\n");
   fprintf(fp_out2," F:=VECTOR OF ERROR FOR EIGENPAIRS OF B\n");
   fprintf(fp_out2," ---------------------------------------------\n");
   fprintf(fp_out2," M         S       G       H         F     \n");
   fprintf(fp_out2," ---------------------------------------------");

   /* read data */
    fscanf (fp_in2,"%72c%*s%*s%*s%ld%ld%ld%*d",
	    title, &nrow, &ncol, &nnzero);
    title[73] = '\0';

    /* skip data format line */
    fscanf(fp_in2,"%*s %*s %*s %*s");

    n = ncol; 
 
    if (nnzero>NZMAX)
    {  fprintf(stderr,"sorry, your matrix is too big\n");
       exit(-1);
    }

    fscanf(fp_in1,"%s%ld%ld%ld%lf%s",name, &em, &numextra, &km, &eps, v);
 
    if (!(strcmp(v, "TRUE"))) {
        vectors = 1;

#if  !defined UNIX_CREAT
       if ((fp_out3 = open(out3, O_CREAT | O_RDWR)) == -1) {
          printf("cannot open output file %s \n", out3);
          exit(-1);
       }
#else
       if ((fp_out3 = creat(out3, PERMS )) == -1) {
          printf("cannot open output file %s \n", out3);
          exit(-1);
       }
#endif


     }
    else vectors = 0;

    p= em + numextra;
    em2=em;
    /*******************************************************************
     * allocate memory						       *
     * pointr - column start array of harwell-boeing sparse matrix     *
     *          format                                       (ncol+1)  *
     * rowind - row indices array of harwell-boeing format   (nnzero)  *
     * value  - nonzero values array of harwell-boeing sparse matrix   *
     *          format                                       (nnzero)  *
     * x      - 2 dimensional array of  iteration vectors    (NSIG*n)  *
     * d      - array of eigenvalues of B                    (p)       *
                         (squares of the singular values of A)         *
     * f      - temporary storage array                      (p)       *
     * cx     - temporary storage array (control quantities) (n)       *
     * u      - work array                                   (nrow)    *
     *******************************************************************/
	 size1 = sizeof(double) * (nnzero + p* 2 + n  + nrow + NSIG * n);
	 size2 = sizeof(long) * (ncol + nnzero + 1);
	 if (!(pointr = (long *)   malloc(size2))   ||
	     !(value  = (double *) malloc(size1))){
	     perror("MALLOC FAILED in MAIN()");
	     exit(errno);
         }

         imem=size1 - (nnzero * sizeof(double));

         rowind=pointr + (ncol+1); 

         tptr=value + nnzero;
         for (ii=0;ii<NSIG; ii++){
          x[ii]=tptr;
	  tptr+=n;
         }
         d=tptr;
         tptr+=p;
         f=tptr;
         tptr+=p;
         cx=tptr;
         tptr+=n;
         u  =tptr;
	 
    /* read data */
    for (i = 0; i <= ncol; i++) fscanf(fp_in2, "%ld", &pointr[i]);
    
    for (i = 0; i < ncol; i++) pointr[i] -= 1;
 
    /* define last element of pointr in case it is not */ 
    pointr[i] = nnzero;
    for (i = 0; i < nnzero; i++) fscanf(fp_in2, "%ld", &rowind[i]);
    for (i = 0; i < nnzero; i++) rowind[i] -= 1;
    for (i = 0; i < nnzero; i++) fscanf(fp_in2, "%lf", &value[i]);


    exetime = timer();
    
    /*  call ritzit */
    ritzit(n,p,km,eps,opb,intros, em,x, d, f, cx, u, &imem);

    exetime = timer() - exetime;
    write_data(n,km,em2,em,p,imem,vectors,eps,title,name);
    t0=timer();

    for (j=0;j<em2;j++) {
      opb(n, &x[j][0], cx);
      tmp1=ddot(n,&(x[j][0]), 1, cx, 1);
      daxpy(n, -tmp1, &(x[j][0]), 1, cx, 1);
      tmp1=sqrt(tmp1);
      xnorm=sqrt(ddot(n, cx, 1, cx, 1));
      /* multiply by matrix A to get (scaled) left S-vector */
      opa(n, &x[j][0], u);
      tmp2=ONE/tmp1;
      dscal(nrow, tmp2, u, 1);
      xnorm=xnorm*tmp2;
      f[j]=xnorm;
      d[j]=tmp1;
      if (vectors) {
           write(fp_out3, (char *)u, nrow * sizeof(double));
           write(fp_out3, (char *)&x[j][0],n * sizeof(double));
     }
    }

   exetime=exetime + timer() -t0;

   fprintf(fp_out1,"\n ...... SISVD EXECUTION TIME= %10.2e\n",exetime);
   fprintf(fp_out1," ......\n ......");
   fprintf(fp_out1," MXV   = %12.ld\n ......    COMPUTED SINGULAR VALUES  (RESIDUAL NORMS)\n ......", mxvcount);

   for (ii=0;ii<em2;ii++)
      fprintf(fp_out1,"\n ...... %3.ld   %22.14e  (%11.2e)",
                                   ii+1,     d[ii],    f[ii]);

     
   fprintf(fp_out1,"\n");
   free(value);
   free(pointr);
   fclose(fp_in1);
   fclose(fp_in2);
   if (vectors) close(fp_out3);
   fclose(fp_out2);
   fclose(fp_out1);

}

/*********************************************************************/
void write_data(long n, long km, long em2, long em, long p,  long imem, long vectors, double eps, char *title, char *name)
{
 char svectors;
if (vectors) svectors='T'; else svectors='F';
  fprintf( fp_out1,"\n ...\n ... SOLVE THE [A^TA] EIGENPROBLEM\n");
  fprintf(fp_out1," ... NO. OF EQUATIONS           = %10.ld\n",n);
  fprintf( fp_out1," ... MAX. NO. OF ITERATIONS     = %10.ld\n ",km);
  fprintf(fp_out1,"... NO. OF DESIRED EIGENPAIRS  = %10.ld\n ",em2);
  fprintf(fp_out1,"... NO. OF APPROX. EIGENPAIRS  = %10.ld\n ",em);
  fprintf(fp_out1,"... INITIAL SUBSPACE DIM.      = %10.ld\n ",p);
  fprintf(fp_out1,"... FINAL   SUBSPACE DIM.      = %10.ld\n ",p-em);
  fprintf(fp_out1,"... MEM REQUIRED (BYTES)       = %10.ld\n",imem);
  fprintf(fp_out1," ... WANT S-VECTORS? [T/F]      = %10.c\n",svectors);
  fprintf(fp_out1," ... TOLERANCE                  = %10.2e\n ",eps);
  fprintf(fp_out1,"... NO. OF ITERATIONS TAKEN    = %10.ld\n ", ksteps);
  fprintf(fp_out1,"... MAX. DEG. CHEBYSHEV POLY.  = %10.ld\n ...\n ", maxm);
  fprintf(fp_out1,"%s\n %s\n ... NO. OF TERMS     (ROWS)    = %10.d\n",title, name, nrow);
  fprintf(fp_out1," ... NO. OF DOCUMENTS (COLS)    = %10.ld\n ",ncol);
  fprintf(fp_out1,"... ORDER OF MATRIX B          = %10.ld\n ...\n", n)
;

fflush(fp_out1);
}
/***********************************************************************/

void intros(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]);
  }
}
  
#include <math.h>
#include <errno.h>
#include <stdlib.h>
#include <stdio.h>

extern long mxvcount;
double ddot(long , double *, long, double *, long);
void dscal(long, double, double *, long);
void tred2(long, long, double **, double *, double *, double **);
long tql2(long, long, double *, double *, double **);
void daxpy(long, double, double *, long, double *, long);
double dmax(double, double);
long imax(long, long);

#include "sisc.h"

/****************************************************************
 *                                                              *
 *                        ritzit()                              *
 *                                                              *
 ****************************************************************/
/****************************************************************

 Description:
 ------------

 This subroutine is a translation of the Fortran-77 procedure
 RITZIT from the SVDPACK library written by Michael W. Berry,
 University of Tennessee, Dept. of Computer Science, 107 Ayres
 Hall, Knoxville TN 37919-1301

 This subroutine determines the absolutely largest eigenvalues
 and corresponding eigenvectors of a real symmetric matrix by
 simultaneous iteration.

 External parameters (constants)
 -------------------------------

 Defined and documented in sisg.h (sis1c.h)

 local parameters: 
 -----------------
 
 (input)

 n      the order of the matrix  whose eigenpairs are sought
        (matrix B for the SVD problem)
 kp     the number of simultaneous iteration vectors
 km     the maximum number of iteration steps to be performed. If
        starting values for the iteration vectors are available,
        km should be prefixed with a minus sign.
 eps    the tolerance for accepting eigenvectors
 opb    the name of the subroutine that defines the matrix B. opb
        is called with the parameters (n, u, w) and must 
        compute w=Bu without altering u for the SVD problem
 inf    the name of the subroutine that may be used to obtain 
        information or exert control during execution. inf is
        called with parameters (ks, kg, kh, f, m), where ks
        is the number of the next iteration step, kg is the
        number of already accepted eigenvectors, kh is the number
        already accepted eigenvalues, and f is the array of error
        quantities for the vectors of x. An element of f has the
        value 4.0 until the corresponding eigenvalue of the matrix
        B  has been accepted.
 kem    the number of eigenvalues and corresponding eigenvectors
         of matrix B desired. kem must be less than kp.
 x      contains, if km is negative, the starting values for the
        iteration vectors.

 (output)
 
 km     is unchanged
 kem    is reset to the number of eigenvalues and eigenvectors
        of matrix B actually accepted within the limit of km steps
 imem   the number of bytes needed for this invocation
 x      contains in its first kem columns orthonormal
        eigenvectors of the matrix B,  corresponding to the 
        eigenvalues in array d.  The remaining columns contain 
        approximations to further eigenvectors.
 d      contains in its first kem positions the absolutely
        largest eigenvalues of the matrix B. The remaining positions
        contain approximations to smaller eigenvalues.

 u, w, b, f, cx, x, s and  rq are temporary storage arrays.

 Functions used
 --------------