Mathc matrices/h05b
Installer ce fichier dans votre répertoire de travail.
|
vdotu.h |
|---|
/* ------------------------------------ */
/* 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);
}
/* ------------------------------------ */
/* ------------------------------------ */
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);
}
/* ------------------------------------ */
/* ------------------------------------ */
Déclaration des fichiers h.