Linear Algebra and the C Language/a0d7


Install and compile this file in your working directory.

/* ------------------------------------ */
/*  Save as:  c00c.c                    */
/* ------------------------------------ */
#include "v_a.h"
/* ------------------------------------ */
#define   RCA        RC4
/* ------------------------------------ */
/* ------------------------------------ */
double f(
double x)
{
        return(1.0/(x));
}
char  feq[] =  "1.0/(x)";
/* ------------------------------------ */
/* ------------------------------------ */
void fun(void)
{
double **A         = rsymmetric_mR(            i_mR(RCA,RCA),9.);
double **EVector   =     eigs_V_mR(A,          i_mR(RCA,RCA));
double **EVectorT  =  transpose_mR(EVector,    i_mR(RCA,RCA)); 
double **T         =        mul_mR(EVectorT,A, i_mR(RCA,RCA));        
double **EValue    =        mul_mR(T,EVector,  i_mR(RCA,RCA));
	
double **invA      =                           i_mR(RCA,RCA);   
double **invEValue =                           i_mR(RCA,RCA);

  clrscrn();
  printf(" A:");
  p_mR(A,S10,P2,C6);  
  
/* EValue = EVectorT A EVector           */  
  f_eigs_mR(f,EValue,invEValue);

/*     invA = EVector invEValue EVectorT */ 
  mul_mR(EVector,invEValue,T);
  mul_mR(T,EVectorT,invA);
  
  printf(" invA = EVector invEValue EVectorT");  
  pE_mR(invA,S10,P4,C6); 
  
  printf(" Ide = A invA "); 
  p_mR(mul_mR(A,invA,T),S10,P5,C10);     
   
  f_mR(A);
  f_mR(invA);
    
  f_mR(EVector);
  f_mR(EVectorT);
  
  f_mR(EValue);
  f_mR(invEValue);    
     
  f_mR(T);
}
/* ------------------------------------ */
int main(void)
{
time_t t;

  srand(time(&t));

do
{
  fun();

} while(stop_w());

  return 0;
}
/* ------------------------------------ */
/* ------------------------------------ */

Screen output example:

                                                                                       
 A:
     +8.00      -2.00      +4.00      -6.00 
     -2.00      +3.00      +6.00      -6.00 
     +4.00      +6.00      +2.00      -5.00 
     -6.00      -6.00      -5.00      -1.00 

 invA = EVector invEValue EVectorT
+6.4696e-02 -8.1470e-02 +2.7157e-02 -3.5144e-02 
-8.1470e-02 -8.5197e-03 +1.1395e-01 -2.9819e-02 
+2.7157e-02 +1.1395e-01 -1.4909e-01 -1.0117e-01 
-3.5144e-02 -2.9819e-02 -1.0117e-01 -1.0437e-01 

 Ide = A invA 
  +1.00000   -0.00000   +0.00000   +0.00000 
  +0.00000   +1.00000   +0.00000   +0.00000 
  +0.00000   +0.00000   +1.00000   +0.00000 
  -0.00000   -0.00000   +0.00000   +1.00000 


 Press   return to continue
 Press X return to stop