Linear Algebra and the C Language/a0lk

/* ------------------------------------ */
/*  Save as :   c03d.c                  */
/* ------------------------------------ */
#include "v_a.h"
/* ------------------------------------ */
/* ------------------------------------ */
#define    RB R3
#define    CB C1
/* ------------------------------------ */
/* ------------------------------------ */
int main(void)
{
double b[RB*(CB)]={
  +0.333333333333, 
  +0.666666666667, 
  +1.000000000000 
};

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

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

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

double **X           = ca_A_mR(x,i_mR(RB,C1));
double **VX          =           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(transpose_mR(B,BT),S5,P4,C7);
  
  printf(" BT B:");
  p_mR(mul_mR(BT,B,BTB),S5,P4,C7);
  
  printf(" inv(BTB):");
  invBTB[R1][C1] = 1./BTB[R1][C1];
  p_mR(invBTB,S5,P4,C7);  
  
  printf(" inv(BTB) BT :");
  p_mR(mul_mR(invBTB,BT,invBTB_BT),S5,P4,C7); 
  
  printf(" B inv(BTB)BT:");
  p_mR(mul_mR(B,invBTB_BT,B_invBTB_BT),S5,P4,C7); 
   
  printf(" V = Id-(B inv(BTB)BT)");
  p_mR(sub_mR(Id,B_invBTB_BT,V),S5,P4,C7);    
  stop();  
  
  clrscrn();
  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(" X:");
  p_mR(X,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(mul_mR(V,X,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 R3  :                           

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

 B:         See the previous c file
+0.3333 
+0.6667 
+1.0000 

 Press return to continue. 


 BT:
+0.3333 +0.6667 +1.0000 

 BT B:
+1.5556 

 inv(BTB):
+0.6429 

 inv(BTB) BT :
+0.2143 +0.4286 +0.6429 

 B inv(BTB)BT:
+0.0714 +0.1429 +0.2143 
+0.1429 +0.2857 +0.4286 
+0.2143 +0.4286 +0.6429 

 V = Id-(B inv(BTB)BT)
+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  

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

 X:
-1.0000 
+2.0000 
-3.0000 

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

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

 Press return to continue.