Linear Algebra and the C Language/a0dq

/* ------------------------------------ */
/*  Save as :   c00c.c                  */
/* ------------------------------------ */
#include "v_a.h"
/* ------------------------------------ */
/* ------------------------------------ */
void fun(int r,int c)
{
double **B       =         r_mR(     i_mR(c,r),9.);
double **B_T     = transpose_mR(B,   i_mR(r,c));

double **U       =  svd_U_Rn_mR(B_T, i_mR(c,r));  
double **U_T     = transpose_mR(U,   i_mR(r,c));  
double **U_TB    =                   i_mR(r,r);       /* rc cr        */
double **V       =  svd_V_Rn_mR(B_T, i_mR(r,r));   
double **VS      =                   i_mR(r,r);       /* rr rr        */

double **S       =                   i_mR(r,r);
double **invS    =                   i_mR(r,r);
                                               /* Pinv = V  invS U_T  */        
double **Pinv    =                   i_mR(r,c);       /* rr  rr   rc  */
double **PinvB   =                   i_mR(r,r);       /* rr rr rc cr  */

/* S = U_T  B  V     */ 
  mul_mR(U_T,B,U_TB);
  mul_mR(U_TB,V,S);

/* invS              */  
  inv_svd_mR(S,invS);
  
/* Pinv = V invS U_T */    
  mul_mR(V,invS,VS);
  mul_mR(VS,U_T,Pinv);    

/* Ide = Pinv B */ 
  mul_mR(Pinv,B,PinvB);

  clrscrn();
  printf(" Copy/Paste into the octave windows \n\n\n");
  p_Octave_mR(B,"B",P2);
  printf("pinv(B)\n\n\n");
  stop();

  clrscrn();    
  printf(" Pinv = V invS U_T "); 
  p_mR(Pinv,S10,P6,C10);  
  
  printf(" Ide = Pinv B ");      
  p_mR(PinvB,S10,P6,C10);   
    
  f_mR(B);
  f_mR(B_T);
  f_mR(U);
  f_mR(U_T);
  f_mR(U_TB);  
  f_mR(V);  
  f_mR(VS); 
  f_mR(S); 
  f_mR(invS);
  f_mR(Pinv);
  f_mR(PinvB); 
}
/* ------------------------------------ */
int main(void)
{
time_t t;

  srand(time(&t));

int i;

do
{
  i = rp_I(R3)+R1;
  
  fun(i,i+R2);

} while(stop_w());

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

                                                                                       
 Copy/Paste into the octave windows 


 B=[
+7.00,+2.00,-6.00;
+9.00,+7.00,-5.00;
+1.00,-8.00,-6.00;
-9.00,-1.00,+3.00;
-8.00,+5.00,-1.00]

pinv(B)


 Press return to continue. 


 Pinv = V invS U_T 
 +0.002456  +0.015315  -0.022833  -0.037373  -0.066431 
 +0.010749  +0.043518  -0.053964  +0.002403  +0.048913 
 -0.052801  -0.027851  -0.083505  -0.010925  -0.075681 

 Ide = Pinv B 
 +1.000000  +0.000000  +0.000000 
 -0.000000  +1.000000  -0.000000 
 +0.000000  +0.000000  +1.000000 


 Press   return to continue
 Press X return to stop