tms1.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 "tmsg.h"
#include "tmsc.h"

/* #define  UNIX_CREAT */

#ifdef UNIX_CREAT
#define PERMS 0664
#endif

long tsvd1(FILE *, long, long, long, long, double, double,
           double *, long, long *, long *, long *, double *, double *,
           long *, double **, double **, double **, double **,
           double **, double **, long **, long *);
long   check_parameter(long, FILE *, long, long, long, long, long, long);
double norm_1();
float  timer(void);


/***********************************************************************
 *                                                                     *
 *                        main()                                       *
 *     Sparse SVD Via Eigensystem of Equivalent 2-Cyclic Matrix        *
 *                  (double precision)                                 *
 *                                                                     *
 ***********************************************************************

   Description
   -----------
   This sample program tms1 uses  subroutine tsvd1 to compute the p-largest
   singular triplets of A via the equivalent symmetric eigensystem of

       [alpha*I   A]
   B = [           ]  , where A is m (nrow) by n (ncol) (nrow >> ncol),
       [A'  alpha*I]

   so that {u, abs(lambda-alpha), v} is a singular triplet of A.
   alpha is chosen so that B is (symmetric) positive definite. In
   the calling program,alpha is determined as the matrix 1-norm of
   A. The user can set alpha to be any known upper bound for the 
   largest singular value of the matrix A.
   (A' = transpose of A)

   User supplied routines:  opb,opat,timer

   opat(x,y) takes an m-vector x and should return A'*x in y.
   opb(n,x,y,shift) takes an n-vector x and should return D*x in y,
                where D=[B-shift*I],I is the n-th order identity matrix.

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

   tms1 utilizes ritz-shifting for acceleration of convergence.

   External parameters
   -------------------
   Defined and documented in tmsg.h

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

   p     : number of singular triplets (largest) desired
   s     : initial subspace dimension
   job   : controls use of ritz shifting (0=no , 1=yes)
   tol   : user-specified tolerance for residuals
   red   : residual norm reduction factor for initiating shifting
   v     : output singular vectors ? (boolean)
   n     : dimension of the eigenproblem for matrix B(nrow + ncol)

  (output)

  sig    : array of singular approximate values
  time   : time breakdown
  iwall  : wall clock time   

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

   BLAS         daxpy,dscal,ddot,dcopy,dswap,dgemm,dgemv,dasum
   USER         opb,opat,timer
   MISC         write_header,check_parameters
   TMS1         tsvd1

   Precision
   ---------

   All floating-point calculations use double precision,
   variables are declared as long or double.

   TMS1 development
   ----------------

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

   August  24 1992:  Date written
   October 23 1996:  Updated main() to produce correct s-vector output
                     to fp_out2

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

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

main()
{
   double  time[5], red, exetime, *sig, *res;
   double  *tempptr1, *workptr1, total, **tempptr2, **workptr2, meps;

   long    *workptr3, **workptr4, *tempptr3, **tempptr4,length,memsize;
   long    k, i, j, id, p, size1, size2, s, job, mem, vec, nnzero;
   long    n,iwall, mxv[3], itime[4], *itcgt, *titer, maxi;

   char    title[73], name[41], v[6], *in1, *in2, *out1, *out2;
   FILE    *fp_in1, *fp_in2;

   in1 = "tmp1";
   in2 = "matrix";
   out1 = "tmo1";
   out2 = "tmv1";

   /* 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);
   }

   /* read data */
   fscanf (fp_in2,"%72c%*s%*s%*s%ld%ld%ld%*d",
           title, &nrow, &ncol, &nnzero);
   title[73] = '\0';
   fscanf (fp_in1,"%s %ld %ld %ld %lf %lf %s %ld", 
           name, &p, &s, &job, &tol,&red,v,&maxi);
   if(!(strcmp(v, "TRUE"))) {
     vec = 1;
#if  !defined UNIX_CREAT
       if ((fp_out2 = open(out2, O_CREAT | O_RDWR)) == -1) {
          printf("cannot open output file %s \n", out2);
          exit(-1);
       }
#else
       if ((fp_out2 = creat(out2, PERMS )) == -1) {
          printf("cannot open output file %s \n", out2);
          exit(-1);
       }
#endif
   }
   else vec = 0;

   n = nrow + ncol;

   /* even though the validity of the parameters will be checked in the
    * SVD code, some parameter checking should also be done before
    * allocating memory to ensure that they are non-negative */

   if(check_parameter(n, fp_out1, maxi, NMAX, NZMAX, p, vec, nnzero)) {
     fclose(fp_in1);
     fclose(fp_in2);
     fclose(fp_out1);
     if(vec) close(fp_out2);
     exit(-1);
   }
   /*******************************************************************
    * 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)  *
    *******************************************************************/
        size1 = sizeof(double) * (nnzero);
        size2 = sizeof(long) * (ncol + nnzero + 1);
        if(!(pointr = (long *)   malloc(size2))   ||
           !(value  = (double *) malloc(size1))) {
            perror("FIRST MALLOC FAILED IN TMS1()");
            exit(errno);
        }

   /* calculate memory allocated */
   mem = size1 + size2;

   rowind = pointr + (ncol + 1);

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

   /* 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]);

  /* allocate memory for arrays */

   memsize = sizeof(double) * ((n * (s + s + 3 + s)) + ( s*(4 + s) ) 
                                + s + s  );

   mem   += memsize;

   /*******************************************************************
    *       Allocate work area and initialize pointers                *
    *                                                                 *
    *       pointer                size                               *
    *                                                                 *
    *    w1   (work1)           n  ( n x  s)                          *
    *    w2   (work2)           n  ( n x  s)                          *
    *    w3   (work3)           s  ( s x  s)                          *
    *    w4   (work4)           n  ( n x  3)                          *
    *    w5   (work5)           s  ( s x  4)                          *
    *    yy   (y     )          n  ( n x  s)                          *
    *    sig                    s                                     *
    *    res                    s                                     *
    *                                                                 *
    *    (memory allocated separately)                                *
    *                                                                 *
    *    iwork(iw   )           s  ( s x  2)                          *
    *    lwork                  s                                     *
    *    titer                  s                                     *
    *    itcgt                  s                                     *
    *                                                                 *
    *******************************************************************/
   if(!(workptr1 = (double *)malloc(memsize))) {
      perror("SECOND MALLOC FAILED IN TMS1()");
     exit(errno);
   }
  
   memsize = sizeof(double *) * (s + s + s + 3 + 4 + s);
   mem   += memsize;

   if(!(workptr2 = (double **)malloc(memsize))) {
     perror("THIRD MALLOC FAILED IN TMS1()");
     exit(errno);
   }
   tempptr1 = workptr1;
   tempptr2 = workptr2;
 
   length    = n * s;
   w1        = tempptr1;
   tempptr1 += length;
   work1     = tempptr2;
   tempptr2 += s;
   j=0;
   for(i=0; i<length; i+=n) work1[j++] = &w1[i];
 
   length    = n * s;
   w2        = tempptr1;
   tempptr1 += length;
   work2     = tempptr2;
   tempptr2 += s;
   j=0;
   for(i=0; i<length; i+=n) work2[j++] = &w2[i];
 
   length    = s * s;
   w3        = tempptr1;
   tempptr1 += length;
   work3     = tempptr2;
   tempptr2 += s;
   j=0;
   for(i=0; i<length; i+=s) work3[j++] = &w3[i];
 
   length    = n * 3;
   w4        = tempptr1;
   tempptr1 += length;
   work4     = tempptr2;
   tempptr2 += 3;
   j=0;
   for(i=0; i<length; i+=n) work4[j++] = &w4[i];

   length    = s * 4;
   w5        = tempptr1;
   tempptr1 += length;
   work5     = tempptr2;
   tempptr2 += 4;
   j=0;
   for(i=0; i<length; i+=s) work5[j++] = &w5[i];
 
   length    = n * s;
   yy        = tempptr1;
   tempptr1 += length;
   y         = tempptr2;
   tempptr2 += s;
   j=0;
   for(i=0; i<length; i+=n) y[j++] = &yy[i];
 
   sig       = tempptr1;
   tempptr1 += s;
 
   res       = tempptr1;
   tempptr1 += s;
 
   /* Allocate memory for logical array lwork (long 0= false and 1=true)
      and integer array iwork                                    */
 
   memsize = sizeof(long) * ( s * (2 + 1) + s + s );
   mem   += memsize;
 
   if(!(workptr3 = (long *)malloc(memsize)) ) {
      perror("FOURTH MALLOC FAILED IN TMS1()");
      exit(errno);
   }
 
   memsize = sizeof(long *) * 2;
   mem   += memsize;
 
   if(!(workptr4 = (long **)malloc(memsize))) {
      perror("FIFTH MALLOC FAILED IN TMS1()");
      exit(errno);
   }
 
   tempptr3 = workptr3;
   tempptr4 = workptr4;
 
   length    = s * 2;
   iw        = tempptr3;
   tempptr3 += length;
   iwork     = tempptr4;
   tempptr4 += 2;
   j=0;
   for(i=0; i<length; i+=s) iwork[j++] = &iw[i];
 
   lwork     = tempptr3;
   tempptr3 += s;
 
   titer     = tempptr3;
   tempptr3 += s;
 
   itcgt     = tempptr3;
   tempptr3 += s;
 

 
   /* estimate alpha (via matrix 1-norm)
    * used for upper bound on max singular value of matrix A */

   alpha = norm_1();
   exetime = timer();

   /* make a trace min. run; exit upon error */ 
   ierr = tsvd1(fp_out1, n, p, s, job, tol, red, sig, maxi,
                &mem, itcgt, titer, time, res,  mxv, work1,
                work2, work3, work4, work5, y, iwork, lwork);
   if(ierr) {
     free(value);
     free(pointr);
     fclose(fp_in1);
     fclose(fp_in2);
     fclose(fp_out1);
     if(vec) close(fp_out2);
     free(workptr1);
     free(workptr2);
     free(workptr3);
     free(workptr4);
     exit(-1);
   } 
   exetime  = timer() - exetime;
   iwall    = (int) (fmod(exetime+1000000.0,1000000.0));
   itime[0] = (int) (100*(time[0]/time[4]));
   itime[1] = (int) (100*(time[1]/time[4]));
   itime[2] = (int) (100*(time[2]/time[4]));
   itime[3] = (int) (100*(time[3]/time[4])); 

   /* write output */

   fprintf(fp_out1, " ... \n");
   fprintf(fp_out1, " ... TRACE MINIMIZATION ON THE  \n");
   fprintf(fp_out1, " ... EQUIVALENT CYCLIC   E_PROBLEM\n");
   fprintf(fp_out1, " ... NO. OF EQUATIONS            =%10ld\n", n);
   fprintf(fp_out1, " ... MAX. NO. OF ITERATIONS      =%10ld\n", maxi);
   fprintf(fp_out1, " ... NO. OF ITERATIONS TAKEN     =%10ld\n",
           titer[p-1]);
   fprintf(fp_out1, " ... NO. OF DESIRED EIGENPAIRS   =%10ld\n",p);
   fprintf(fp_out1, " ... INITIAL SUBSPACE DIM.       =%10ld\n",s);
   fprintf(fp_out1, " ... FINAL   SUBSPACE DIM.       =%10ld\n",s-p);
   fprintf(fp_out1, " ... MEMORY REQUIRED (BYTES)     =%10ld\n",mem);
   fprintf(fp_out1, " ... JOB: 0=NO SHIFT, 1=SHIFT    =%10ld\n",job);
   if (vec)
     fprintf(fp_out1, " ... WANT S-VECTORS?   [T/F]   =           T\n");
   else
     fprintf(fp_out1, " ... WANT S-VECTORS?   [T/F]   =           F\n");
   fprintf(fp_out1, " ... ALPHA                       =%10.2E\n",alpha);
   fprintf(fp_out1, " ... FINAL RESIDUAL TOLERANCE    =%10.2E\n" ,tol);
   fprintf(fp_out1, " ... RESIDUAL REDUCTION TOL.     =%10.2E\n",red);
   fprintf(fp_out1, " ... MULTIPLICATIONS BY A        =%10ld\n",mxv[0]);
   fprintf(fp_out1, " ... MULTIPLICATIONS BY A^T      =%10ld\n",mxv[1]);
   fprintf(fp_out1, " ... TOTAL NUMBER OF MULT.       =%10ld\n",mxv[2]);
   fprintf(fp_out1, " ... IERR FROM SUBR. TSVD1       =%10ld\n\n",ierr);
   fprintf(fp_out1,"... CPU TIME BREAKDOWN:\n");
   fprintf(fp_out1,"... GRAM-SCHMIDT ORTHOG.          =%10.2E (%3ld%%) \n",
           time[0],itime[0]);
   fprintf(fp_out1,"... SPECTRAL DECOMPOSITION        =%10.2E (%3ld%%) \n",
           time[1],itime[1]);
   fprintf(fp_out1,"... CONVERGENCE CRITERIA          =%10.2E (%3ld%%) \n",
           time[2],itime[2]);
   fprintf(fp_out1,"... CONJUGATE GRADIENT            =%10.2E (%3ld%%) \n",
           time[3], itime[3]);
   fprintf(fp_out1,"... TOTAL CPU  TIME (SEC)         =%10.2E\n", time[4]);
   fprintf(fp_out1,"... WALL-CLOCK TIME (SEC)         =%10ld\n\n", iwall);
   fprintf(fp_out1, " %s\n", title);
   fprintf(fp_out1, "           %s\n", name);
   fprintf(fp_out1, " ... NO. OF TERMS     (ROWS)   = %10ld\n", nrow);
   fprintf(fp_out1, " ... NO. OF DOCUMENTS (COLS)   = %10ld\n", ncol);
   fprintf(fp_out1, " ... ORDER OF MATRIX A         = %10ld\n", n);
   fprintf(fp_out1,"......\n");
   fprintf(fp_out1, "...... COMPUTED SINGULAR VALUES  (RESIDUAL NORMS)  T-MIN STEPS CG STEPS\n"); 
   fprintf(fp_out1,"......\n");
   for(i=0; i<p; i++)
   fprintf(fp_out1, "... %2ld  %24.14E ( %10.2E)       %ld           %ld \n",
           i+1, sig[i], res[i],titer[i],itcgt[i]); 
   if (vec) {
       for (i = 0;  i < p;  i++)
          write(fp_out2, (void *) &y[i][0], sizeof(double) * n); }

   free(value);
   free(pointr);
   fclose(fp_in1);
   fclose(fp_in2);
   fclose(fp_out1);
   if (vec) close(fp_out2);
   free(workptr1);
   free(workptr2);
   free(workptr3);
   free(workptr4);
   exit(0);