Linear Algebra and the C Language/a0lp

/* ------------------------------------ */
/*  Save as:   c04d.c                   */
/* ------------------------------------ */
#include "v_a.h"
/* ------------------------------------ */
/* ------------------------------------ */
#define   RB R4
#define   CB C2
/* ------------------------------------ */
/* ------------------------------------ */
int main(void)
{
double b[RB*(CB)]={
 +0.000000000000, -2.000000000000, 
 +0.000000000000, -1.666666666667, 
 +1.000000000000, +0.000000000000, 
 +0.000000000000, +1.000000000000
};

double x[RB*(C1)]={
 -1, 
  2, 
 -3,
 -3,
};

double **B           =      ca_A_mR(b,              i_mR(RB,CB));
double **BT          = transpose_mR(B,              i_mR(CB,RB));
double **BTB         =       mul_mR(BT,B,           i_mR(CB,CB));  
double **invBTB      =       inv_mR(BTB,            i_mR(CB,CB)); 
double **invBTB_BT   =       mul_mR(invBTB,BT,      i_mR(CB,RB));   
double **B_invBTB_BT =       mul_mR(B,invBTB_BT,    i_mR(RB,RB));  

double **Id          =       eye_mR(                i_mR(RB,RB));
double **V           =       sub_mR(Id,B_invBTB_BT, i_mR(RB,RB)); 
                                 /* Id-(B inv(BT B) BT) */

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

  clrscrn();
  printf(" B is a basis for the orthogonal complement of AT: \n\n"
         " Find a transformation matrix for                    \n"
         " a projection onto R%d:                            \n\n"
         " Proj(x) = [Id-(B inv(BTB)BT)]                     \n\n"
         " B:         See the previous c file",RB);
  p_mR(B,S5,P4,C7);
  stop();
  
  clrscrn();
  printf(" BT:");
  p_mR(BT,S5,P4,C7);
  printf(" BT B:");
  p_mR(BTB,S5,P4,C7);
  
  printf(" inv(BTB):");
  p_mR(invBTB,S5,P4,C7);   
  
  printf(" inv(BTB) BT:");
  p_mR(invBTB_BT,S5,P4,C7); 
  
  printf(" B inv(BTB)BT:");
  p_mR(B_invBTB_BT,S5,P4,C7);
  stop();  
  
  clrscrn();    
  printf(" V = Id - (B inv(BT B) BT)");
  p_mR(V,S5,P4,C7);    

  printf(" V is transformation matrix for     \n"
         " a projection onto a R%d subspace \n\n"
         " V:",RB);
  p_mR(V,S5,P4,C7); 
  
  printf(" Proj(x) = [Id-(B inv(BTB)BT)] x \n\n" 
         " Proj(x) =  V x     Verify the result with the first C file");
  p_mR(VX,S5,P4,C7); 
  stop();    
  
  
  f_mR(B);
  f_mR(BT);
  f_mR(BTB);       
  f_mR(invBTB);    
  f_mR(invBTB_BT); 
  f_mR(V);       /* Id-(B inv(BT B) BT) */   
  
  f_mR(X); 
  f_mR(VX);         

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

                                                                                       
 B is a basis for the orthogonal complement of AT: 

 Find a transformation matrix for                    
 a projection onto R4:                            

 Proj(x) = [Id-(B inv(BTB)BT)]                     

 B:         See the previous c file
+0.0000 -2.0000 
+0.0000 -1.6667 
+1.0000 +0.0000 
+0.0000 +1.0000 

 Press return to continue. 


 BT:
+0.0000 +0.0000 +1.0000 +0.0000 
-2.0000 -1.6667 +0.0000 +1.0000 

 BT B:
+1.0000 +0.0000 
+0.0000 +7.7778 

 inv(BTB):
+1.0000 -0.0000 
-0.0000 +0.1286 

 inv(BTB) BT:
+0.0000 +0.0000 +1.0000 +0.0000 
-0.2571 -0.2143 +0.0000 +0.1286 

 B inv(BTB)BT:
+0.5143 +0.4286 +0.0000 -0.2571 
+0.4286 +0.3571 +0.0000 -0.2143 
+0.0000 +0.0000 +1.0000 +0.0000 
-0.2571 -0.2143 +0.0000 +0.1286 

 Press return to continue. 


 V = Id - (B inv(BT B) BT)
+0.4857 -0.4286 +0.0000 +0.2571 
-0.4286 +0.6429 +0.0000 +0.2143 
+0.0000 +0.0000 +0.0000 +0.0000 
+0.2571 +0.2143 +0.0000 +0.8714 

 V is transformation matrix for     
 a projection onto a R4 subspace  

 V:
+0.4857 -0.4286 +0.0000 +0.2571 
-0.4286 +0.6429 +0.0000 +0.2143 
+0.0000 +0.0000 +0.0000 +0.0000 
+0.2571 +0.2143 +0.0000 +0.8714 

 Proj(x) = [Id-(B inv(BTB)BT)] x 

 Proj(x) =  V x     Verify the result with the first C file
-2.1143 
+1.0714 
+0.0000 
-2.4429 

 Press return to continue.