Mathc matrices/h12b
Apparence
Installer ce fichier dans votre répertoire de travail.
vdsvdrn.h |
---|
/* ------------------------------------ */
/* Save as : vdsvdrn.h */
/* ------------------------------------ */
/* ------------------------------------ */
double **svd_V_Rn_mR(
double **A,
double **V
)
{
int i;
int r = rsize_R(A);
int c = csize_R(A);
double **A_T = i_mR(c,r);
double **AA_T = i_mR(r,r);
double **EigsValue = i_mR(r,C1);
double **eignV = i_mR(r,r);
double **Ab;
double **b = m0_mR(i_mR(r,C1));
double **Ide = eye_mR(i_mR(r,r));
double **sIde = i_mR(r,r);
double **AA_TmnssIde = i_mR(r,r);
/* ------------------------------------
eigs_value_AA_T_mR(A, eignValue); if(c > r)
eigs_value_AA_T_mR(A_T,eignValue); if(r > c)
------------------------------------ */
eigs_value_AA_T_mR(A,EigsValue);
transpose_mR(A,A_T) ;
mul_mR(A,A_T,AA_T);
for(i = R1; i <= rsize_R(EigsValue); i++)
{
smul_mR(EigsValue[i][C1],Ide,sIde);
MmnsD_mR(AA_T,sIde,AA_TmnssIde);
Ab = c_A_b_Ab_mR(AA_TmnssIde,b,i_Abr_Ac_bc_mR(r,r,C1));
GJ_PP_FreeV_mR(Ab,eignV);
c_c_mR(eignV,C2,V,i);
}
Normalize_mR(V);
f_mR(EigsValue);
f_mR(eignV);
f_mR(A_T);
f_mR(AA_T);
f_mR(Ab);
f_mR(b);
f_mR(Ide);
f_mR(sIde);
f_mR(AA_TmnssIde);
return(V);
}
/* ------------------------------------ */
/* ------------------------------------ */
double **svd_U_Rn_mR(
double **A,
double **U
)
{
int i;
int r = rsize_R(A);
int c = csize_R(A);
double **A_T = i_mR(c,r);
double **A_TA = i_mR(c,c);
double **EigsValue = i_mR(r,C1);
double **eignV = i_mR(c,r);
double **Ab;
double **b = m0_mR(i_mR(c,C1));
double **Ide = eye_mR(i_mR(c,c));
double **sIde = i_mR(c,c);
double **A_TAmnssIde = i_mR(c,c);
/* ------------------------------------
eigs_value_AA_T_mR(A, eignValue); if(c > r)
eigs_value_AA_T_mR(A_T,eignValue); if(r > c)
------------------------------------ */
eigs_value_AA_T_mR(A,EigsValue);
transpose_mR(A,A_T);
mul_mR(A_T,A,A_TA);
for(i = R1; i <= rsize_R(EigsValue); i++)
{
smul_mR(EigsValue[i][C1],Ide,sIde);
MmnsD_mR(A_TA,sIde,A_TAmnssIde);
Ab = i_Abr_Ac_bc_mR(rsize_R(A_TAmnssIde),
csize_R(A_TAmnssIde),
csize_R(b));
c_A_b_Ab_mR(A_TAmnssIde,b, Ab);
GJ_PP_FreeV_mR(Ab,eignV);
c_c_mR(eignV,C2,U,i);
}
Normalize_mR(U);
f_mR(EigsValue);
f_mR(eignV);
f_mR(A_TA);
f_mR(Ab);
f_mR(b);
f_mR(Ide);
f_mR(sIde);
f_mR(A_TAmnssIde);
return(U);
}
/* ------------------------------------ */
/* ------------------------------------ */
double **pseudo_Rn_mR(
double **A,
double **Pinv,
double factor_E
)
{
int r = rsize_R(A);
int c = csize_R(A);
double **a = smul_mR(factor_E,A,
i_mR(r,c));
double **a_T = transpose_mR(a,
i_mR(c,r));
double **U = i_mR(r,c);
double **U_T = i_mR(c,r);
double **V = i_mR(c,c);
double **U_TA = i_mR(c,c);
double **S = i_mR(c,c);
double **S_T = i_mR(c,c);
double **invS_T = i_mR(c,c);
double **VS_T = i_mR(c,c);
svd_U_Rn_mR(a_T,U);
svd_V_Rn_mR(a_T,V);
/* S = U_T * A * V */
transpose_mR( U,U_T);
mul_mR(U_T,A,U_TA);
mul_mR(U_TA,V,S);
/* invS_T = */
transpose_mR(S,S_T);
inv_svd_mR(S_T,invS_T);
/* PseudoInverse = V * invS_T * U_T */
mul_mR(V,invS_T,VS_T);
mul_mR(VS_T,U_T,Pinv);
f_mR(U);
f_mR(U_T);
f_mR(V);
f_mR(a);
f_mR(a_T);
f_mR(U_TA);
f_mR(S);
f_mR(S_T);
f_mR(invS_T);
f_mR(VS_T);
return(Pinv);
}
/* ------------------------------------ */
/* ------------------------------------ */