Linear Algebra and the C Language/a0lh

/* ------------------------------------ */
/*  Save as:   c03a.c                   */
/* ------------------------------------ */
#include "v_a.h"
/* ------------------------------------ */
/* ------------------------------------ */
#define   RA R3
#define   CA C2
/* ------------------------------------ */
/* ------------------------------------ */
int main(void)
{
double a[RA*(CA)]={
 -2, -3,
  1,  0,
  0,  1
};

double x[RA*C1]={
 -1, 
  2, 
 -3
};

double **A         =      ca_A_mR(a,          i_mR(RA,CA));
double **AT        = transpose_mR(A,          i_mR(CA,RA));
double **ATA       =       mul_mR(AT,A,       i_mR(CA,CA)); 
double **invATA    =     invgj_mR(ATA,        i_mR(CA,CA));
double **invATA_AT =       mul_mR(invATA,AT,  i_mR(CA,RA));
double **V         =       mul_mR(A,invATA_AT,i_mR(RA,RA));
                               /* A inv(AT A) AT */ 

double **X         =      ca_A_mR(x,          i_mR(RA,C1));
double **VX        =       mul_mR(V,X,        i_mR(RA,C1));

  clrscrn();
  printf(" A is subspace of R%d             \n\n"
         " Find a transformation matrix for   \n"
         " a projection onto R%d.           \n\n"
         " Proj(x) =  A inv(AT A) AT        \n\n",RA,RA);
         
  printf(" A:");
  p_mR(A,S5,P1,C7);
  stop();
  
  clrscrn();
  printf(" AT:");
  p_mR(AT,S5,P1,C7);
  
  printf(" AT A:");
  p_mR(ATA,S5,P1,C7);
  
  printf(" inv(ATA):");
  p_mR(invATA,S5,P4,C7); 
   
  printf(" inv(ATA)AT:");
  p_mR(invATA_AT,S5,P4,C7); 
  
  printf(" V = A inv(ATA)AT");
  p_mR(V,S5,P4,C7);    
  stop();  
  
  clrscrn();
  printf(" V is transformation matrix for     \n"
         " a projection onto a R%d subspace \n\n"
         " V:",RA);
  p_mR(V,S5,P4,C7); 
  
  printf(" X:");
  p_mR(X,S5,P1,C7);
  
  printf(" Proj(x) =  A inv(ATA)AT x \n\n"  
         " Proj(x) =  V x      Verify the result with the last C file");  
  p_mR(VX,S5,P4,C7); 
  stop();    
  
  f_mR(A);
  f_mR(AT);
  f_mR(ATA);       
  f_mR(invATA);    
  f_mR(invATA_AT); 
  f_mR(V);         /* A inv(AT A) AT */   
  
  f_mR(X); 
  f_mR(VX);         

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

                                                                                       
 A is subspace of R3             

 Find a transformation matrix for   
 a projection onto R3 :          

 Proj(x) =  A inv(AT A) AT        

 A:
 -2.0  -3.0 
 +1.0  +0.0 
 +0.0  +1.0 

 Press return to continue. 


 AT:
 -2.0  +1.0  +0.0 
 -3.0  +0.0  +1.0 

 AT A:
 +5.0  +6.0 
 +6.0 +10.0 

 inv(ATA):
+0.7143 -0.4286 
-0.4286 +0.3571 

 inv(ATA)AT:
-0.1429 +0.7143 -0.4286 
-0.2143 -0.4286 +0.3571 

 V = A inv(ATA)AT
+0.9286 -0.1429 -0.2143 
-0.1429 +0.7143 -0.4286 
-0.2143 -0.4286 +0.3571 

 Press return to continue. 


 V is transformation matrix for      
 a projection onto a R3 subspace.


+0.9286 -0.1429 -0.2143 
-0.1429 +0.7143 -0.4286 
-0.2143 -0.4286 +0.3571 

 X:
 -1.0 
 +2.0 
 -3.0 

 Proj(x) =  A inv(ATA)AT x 

 Proj(x) =  V x      Verify the result with the last C file
-0.5714 
+2.8571 
-1.7143 

 Press return to continue.