Linear Algebra and the C Language/a0bw

/* ------------------------------------ */
/*  Save as :  c00b.c                   */
/* ------------------------------------ */
#include "v_a.h"
/* ------------------------------------ */
#define   RCA        RC4
#define   Pn           3
/* ------------------------------------ */
double f(
double x)
{  
        return(pow(x,1./Pn));
}
char  feq[] =  "x**Pn";
/* ------------------------------------ */
/* ------------------------------------ */
void fun(void)
{
double **A             = rPsymmetric_mR(i_mR(RCA,RCA),99.);
double **SqrtA          =               i_mR(RCA,RCA);

double **EigsVector    =               i_mR(RCA,RCA);
double **T_EigsVector  =               i_mR(RCA,RCA);

double **EigsValue     =               i_mR(RCA,RCA);
double **f_EigsValue   =               i_mR(RCA,RCA);

double **T1            =               i_mR(RCA,RCA);
  
  clrscrn();
  printf(" A :");
  p_mR(A, S7,P0, C8); 
  
  eigs_V_mR(A,EigsVector); 
  transpose_mR(EigsVector,T_EigsVector);
      
/* EigsValue : T_EigsVector * A * EigsVector */
  mul_mR(T_EigsVector,A,T1);
  mul_mR(T1,EigsVector,EigsValue);      
     
  printf(" A**(1/%d) = EigsVector * EigsValue**(1/%d)"
         " * T_EigsVector\n", Pn, Pn);
  f_eigs_mR(f,EigsValue,f_EigsValue);      
     mul_mR(EigsVector,f_EigsValue,T1);
     mul_mR(T1,T_EigsVector,SqrtA);
       p_mR(SqrtA,S10,P4,C8);    
  
  printf(" A = A**%d ",Pn);  
  pow_mR(Pn,SqrtA,A);
  p_mR(A, S7,P0, C8);   

  f_mR(A);
  f_mR(SqrtA);
    
  f_mR(EigsVector);
  f_mR(T_EigsVector);
  
  f_mR(EigsValue);
  f_mR(f_EigsValue);    
     
  f_mR(T1);
}
/* ------------------------------------ */
int main(void)
{
time_t t;

  srand(time(&t));

do
{
  fun();

} while(stop_w());

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

                                                                                       
 A :
 +10434     +35     +65     -32 
    +35  +10323     -85     -19 
    +65     -85   +8658     +53 
    -32     -19     +53   +3330 

 A**(1/3) = EigsVector * EigsValue**(1/3) * T_EigsVector

  +21.8514    +0.0246    +0.0484    -0.0313 
   +0.0246   +21.7736    -0.0633    -0.0183 
   +0.0484    -0.0633   +20.5335    +0.0557 
   -0.0313    -0.0183    +0.0557   +14.9327 

 A = A**3 
 +10434     +35     +65     -32 
    +35  +10323     -85     -19 
    +65     -85   +8658     +53 
    -32     -19     +53   +3330 


 Press   return to continue
 Press X return to stop