Linear Algebra and the C Language/a0dy

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

double **U      =  svd_U_Cn_mR(A, i_mR(r,r)) ;
double **U_T    = transpose_mR(U, i_mR(r,r));
double **U_TA   =                 i_mR(r,c);

double **V      =  svd_V_Cn_mR(A, i_mR(c,c));   
                                                /*    V   S_T         */
double **VS_T   =                 i_mR(c,r);    /*    cc  cr          */

double **S      =                 i_mR(r,c);
double **S_T    =                 i_mR(c,r);
double **invS_T =                 i_mR(c,r); 
                                                /*    V  invS_T U_T   */
double **Pinv   =                 i_mR(c,r);    /*    cc cr     rr    */

                                                /* A  V  invS_T U_T   */
double **Ide    =                 i_mR(r,r);    /* rc cc cr     rr    */

/* S = U_T A V         */ 
  mul_mR(U_T,A,U_TA);
  mul_mR(U_TA,V,S);

/* Pinv = V invS_T U_T */
  transpose_mR(S,S_T);
    inv_svd_mR(S_T,invS_T);
        mul_mR(V,invS_T,VS_T);
        mul_mR(VS_T,U_T,Pinv);  
  
/* Ide = A Pinv        */
  mul_mR(A,Pinv,Ide);      
  
  clrscrn();
  printf(" Copy/Paste into the octave windows \n\n\n");
  p_Octave_mR(A,"A",P2);
  printf("pinv(A)\n\n\n");
  stop();

  clrscrn(); 
  printf(" Pinv = V invS_T U_T "); 
  pE_mR(Pinv,S10,P6,C10); 
  
  printf(" Ide = A Pinv ");
  p_mR(Ide,S10,P6,C10);  
    
  f_mR(A); 
  f_mR(U);
  f_mR(U_T); 
  f_mR(U_TA);
  f_mR(V); 
  f_mR(VS_T);   
  f_mR(S); 
  f_mR(S_T);
  f_mR(invS_T);
  f_mR(Pinv); 
  f_mR(Ide); 
}
/* ------------------------------------ */
int main(void)
{
time_t t;

  srand(time(&t));

int i;

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

} while(stop_w());

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

                                                                                       
 Copy/Paste into the octave windows 


 A=[
-9.00,-2.00,-6.00,-7.00,+6.00,-4.00;
+4.00,+7.00,-8.00,-2.00,-1.00,-4.00;
-4.00,-4.00,+6.00,-6.00,+1.00,-5.00]

pinv(A)


 Press return to continue. 


 Pinv = V invS_T U_T 
-5.107413e-02 +4.175497e-02 +1.868221e-02 
-1.395237e-02 +4.883623e-02 +3.030004e-04 
-4.690935e-02 -2.043016e-02 +6.399132e-02 
-8.021914e-03 -3.586208e-02 -5.829170e-02 
+3.773916e-02 -2.230133e-02 -2.483506e-02 
+1.290410e-02 -5.841493e-02 -7.341556e-02 

 Ide = A Pinv 
 +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