Linear Algebra and the C Language/a0id

Install and compile this file in your working directory.

/* ------------------------------------ */
/*  Save as :  c00e.c                   */
/* ------------------------------------ */
#include "v_a.h"
/* ------------------------------------ */
void fun(int r)
{
double **A        =  rdefinite_positive_mR( i_mR(r,r), 99);
double **SqrtA    =                         i_mR(r,r);
double **invSqrtA =                         i_mR(r,r);

double **V        =                         i_mR(r,r);
double **VT       =                         i_mR(r,r);

double **EValue   =                         i_mR(r,r);
double **f_EValue =                         i_mR(r,r);

double **T1       =                         i_mR(r,r);
  
  clrscrn();
  printf(" A:");
  p_mR(A, S10,P4, C8); 

/* Eigenvector */   
  eigs_V_mR(A,V); 
  transpose_mR(V,VT);
      
/* EigenValue = VT * A * V */
  mul_mR(VT,A,T1);
  mul_mR(T1,V,EValue);     
     
  printf(" sqrt(A) = V * sqrt(EValue) * VT\n");
  f_eigs_mR(sqrt,EValue,f_EValue);
     mul_mR(V,f_EValue,T1);
     mul_mR(T1,VT, SqrtA); 
       p_mR(SqrtA, S10,P4, C8);    
  
  printf(" inv(sqrt(A)):");  
  inv_mR(SqrtA,invSqrtA);
  p_mR(invSqrtA, S10,P4, C8);   
  
  printf(" IDE = sqrt(A) inv(sqrt(A))");  
  mul_mR(SqrtA,invSqrtA,T1);
  p_mR(T1, S10,P4, C8);  
   
  f_mR(A);
  f_mR(SqrtA);
  f_mR(invSqrtA);
      
  f_mR(V);
  f_mR(VT);
  
  f_mR(EValue);
  f_mR(f_EValue);    
     
  f_mR(T1);
}
/* ------------------------------------ */
int main(void)
{
time_t t;

  srand(time(&t));

do
{
  fun(RC3);

} while(stop_w());

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

A positive-matrix is positive definite if and only if its positive square root is invertible.

Screen output example:

                                                                                       
 A:
  +82.6521    -3.8208    +5.0643 
   -3.8208   +45.6841   +21.6262 
   +5.0643   +21.6262   +68.6638 

 sqrt(A) = V * sqrt(EValue) * VT

   +9.0817    -0.2736    +0.3173 
   -0.2736    +6.5908    +1.4732 
   +0.3173    +1.4732    +8.1482 

 inv(sqrt(A)):
   +0.1105    +0.0058    -0.0053 
   +0.0058    +0.1584    -0.0289 
   -0.0053    -0.0289    +0.1282 

 IDE = sqrt(A) inv(sqrt(A))
   +1.0000    -0.0000    -0.0000 
   -0.0000    +1.0000    -0.0000 
   +0.0000    +0.0000    +1.0000 


 Press   return to continue
 Press X return to stop