Linear Algebra and the C Language/a02x

__NOTOC__ Install this file in your working directory. /* ------------------------------------ */ /* Save as : vdsvdcn.h */ /* ------------------------------------ */ /* ------------------------------------ */ 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); } /* ------------------------------------ */ /* ------------------------------------ */ 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],Ide,sIde); EValue[i][C1] = 0 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(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 **Pinv_Cn_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 **U = i_mR(r,r); double **tU = i_mR(r,r); double **V = i_mR(c,c); double **tUA = i_mR(r,c); double **EValue = i_mR(r,c); double **tEValue = i_mR(c,r); double **invtEValue = i_mR(c,r); double **VinvtEValue = i_mR(c,r); svd_U_Cn_mR(FA,U); svd_V_Cn_mR(FA,V); /* EValue = tU A V */ transpose_mR(U,tU); mul_mR(tU,A,tUA); mul_mR(tUA,V,EValue); transpose_mR(EValue,tEValue); inv_svd_mR(tEValue,invtEValue); /* PseudoInverse = V invtEValue tU */ mul_mR(V,invtEValue,VinvtEValue); mul_mR(VinvtEValue,tU,Pinv); f_mR(FA); f_mR(U); f_mR(V); f_mR(tU); f_mR(tUA); f_mR(EValue); f_mR(tEValue); f_mR(invtEValue); f_mR(VinvtEValue); return(Pinv); } /* ------------------------------------ */ /* ------------------------------------ */