411 REFERENCE · Wikibooks
Linear Algebra and the C Language/a0e0
Reference material by Wikibooks contributors. This retains its source identity and is not 411 news reporting.
__NOTOC__
Code study: svd_V_Cn_mR(); svd_U_Cn_mR();
+ svd_V_Cn_mR(); svd_U_Cn_mR();
/* -------------------------------- */
/* -------------------------------- */
double **svd_V_Cn_mR(
double **A,
double **V
)
{
int i;
int j;
int r = rsize_R(A);
int c = csize_R(A);
double **tA = i_mR(c,r);
double **tAA = i_mR(c,c);
double **AtA = i_mR(r,r);
double **EValue = i_mR(c,C1);
double **EV = i_mR(c,c);
double **Ab;
double **Ab0;
double **b = i_mR(c,C1);
double **Ide = i_mR(c,c);
double **sIde = i_mR(c,c);
double **tAAsIde = i_mR(c,c);
transpose_mR(A,tA) ;
mul_mR(A,tA,AtA);
eigs_mR(AtA,EValue);
mul_mR(tA,A,tAA);
eye_mR(Ide);
for(i=R1; i<= rsize_R(EValue); i++)
{
if(EValue[i][C1])
{
smul_mR(EValue[i][C1],Ide,sIde);
sub_mR(tAA,sIde,tAAsIde);
Ab = c_A_b_Ab_mR(tAAsIde,b,
i_Abr_Ac_bc_mR(c,c,C1));
GJ_PP_FreeV_mR(Ab,EV);
c_c_mR(EV,C2,V,i);
}
else
{
/* smul_mR(EValue[i][C1]==0,Ide,sIde);
sub_mR(tAA,sIde,tAAsIde); */
Ab0 = c_A_b_Ab_mR(tAA,b,
i_Abr_Ac_bc_mR(c,c,(C1+c-r)));
GJ_PP_FreeV_Value0_mR(Ab0,EV);
for(j=C2; i<=rsize_R(EValue); i++,j++)
c_c_mR(EV,j,V,i);
}
}
Normalize_mR(V);
f_mR(EValue);
f_mR(EV);
f_mR(tA);
f_mR(tAA);
f_mR(AtA);
f_mR(Ab);
f_mR(Ab0);
f_mR(b);
f_mR(Ide);
f_mR(sIde);
f_mR(tAAsIde);
return(V);
}
/* -------------------------------- */
/* -------------------------------- */
/* -------------------------------- */
/* -------------------------------- */
double **svd_U_Cn_mR(
double **A,
double **U
)
{
int i;
int r = rsize_R(A);
int c = csize_R(A);
double **tA = i_mR(c,r);
double **AtA = i_mR(r,r);
double **EValue = i_mR(r,C1);
double **EV = i_mR(r,r);
double **Ab;
double **b = i_mR(r,C1);
double **Ide = i_mR(r,r);
double **sIde = i_mR(r,r);
double **AtAsIde = i_mR(r,r);
transpose_mR(A,tA) ;
mul_mR(A,tA,AtA);
eigs_mR(AtA,EValue);
eye_mR(Ide);
for(i=R1; i<=rsize_R(EValue); i++)
{
smul_mR(EValue[i][C1],Ide,sIde);
sub_mR(AtA,sIde,AtAsIde);
Ab = c_A_b_Ab_mR(AtAsIde,b,
i_Abr_Ac_bc_mR(r,r,C1));
GJ_PP_FreeV_mR(Ab,EV);
c_c_mR(EV,C2,U,i);
}
Normalize_mR(U);
f_mR(EValue);
f_mR(EV);
f_mR(tA);
f_mR(AtA);
f_mR(Ab);
f_mR(b);
f_mR(Ide);
f_mR(sIde);
f_mR(AtAsIde);
return(U);
}
/* -------------------------------- */
/* -------------------------------- */