Linear Algebra and the C Language/a02t

__NOTOC__ Install this file in your working directory. /* ------------------------------------ */ /* 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 **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,V,i); } Normalize_mR(V); 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(V); } /* ------------------------------------ */ /* ------------------------------------ */ double **svd_U_Rn_mR( double **A, double **U ) { int i; 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(r,C1); double **EV = i_mR(c,r); double **Ab; 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); eye_mR(Ide); mul_mR(tA,A,tAA); for(i=R1; i<=rsize_R(EValue); i++) { 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,U,i); } Normalize_mR(U); f_mR(EValue); f_mR(EV); f_mR(tA); f_mR(tAA); f_mR(AtA); f_mR(Ab); f_mR(b); f_mR(Ide); f_mR(sIde); f_mR(tAAsIde); return(U); } /* ------------------------------------ */ /* ------------------------------------ */ double **Pinv_Rn_mR( double **A, double **Pinv, double FACTOR ) { int r = rsize_R(A); int c = csize_R(A); double **FA = smul_mR(FACTOR,A, i_mR(r,c)); double **tFA = transpose_mR(FA, i_mR(c,r)); double **U = i_mR(r,c); double **tU = i_mR(c,r); double **V = i_mR(c,c); double **tUA = i_mR(c,c); double **EValue = i_mR(c,c); double **invEValue = i_mR(c,c); double **VinvEValue = i_mR(c,c); svd_U_Rn_mR(tFA,U); svd_V_Rn_mR(tFA,V); /* EValue = tU A V */ transpose_mR(U,tU); mul_mR(tU,A,tUA); mul_mR(tUA,V,EValue); inv_svd_mR(EValue,invEValue); /* PseudoInverse = V invEValue tU */ mul_mR(V,invEValue,VinvEValue); mul_mR(VinvEValue,tU,Pinv); f_mR(U); f_mR(tU); f_mR(V); f_mR(FA); f_mR(tFA); f_mR(tUA); f_mR(EValue); f_mR(invEValue); f_mR(VinvEValue); return(Pinv); } /* ------------------------------------ */ /* ------------------------------------ */