codeeig.c

快速傅立叶变换程序代码,学信号的同学,可要注意了

/* 
   codeeig.c                                   (c) DJCM 95 Sep

   97 Aug

   finds the second eig of a graph assuming it is regular. 
 example usage

 codeeig -ai GHC/n17136 -v 2 -l 200

 This only does one random start so it should be repeated a few times
 with different seeds and a max taken.

*/

#include "./ansi/r.h"
#include "./ansi/rand2.h"
#include "./ansi/mynr.h"
#include "./ansi/cmatrix.h"

typedef struct { /* codeeig_control */
#include "codeeig_var_str.h"
} codeeig_control ;

static void print_usage ( char ** , FILE *  , codeeig_control * ) ;
static void c_defaults ( codeeig_control * ) ; 
static int process_command ( int  , char ** , codeeig_control * ) ;
 
void   main ( int , char ** ) ;

/*
        MAIN
                     */
void main ( int argc, char *argv[] )
{
  alist_matrix a ;
  double *t , *x , sf , ss , isss , sss , s ;
  int n , N , M ;
  int verbose , status=0 , l ;
  codeeig_control c ;

  strcpy ( c.afile , "default" ) ; 
  c_defaults ( &c ) ; 

  if ( process_command (argc, argv, &c ) < 0 ) exit (0) ;
  if ( read_allocate_alist ( &a , c.afile ) < 0 ) exit (0) ; 

  N = a.N ; M = a.M ; verbose = c.verbose ; 

/* start from a random vector (real) satisfying x.n = 0 */
/* propagate it through to t and back; rescale; repeat; */
/* keep track of rescalings. (and whether x.n = 0 )     */

  x =  dvector ( 1 , a.N ) ; 
  t =  dvector ( 1 , a.M ) ; 

  if ( c.zss ) {
    isss = 1.0/sqrt((double)(2*(N/2))) ; 
    for ( n = 1 ; n <= N ; n ++ ) { /* set x arbitrarily */
      x[n] = (n%2) ? isss : -isss ; 
    } 
    if (N%2) {  
      x[N] = 0.0 ; 
    } /* ensure zero sum. */
  } else {
    if ( c.seed ) {
      ran_seed ( c.seed ) ; 
      for ( n = 1 ; n <= N ; n ++ ) { /* set x arbitrarily */
	x[n] = rann()  ; 
      } 
      
    } else {
      isss = 1.0/sqrt((double)((N))) ; 
      for ( n = 1 ; n <= N ; n ++ ) { /* set x arbitrarily */
	x[n] = isss  ; 
      } 
    }
  }
  for ( l = 1 ; l <= c.L ; l ++ ) { /* main loop */
    ss = 0.0 ;  s = 0.0 ; 
    for ( n = 1 ; n <= N ; n ++ ) { 
      ss += x[n] * x[n] ; 
      s += x[n]  ; 
    }
    sf = s / (double)(N) ; 
    if ( c.rzs ) {
      for ( n = 1 ; n <= N ; n ++ ) { 
	x[n] -= sf ; 
      }
    }
    sss = sqrt( ss ) ;  isss = 1.0 / sss ;
    for ( n = 1 ; n <= N ; n ++ ) { 
      x[n] *= isss ; 
    }
    if ( ( verbose && (!(l%c.period)) ) || ( l == c.L ) )  
      printf ( "%-3d\ts: %10.4g\tss: %8.4g\tsss: %8.4g\n" , l , s , ss , sss ) ;
    
    alist_times_dvector ( &a , x , t ) ; 
    alist_transpose_dvector ( &a , t , x ) ; 

  } /* end of r loop */
  if ( verbose >= 2 ) {
    pause_for_return();
    for ( n = 1 ; n <= N ; n ++ ) { 
      printf ( "%10.4g" , x[n] ) ;
      if ( !(n%5) ) printf ( "\n" ) ; 
      else printf ( "\t" ) ; 
    }
    printf ( "\n" ) ; 
  }

}

static void c_defaults ( codeeig_control *c ) 
{
#include "codeeig_var_def.c"
}

static int process_command ( int argc , char **argv , codeeig_control *c ) {
  int p_usage = 0 ;
  int status = 0 ;
  int cs , i ;

  if ( argc < 1 )     {
    p_usage = 1 ; 
    status -- ;
  }

#define ERROR1 fprintf ( stderr , "arg to `%s' missing\n" , argv[i] ) ; \
               status --
#define ERROR2 fprintf ( stderr , "args to `%s' missing\n" , argv[i] ) ; \
               status --
#define ERRORREG fprintf ( stderr , "regtype must be defined before `%s'\n" , argv[i] ) ; \
               status --
  for (i = 1 ; i < argc; i++)    {
    cs = 1 ;
    if ( strcmp (argv[i], "-V") == 0 )        {
      c->verbose = 1;
    }
#include "codeeig_var_clr.c"
    else {
      fprintf ( stderr , "arg `%s' not recognised\n" , argv[i] ) ; 
      p_usage = 1 ;
      status -- ;
    }
    if ( cs == 0 ) {
      fprintf ( stderr , "arg at or before `%s' has incorrect format\n" , 
	       argv[i] ) ;
      p_usage = 1 ;
      status -- ;
    }
  }
  if ( p_usage ) print_usage ( argv , stderr , c ) ;
  return ( status ) ;
}
#undef ERROR1
#undef ERROR2
#undef ERRORREG

#define DNT fprintf( fp, "\n        ")
#define NLNE  fprintf( fp, "\n")

static void print_usage ( char **argv , FILE * fp , codeeig_control *c )
{
  
  fprintf( fp, "Usage: %s ",argv[0]);
  fprintf( fp, " [optional arguments]");

#include "codeeig_var_usg.c"
  return ;
}

/*
<!-- hhmts start -->
Last modified: Mon Dec  9 23:14:05 1996
<!-- hhmts end -->
*/