lu.c

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

/* this used to be clib/matrix.c, is now lu.c
  several routines have been left behind in clib */

#include "r.h"
#include "mynr.h"

int allocate_dmatrix_family ( dmatrix_family *s , int n , int control ) {
  int status = 0 ; 
  s->n = n ; 
  s->m  = dmatrix ( 1 , n ,  1 , n ) ; 
  if ( control ) { s->in = dmatrix ( 1 , n ,  1 , n ) ; 
		   s->in_allocated = 1 ; }
  else { s->in_allocated = 0 ; }
  s->lu = dmatrix(1,n,1,n);
  s->indx=ivector(1,n);
  s->in_computed = 0 ; 
  s->lu_computed = 0 ; 
  s->det = 0.0 ; 
  return status ;
}

int free_dmatrix_family ( dmatrix_family *s ) {
  int status = 0 ; 
  free_dmatrix ( s->m , 1 , s->n ,  1 , s->n ) ; 
  if ( s->in_allocated ) free_dmatrix ( s->in , 1 , s->n ,  1 , s->n ) ; 
  free_dmatrix ( s->lu , 1 , s->n ,  1 , s->n ) ; 
  free_ivector ( s->indx , 1 , s->n ) ; 
  return status ;
}

int lu_invert ( dmatrix_family *s , double *v ) {
  int status = 0 ; 
  lubksb(s->lu,s->n,s->indx,v);
  return status ;
}

void show_dmatrix_family ( dmatrix_family *s , FILE * fp ) 
{
  int i,j, n=s->n ;
  double **m = s->m ;
  double **in = s->in ;

  for (i=1; i<=n; i++){
    for (j=1; j<=n; j++){
      fprintf(fp,"%6.1g ",m[i][j]);
    }
    fnewline;
  }     
  if ( s->in_computed ) {
    for (i=1; i<=n; i++){
      for (j=1; j<=n; j++){
	fprintf(fp,"%6.1g ",in[i][j]);
      }
      fnewline;
    }     
  }
}

int write_lu ( dmatrix_family *s , char *file ) 
/* writes the lu portion only */
/* was void writelumatrix (m,n,indx,file) */
{
  int i,j,status = 0 , n=s->n ;
  double **m = s->lu ;
  int *indx = s->indx ; 
  FILE    *fp;

  fp = fopen( file, "w" );
  if( !fp )   fprintf( stderr, "No such file: %s\n", file ), status -- ;
  else {
    
    for (i=1; i<=n; i++){
      for (j=1; j<=n; j++){
	fprintf(fp,"%9g ",m[i][j]);
      }
      fnewline;
    }     
    for (i=1; i<=n; i++)
      fprintf(fp,"%d ",indx[i]);
    fnewline;
    fclose(fp);              
  }
  return status ; 
}


int readinlumatrix /* compatible with old programs */
/* ( double **b, int *indx, int n, char *file) */
( dmatrix_family *s , char *file ) {
  int	i, j , status = 0 ;
  FILE    *fp;
  double **b = s->lu ; 
  int *indx = s->indx ; 
  int n = s->n ; 
  
  fp = fopen( file, "r" );
  if( !fp ) {
    fprintf( stderr, "No such file: %s\n", file ), status -- ;
    s->lu_computed = 0 ; }
  else {
    printf( "reading in matrix\n" );
    for (i=1; i<=n; i++){
      for (j=1; j<=n; j++){
	if ( fscanf(fp,"%lf ",&b[i][j]) == EOF ) {
	  status = -1 ;
	  break ;
	}
      }
    }
    if ( status == 0 ) {
      for (i=1; i<=n; i++){
	if ( fscanf(fp,"%d ",&indx[i]) == EOF ) {
	  status = -1 ;
	  break ; 
	}
      }
    }
    fclose( fp );
    if ( status == 0 ) { printf( "lu matrix in\n" ); s->lu_computed = 1 ; }
    else fprintf( stderr, 
		 "Warning: readinlumatrix failed at component %d\n",i);
  }
  return status;
}

double lu_quadratic_form ( double *v1, dmatrix_family *s ,double *v2 ) {
 /* NB v2 must be scratchable */
  return lumatrixproduct( v1, s->lu , v2 , s->n , s->indx ) ;
}

double	lumatrixproduct(double *v1,double **m,double *v2,int n,int *indx)
 /* NB v2 must be scratchable */
{
  int i;
  double d=0.0;
  
  lubksb(m,n,indx,v2);
  for(i=1;i<=n;i++) 
    d+= v1[i]*v2[i];
  return d;
}

#define	float	double
#define	free_vector free_dvector 
#define	vector dvector 
#include "ludcmp.c"
#include "lubksb.c"
#undef	float
#undef	free_vector
#undef	vector 

int	invert_dmatrix_family( dmatrix_family *mp , int control )
 /* returns the determinant of m in mp->det */ 
{
  int n;
  int i,j;
  double d=1.0,*col, **in=mp->in;
  int status = 0 ; 

  n=mp->n;
  col=dvector(1,n);
  dmatrixfromdmatrix(mp->lu,1,n,1,n,mp->m);
  
  status += ludcmp(mp->lu,n,mp->indx,&d); /* see NR page 46 */
  if ( status >= 0 ) {
    mp->lu_computed  =  1 ; 

    for(j=1;j<=n;j++) d*= mp->lu[j][j];
    
    if ( mp->in_allocated && ( control > 0 )) {
      for(j=1;j<=n;j++) {
	for(i=1;i<=n;i++) col[i]=0.0;
	col[j]=1.0;
	lubksb(mp->lu,n,mp->indx,col);
	for(i=1;i<=n;i++) in[i][j]=col[i];
      }
      mp->in_computed  =  1 ; 
    }
    mp->det = d;
  } else {
    mp->det = 0.0;
  }

  free_dvector(col,1,n);
  return status ;
}


void	symmetrize_dmatrix(double **m, int n)
{
  int	i,j;
  
  for(j=2;j<=n;j++) 
    for(i=1;i<j;i++) 
      m[i][j]=m[j][i]=0.5*(m[i][j]+m[j][i]);
}