tms2.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   tsvd2(FILE *, long, long, long, long, double, double,
             double *, long, long *, long *, long *, double *, double *,
             long *, double **, double **, double **, double **,
             double **, double **, long **, long *, long *); 
long   check_parameter(FILE *, long, long, long, long, long, long, long);
double norm_1();
long   imin(long, long);
float  timer(void);

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

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

   This sample program tms2 uses subroutine tsvd2 to compute the p-largest
   singular triplets of A via the equivalent symmetric eigensystem of

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

   so that sqrt(alpha*alpha-lambda) is a singular value 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.

   tms2 utilizes ritz-shifting and chebyshev polynomials 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 and poly acceleration
           (0=no acceleration, 1=ritz-shifting only, 2=both acceleration 
                                                       techniques)
   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 approximate singular values
  time   : cpu 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
   TMS2         tsvd2

   Precision
   ---------

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

   TMS2 development
   ----------------

   TMS2 is a C translation of the Fortran-77 TMS2 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,*workptr5,**workptr6;
   double  **tempptr2,**workptr2,*tempptr5,**tempptr6;

   long  *workptr3,**workptr4,*tempptr3,**tempptr4,memsize,length;

   long n,k,i,j, id, p, size1, size2, s, job,mem, vec,nnzero;
   long iwall,mxv[3], itime[4],*itcgt,*titer,maxd, maxi;

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

   in1 = "tmp2";
   in2 = "matrix";
   out1 = "tmo2";
   out2 = "tmv2";

   /* 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 = imin(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 nonnegative */

   if(check_parameter(fp_out1, maxi, NMAX, NZMAX, p, vec, nnzero, n)) {
     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 TMS2()");
            exit(errno);
        }

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

   size1 = sizeof(double) * s;
   size2 = sizeof(long) * s;
   mem = 2 * (size1 + size2);

   if(!(sig = (double *) malloc(size1)) || 
      !(res = (double *) malloc(size1)) ||
      !(itcgt = (long *) malloc(size2)) || 
      !(titer = (long *) malloc(size2))) {
            perror("SECOND MALLOC FAILED IN TMS2()");
            exit(errno);
        }
   /* allocate memory for arrays */
  memsize = sizeof(double) *
            ((n * (s + s + 3)) + ( s*(4 + s) ) + s*(nrow+ncol) );

  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  ( nrow+ncol x  s)
*
*                                                                       *
*        (memory allocated separately)                                  *
*                                                                       *
*        iwork(iw   )           s  ( s x  2)                            *
*        lwork                  s                                       *
*                                                                       *
***********************************************************************/

  if(!(workptr1 = (double *)malloc(memsize))) {
    perror("THIRD MALLOC FAILED IN TMS2()");
    exit(errno);
  }

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

  if(!(workptr2 = (double **)malloc(memsize))) {
    perror("FOURTH MALLOC FAILED IN TMS2()");
    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    = (nrow+ncol) * s;
  yy        = tempptr1;
  tempptr1 += length;
  y         = tempptr2;
  tempptr2 += s;
  j=0;
  for(i=0; i<length; i+=(nrow+ncol)) y[j++] = &yy[i];


/* Allocate memory for logical array lwork ( long 0= false and 1=true) *
   and integer array iwork                                            */

  memsize = sizeof(long) * ( s * (2 + 1));
  mem   += memsize;

  if(!(workptr3 = (long *)malloc(memsize)) ) {
     perror("FIFTH MALLOC FAILED IN TMS2()");
     exit(errno);
    }

  memsize = sizeof(long *) * 2;
  mem   += memsize;

  if(!(workptr4 = (long **)malloc(memsize))) {
     perror("SIXTH MALLOC FAILED IN TMS2()");
     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;

   /* estimate alpha (via matrix 1-norm)
    * used for uppper bound on max singular value of matrix A */
   alpha = norm_1();
   exetime = timer();

   /* make a trace min. run; exit upon error */ 
   ierr = tsvd2(fp_out1,n,p,s,job,tol,red, sig, maxi, &mem, itcgt,
                titer, time, res, mxv, work1, work2, work3, work4, work5,
                y, iwork, lwork,&maxd);
   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])); 

   if(job!=2) maxd = 0;

   /* write output */

   fprintf(fp_out1, " ... \n");
   fprintf(fp_out1, " ... TRACE MINIMIZATION ON THE  \n");
   fprintf(fp_out1, " ... EQUIVALENT A^TA  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, " ... MAX. DEG. CHEBYSHEV POLY.   =%10ld\n",maxd);
   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. TSVD2       =%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]);