test_lu.c

来自「快速傅立叶变换程序代码,学信号的同学,可要注意了」· C语言 代码 · 共 79 行

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

/* 
   test program for matrix inversion.
*/
   
void main(int argc, char *argv[])
{
  dmatrix_family param , other ;
  double *x , *b , *b2 ;
  int i , n ;
  double d ;
  char file[100]="tmpfile" ; 

  /* Load up the parameters of the function that you want to optimize */
  printf("Solving A x = b\nDimension of A (eg 2)?\n");
  inputi(&n);
  allocate_dmatrix_family ( &param , n , 1 ) ;

  b=dvector(1,n);
  b2=dvector(1,n);
  x=dvector(1,n);

  typeindmatrix(param.m,1,n,1,n);
  printf("b vector?\n");
  typeindvector(b,1,n);

  d = invert_dmatrix_family ( &param, 1 ) ; /* control integer = 1
					       causes inverse to be computed */
  
  /* to solve x = A^-1 b, copy b into x */

  dvectorfromdvector( x , 1 , n , b ) ; 
  lu_invert ( &param , x ) ; 

  printf("Solution:\n");
  for ( i = 1 ; i <= n ; i ++ ) {
    printf ( "%10.5g " , x[i] ) ; 
  }
  printf ( "\n" ) ; 

  printf ( "writing to %s\n", file ) ; 
  write_lu ( &param , file ) ; 
  printf ( "done\nreading again\n" ) ; 
  readinlumatrix ( &param , file ) ; 
  printf ( "\ndone\n" ) ; 

  dvectorfromdvector( x , 1 , n , b ) ; 
  lu_invert ( &param , x ) ; 

  printf("Solution:\n");
  for ( i = 1 ; i <= n ; i ++ ) {
    printf ( "%10.5g " , x[i] ) ; 
  }
  printf ( "\n" ) ; 

  printf ("Now going back\n" ) ; 
  allocate_dmatrix_family ( &other , n , 1 ) ;
  dmatrixfromdmatrix ( other.m , 1 , n , 1 , n , param.in ) ; 
  d = invert_dmatrix_family ( &other , 1 ) ; /* control integer = 1
					       causes inverse to be computed */
  lu_invert ( &other , x ) ; /* x should come back to b */
  printf("Solution should equal b:\nx     b\n");
  for ( i = 1 ; i <= n ; i ++ ) {
    printf ( "%10.5g " , x[i] ) ; 
    printf ( "%10.5g " , b[i] ) ; 
    printf ( "\n" ) ; 
  }
  printf("here is the inverse \n" );
  printoutdmatrix ( param.in , 1,n,1,n, 94 ) ;
    printf ( "\n" ) ; 
  printf("here is the inverse of the inverse - should equal A \n" );
  printoutdmatrix ( other.in , 1,n,1,n, 94 ) ;
    printf ( "\n" ) ; 

  exit(0);
}