Mathc matrices/h20b
Installer ce fichier dans votre répertoire de travail.
dot_diag.h |
---|
/* ------------------------------------ */
/* Save as : dot_diag.h */
/* ------------------------------------ */
/* ------------------------------------ */
/* ------------------------------------ */
double **f_sqrtdiag_mR(
double **Diag,
double **SqrtDiag
)
{
int r;
dif_sizes_mR(Diag,SqrtDiag,"f_sqrtdiag_mR;","Diag SqrtDiag");
for ( r=R1; r<Diag[R_SIZE][C0]; r++)
SqrtDiag[r][r] = sqrt(Diag[r][r]);
return(SqrtDiag);
}
/* ------------------------------------ */
/* ------------------------------------ */
/* ------- <u.v> = v^t SqrtA_T SqrtA u */
/* */
/* SqrtA_T == SqrtA */
/* ------------------------------------ */
double dot_Diag_R(
double **A,
double **u,
double **v
)
{
double **SqrtA = i_mR(csize_R(A), rsize_R(A));
double **vt = i_mR(R1, rsize_R(v));
double **T_vtAt = i_mR(R1, rsize_R(v));
double **T_vtAtA = i_mR(R1, rsize_R(v));
double **T_vtAtAu = i_mR(R1, C1);
double dotDiag;
transpose_mR(v,vt);
f_sqrtdiag_mR(A,SqrtA);
mul_mR(vt,SqrtA,T_vtAt);
mul_mR(T_vtAt,SqrtA,T_vtAtA);
mul_mR(T_vtAtA,u,T_vtAtAu);
dotDiag = T_vtAtAu[R1][C1];
f_mR(vt);
f_mR(T_vtAt);
f_mR(T_vtAtA);
f_mR(T_vtAtAu);
return (dotDiag);
}
/* ------------------------------------ */
/* ------ ||u|| = <u.u>^(1/2) ------- */
/* ------------------------------------ */
double norm_Diag_R(
double **A,
double **u
)
{
return ( pow( dot_Diag_R(A,u,u),(1./2.)) );
}
/* ------------------------------------ */
/* ------ d(u,v) = ||u-v||^(1/2) ------- */
/* ------------------------------------ */
double dist_Diag_R(
double **A,
double **u,
double **v
)
{
double **umnsv = i_mR(rsize_R(v),csize_R(v));
double d;
sub_mR(u,v,umnsv);
/* d(u,v) = ||u-v||^(1/2) */
d = norm_Diag_R(A,umnsv);
f_mR(umnsv);
return (d);
}
/* ------------------------------------ */
/* ------- <u.v> = v^t SqrtA_T SqrtA u */
/* ------------------------------------ */
double **proj_Diag_mR(
double **A,
double **u,
double **v,
double **projuv
)
{
/* <u,v> / ||v||^2 */
double s = dot_Diag_R(A,u,v) / dot_Diag_R(A,v,v);
/* (<u,v>/||v||^2) v */
smul_mR(s,v,projuv);
return(projuv);
}
/* ------------------------------------ */
/* ------- <u.v> = v^t SqrtA_T SqrtA u */
/* ------------------------------------ */
double **orth_Diag_mR(
double **A,
double **U,
double **Orth
)
{
double **Vect_U = i_mR(rsize_R(U),C1);
double **Vect_T = i_mR(rsize_R(U),C1);
double **Vect_Orth = i_mR(rsize_R(U),C1);
double **Vect_projuv = i_mR(rsize_R(U),C1);
int c_U;
int c_Orth;
c_c_mR(U,C1, Orth,C1);
for( c_U = C2; c_U < U[C_SIZE][C0]; c_U++)
{
c_c_mR(U,c_U, Vect_U,C1);
c_c_mR(U,c_U, Vect_T,C1);
for( c_Orth = C1; c_Orth < c_U; c_Orth++)
{
c_c_mR(Orth,c_Orth, Vect_Orth,C1);
proj_Diag_mR(A,Vect_U,Vect_Orth,Vect_projuv);
sub_mR(Vect_T,Vect_projuv,Vect_Orth);
c_mR(Vect_Orth,Vect_T);
}
c_c_mR(Vect_T,C1, Orth,c_Orth);
}
f_mR(Vect_U);
f_mR(Vect_T);
f_mR(Vect_Orth);
f_mR(Vect_projuv);
return(Orth);
}
/* ------------------------------------ */
/* ------- <u.v> = v^t SqrtA_T SqrtA u */
/* ------------------------------------ */
double **Normalize_Diag_mR(
double **A,
double **Orth
)
{
double **u = i_mR(rsize_R(Orth),C1);
double **Nu = i_mR(rsize_R(Orth),C1);
int c_Orth = csize_R(Orth);
int c;
for( c = C1; c <= c_Orth; c++)
{
c_c_mR(Orth,c, u,C1);
smul_mR(1./norm_Diag_R(A,u),u,Nu);
c_c_mR(Nu,C1, Orth,c);
}
f_mR(u);
f_mR(Nu);
return(Orth);
}
/* ------------------------------------ */
/* ------- <u.v> = v^t SqrtA_T SqrtA u */
/* ------------------------------------ */
double **mul_Diag_mR(
double **A,
double **U,
double **V,
double **UV
)
{
double **u = i_mR(rsize_R(U),C1);
double **v = i_mR(rsize_R(U),C1);
int c1;
int c2;
for ( c1=C1; c1<U[C_SIZE][C0]; c1++)
for ( c2=C1; c2<U[C_SIZE][C0]; c2++)
{
c_c_mR(U,c1, u,C1);
c_c_mR(V,c2, v,C1);
UV[c1][c2] = dot_Diag_R(A,u,v);
}
f_mR(u);
f_mR(v);
return(UV);
}
/* ------------------------------------ */
/* ------- <u.v> = v^t SqrtA_T SqrtA u */
/* ------------------------------------ */
void QR_Diag_mR(
double **A,
double **U,
double **Q,
double **R)
{
orth_Diag_mR(A,U,Q);
Normalize_Diag_mR(A,Q);
mul_Diag_mR(A,Q,U,R);
}
/* ------------------------------------ */
/* ------------------------------------ */
/* ------------------------------------ */
double **eigs_value_Diag_mR(
double **A,
double **U,
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(U,T);
for(i=0;i<LOOP_EIGSUV;i++)
{
QR_Diag_mR(A,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 **xeigs_value_Diag_mR(
double **A,
double **U,
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(U,T);
do
{
QR_Diag_mR(A,T,Q,R);
mul_mR(R,Q,T);
clrscrn();
pE_mR(T,S5,P3,C6);
} while(stop_w());
getchar();
for(rc=R1;rc<=r;rc++)
{
EigsValue[rc][C1] = T[rc][rc];
}
f_mR(T);
f_mR(Q);
f_mR(R);
return(EigsValue);
}*/
/* ------------------------------------ */
/* ------------------------------------ */
Déclaration des fichiers h.