Linear Algebra and the C Language/a02k

__NOTOC__ Install this file in your working directory. /* ------------------------------------ */ /* Save as : vdotu.h */ /* ------------------------------------ */ /* ------------------------------------ */ /* ------- u.v = v^t u ---------------- */ /* ------------------------------------ */ double dot_R( double **u, double **v ) { double **vt = i_mR(R1, rsize_R(v)); double **vtu = i_mR(R1, C1); double udotv; transpose_mR(v,vt); mul_mR(vt,u,vtu); udotv = vtu[R1][C1]; f_mR(vt); f_mR(vtu); return (udotv); } /* ------------------------------------ */ /* ------ ||u|| = (u.u)^(1/2) ------- */ /* ------------------------------------ */ double norm_R( double **u ) { return (pow(dot_R(u,u),(1./2.))); } /* ------------------------------------ */ /* -- d(u,v) = ||u-v|| -- */ /* ------------------------------------ */ double dist_R( double **u, double **v ) { double **umnsv = i_mR(rsize_R(v),csize_R(v)); double d; sub_mR(u,v,umnsv); d = norm_R(umnsv); f_mR(umnsv); return (d); } /* ------------------------------------ */ /* (<u,v>/||v||^2) * v */ /* ------------------------------------ */ double **proj_mR( double **u, double **v, double **projuv ) { /* <u,v> / ||v||^2 */ double s = dot_R(u,v) / dot_R(v,v); smul_mR(s,v,projuv); return(projuv); } /* ------------------------------------ */ /* ------------------------------------ */ double **normalize_mR( double **U, double **Q ) { double **u = i_mR(rsize_R(U),C1); double **q = i_mR(rsize_R(U),C1); int c_U = csize_R(U); int c; for( c = C1; c <= c_U; c++) { c_c_mR(U,c, u,C1); smul_mR(1./norm_R(u),u,q); c_c_mR(q,C1, Q,c); } f_mR(u); f_mR(q); return(Q); } /* ------------------------------------ */ /* ------------------------------------ */ double **Normalize_mR( double **Q ) { double **u = i_mR(rsize_R(Q),C1); double **q = i_mR(rsize_R(Q),C1); int c_Q = csize_R(Q); int c; for( c = C1; c <= c_Q; c++) { c_c_mR(Q,c, u,C1); smul_mR(1./norm_R(u),u, q); c_c_mR(q,C1, Q,c); } f_mR(u); f_mR(q); return(Q); } /* ------------------------------------ */ /* cos(alpha) = <p,q> / ||p|| ||q|| */ /* ------------------------------------ */ double cos_R( double **u, double **v ) { return(( dot_R(u,v) / (norm_R(u) * norm_R(v)) )); } /* ------------------------------------ */ /* ------------------------------------ */ double **orth_mR( double **U, double **Q, int Normal /* YES = 1 NO = 0 */ ) { double **u = i_mR(rsize_R(U),C1); double **v = i_mR(rsize_R(U),C1); double **q = i_mR(rsize_R(U),C1); double **projuv = i_mR(rsize_R(U),C1); double **tq = i_mR(rsize_R(U),C1); int c_U; int c_Q; c_c_mR(U,C1, Q,C1); for( c_U = C2; c_U < U[C_SIZE][C0]; c_U++) { c_c_mR(U,c_U, u,C1); c_c_mR(U,c_U, q,C1); for( c_Q = C1; c_Q < c_U; c_Q++) { c_c_mR(Q,c_Q, v,C1); proj_mR(u,v, projuv); /* [<u,v> / ||v||^2] * v */ sub_mR(q, projuv, tq); c_mR( tq, q); } c_c_mR(q,C1, Q,c_Q); } if(Normal) Normalize_mR(Q); f_mR(u); f_mR(v); f_mR(q); f_mR(projuv); f_mR(tq); return(Q); } /* ------------------------------------ */ /* ------------------------------------ */ /* ------------------------------------ */ /* ------------------------------------ double **xorth_mR( double **U, double **Q ) { double **u = i_mR(rsize_R(U),C1); double **v = i_mR(rsize_R(U),C1); double **q = i_mR(rsize_R(U),C1); double **projuv = i_mR(rsize_R(U),C1); int c_U; int c_Q; c_c_mR(U,C1, Q,C1); for( c_U = C2; c_U < U[C_SIZE][C0]; c_U++) { c_c_mR(U,c_U, u,C1); c_c_mR(U,c_U, q,C1); for( c_Q = C1; c_Q < c_U; c_Q++) { c_c_mR(Q,c_Q, v,C1); proj_mR(u,v, projuv); sub_mR(q,projuv, v); c_mR(v,q); } c_c_mR(q,C1, Q,c_Q); } Normalize_mR(Q); f_mR(u); f_mR(v); f_mR(q); f_mR(projuv); return(Q); }*/ /* ------------------------------------ */ /* ------------------------------------ */ void QR_mR( double **U, double **Q, double **R) { double **Q_T = i_mR(csize_R(Q),rsize_R(Q)); orth_mR(U,Q,YES); transpose_mR(Q,Q_T); mul_mR(Q_T,U, R); f_mR(Q_T); } /* ------------------------------------ */ double **eigs_mR( double **A, double **EigsValue ) { int r = rsize_R(A); double **T = i_mR(r,r); double **Q = i_mR(r,r); double **R = i_mR(r,r); int i; int rc; c_mR(A,T); for(i=0;i<LOOP_EIGSUV;i++) { QR_mR(T,Q,R); mul_mR(R,Q,T); } for(rc=R1;rc<=r;rc++) { EigsValue[rc][C1] = T[rc][rc]; } f_mR(T); f_mR(Q); f_mR(R); return(EigsValue); } /* ------------------------------------ */ double **sqrt_svd_mR( double **A, double **sqrtA ) { int r; int c; dif_sizes_mR(A,sqrtA,"svd_sqrt_mR();","(A or sqrtA)"); for ( r=R1; r<A[R_SIZE][C0]; r++) for ( c=C1; c<A[C_SIZE][C0]; c++) sqrtA[r][c] = sqrt(A[r][c]); return(sqrtA); } /* ------------------------------------ */ double **inv_svd_mR( double **S, double **invS ) { int r; int c; dif_sizes_mR(S,invS, "inv_svd_mR();","(S,invS)"); for ( r=R1; r<S[R_SIZE][C0]; r++) for ( c=C1; c<S[C_SIZE][C0]; c++) if(r==c) { if(fabs(S[r][c]) < ERROR_E) invS[r][c] = 0; else invS[r][c] = (1./S[r][c]); } else invS[r][c] = 0; return(invS); } /* ------------------------------------ */ /* ------------------------------------ */ /* ------------------------------------ */ double **svds_mR( double **A, double **SvdValue) { int rA = rsize_R(A); int cA = csize_R(A); double **A_T = i_mR(cA,rA); double **S = i_mR(cA,cA); double **SvdValue_2 = i_mR(cA,C1); transpose_mR(A,A_T); mul_mR(A_T,A, S); eigs_mR(S,SvdValue_2); sqrt_svd_mR(SvdValue_2,SvdValue); f_mR(A_T); f_mR(S); f_mR(SvdValue_2); return(SvdValue); } /* ------------------------------------ */ /* ------------------------------------ */