torture.c
来自「C语言版本的矩阵库」· C语言 代码 · 共 1,053 行 · 第 1/2 页
C
1,053 行
for ( i = 0; i < D->m; i++ )
m_set_val(D,i,i,m_entry(D,i,i)-1.0);
if ( m_norm1(D) >= MACHEPS*m_norm1(Q)*m_norm_inf(Q) )
{
errmesg("QRfactor()/makeQ()/makeR()");
printf("# QR orthogonality error = %g [cf MACHEPS = %g]\n",
m_norm1(D), MACHEPS);
}
MEMCHK();
/* QRCP factorisation */
m_copy(B,A);
notice("QR factor/solve with column pivoting");
pivot = px_resize(pivot,A->n);
QRCPfactor(A,diag,pivot);
z = v_resize(z,A->n);
QRCPsolve(A,diag,pivot,y,z);
/* pxinv_vec(pivot,z,x); */
/* now compute residual (z) vector */
mv_mlt(B,x,z); v_sub(z,y,z);
/* check B^T.z = 0 */
vm_mlt(B,z,u);
if ( v_norm2(u) >= MACHEPS*m_norm1(B)*v_norm2(y) )
{
errmesg("QRCPfactor()/QRsolve()");
printf("# QR solution error = %g [cf MACHEPS = %g]\n",
v_norm2(u), MACHEPS);
}
Q = m_resize(Q,A->m,A->m);
makeQ(A,diag,Q);
makeR(A,A);
m_mlt(Q,A,C);
M_FREE(D);
D = m_get(B->m,B->n);
px_cols(pivot,C,D);
m_sub(B,D,D);
if ( m_norm1(D) >= MACHEPS*m_norm1(Q)*m_norm1(B) )
{
errmesg("QRCPfactor()/makeQ()/makeR()");
printf("# QR reconstruction error = %g [cf MACHEPS = %g]\n",
m_norm1(D), MACHEPS);
}
MEMCHK();
/* Cholesky and LDL^T factorisation */
/* Use these for normal equations approach */
notice("Cholesky factor/solve");
mtrm_mlt(B,B,A);
CHfactor(A);
u = v_resize(u,B->n);
vm_mlt(B,y,u);
z = v_resize(z,B->n);
CHsolve(A,u,z);
v_sub(x,z,z);
if ( v_norm2(z) >= MACHEPS*v_norm2(x)*100 )
{
errmesg("CHfactor()/CHsolve()");
printf("# Cholesky solution error = %g [cf MACHEPS = %g]\n",
v_norm2(z), MACHEPS);
}
/* modified Cholesky factorisation should be identical with Cholesky
factorisation provided the matrix is "sufficiently positive definite" */
mtrm_mlt(B,B,C);
MCHfactor(C,MACHEPS);
m_sub(A,C,C);
if ( m_norm1(C) >= MACHEPS*m_norm1(A) )
{
errmesg("MCHfactor()");
printf("# Modified Cholesky error = %g [cf MACHEPS = %g]\n",
m_norm1(C), MACHEPS);
}
/* now test the LDL^T factorisation -- using a negative def. matrix */
mtrm_mlt(B,B,A);
sm_mlt(-1.0,A,A);
m_copy(A,C);
LDLfactor(A);
LDLsolve(A,u,z);
w = v_get(A->m);
mv_mlt(C,z,w);
v_sub(w,u,w);
if ( v_norm2(w) >= MACHEPS*v_norm2(u)*m_norm1(C) )
{
errmesg("LDLfactor()/LDLsolve()");
printf("# LDL^T residual = %g [cf MACHEPS = %g]\n",
v_norm2(w), MACHEPS);
}
v_add(x,z,z);
if ( v_norm2(z) >= MACHEPS*v_norm2(x)*100 )
{
errmesg("LDLfactor()/LDLsolve()");
printf("# LDL^T solution error = %g [cf MACHEPS = %g]\n",
v_norm2(z), MACHEPS);
}
MEMCHK();
/* and now the Bunch-Kaufman-Parlett method */
/* set up D to be an indefinite diagonal matrix */
notice("Bunch-Kaufman-Parlett factor/solve");
D = m_resize(D,B->m,B->m);
m_zero(D);
w = v_resize(w,B->m);
v_rand(w);
for ( i = 0; i < w->dim; i++ )
if ( v_entry(w,i) >= 0.5 )
m_set_val(D,i,i,1.0);
else
m_set_val(D,i,i,-1.0);
/* set A <- B^T.D.B */
C = m_resize(C,B->n,B->n);
C = mtrm_mlt(B,D,C);
A = m_mlt(C,B,A);
C = m_resize(C,B->n,B->n);
C = m_copy(A,C);
/* ... and use BKPfactor() */
blocks = px_get(A->m);
pivot = px_resize(pivot,A->m);
x = v_resize(x,A->m);
y = v_resize(y,A->m);
z = v_resize(z,A->m);
v_rand(x);
mv_mlt(A,x,y);
BKPfactor(A,pivot,blocks);
printf("# BKP pivot =\n"); px_output(pivot);
printf("# BKP blocks =\n"); px_output(blocks);
BKPsolve(A,pivot,blocks,y,z);
/* compute & check residual */
mv_mlt(C,z,w);
v_sub(w,y,w);
if ( v_norm2(w) >= MACHEPS*m_norm1(C)*v_norm2(z) )
{
errmesg("BKPfactor()/BKPsolve()");
printf("# BKP residual size = %g [cf MACHEPS = %g]\n",
v_norm2(w), MACHEPS);
}
/* check update routines */
/* check LDLupdate() first */
notice("update L.D.L^T routine");
A = mtrm_mlt(B,B,A);
m_resize(C,A->m,A->n);
C = m_copy(A,C);
LDLfactor(A);
s1 = 3.7;
w = v_resize(w,A->m);
v_rand(w);
for ( i = 0; i < C->m; i++ )
for ( j = 0; j < C->n; j++ )
m_set_val(C,i,j,m_entry(C,i,j)+s1*v_entry(w,i)*v_entry(w,j));
LDLfactor(C);
LDLupdate(A,w,s1);
/* zero out strictly upper triangular parts of A and C */
for ( i = 0; i < A->m; i++ )
for ( j = i+1; j < A->n; j++ )
{
m_set_val(A,i,j,0.0);
m_set_val(C,i,j,0.0);
}
if ( m_norm1(m_sub(A,C,C)) >= sqrt(MACHEPS)*m_norm1(A) )
{
errmesg("LDLupdate()");
printf("# LDL update matrix error = %g [cf MACHEPS = %g]\n",
m_norm1(C), MACHEPS);
}
/* BAND MATRICES */
#define COL 40
#define UDIAG 5
#define LDIAG 2
smrand(101);
bA = bd_get(LDIAG,UDIAG,COL);
bB = bd_get(LDIAG,UDIAG,COL);
bC = bd_get(LDIAG,UDIAG,COL);
A = m_resize(A,COL,COL);
B = m_resize(B,COL,COL);
pivot = px_resize(pivot,COL);
x = v_resize(x,COL);
w = v_resize(w,COL);
z = v_resize(z,COL);
m_rand(A);
/* generate band matrix */
mat2band(A,LDIAG,UDIAG,bA);
band2mat(bA,A); /* now A is banded */
bB = bd_copy(bA,bB);
v_rand(x);
mv_mlt(A,x,w);
/* test of bd_mv_mlt */
notice("bd_mv_mlt");
bd_mv_mlt(bA,x,z);
v_sub(z,w,z);
if (v_norm2(z) > v_norm2(x)*sqrt(MACHEPS)) {
errmesg("incorrect vector (bd_mv_mlt)");
printf(" ||exact vector. - computed vector.|| = %g [MACHEPS = %g]\n",
v_norm2(z),MACHEPS);
}
z = v_copy(w,z);
notice("band LU factorization");
bdLUfactor(bA,pivot);
/* pivot will be changed */
bdLUsolve(bA,pivot,z,z);
v_sub(x,z,z);
if (v_norm2(z) > v_norm2(x)*sqrt(MACHEPS)) {
errmesg("incorrect solution (band LU factorization)");
printf(" ||exact sol. - computed sol.|| = %g [MACHEPS = %g]\n",
v_norm2(z),MACHEPS);
}
/* solve transpose system */
notice("band LU factorization for transpose system");
m_transp(A,B);
mv_mlt(B,x,w);
bd_copy(bB,bA);
bd_transp(bA,bA);
/* transposition in situ */
bd_transp(bA,bB);
bd_transp(bB,bB);
bdLUfactor(bB,pivot);
bdLUsolve(bB,pivot,w,z);
v_sub(x,z,z);
if (v_norm2(z) > v_norm2(x)*sqrt(MACHEPS)) {
errmesg("incorrect solution (band transposed LU factorization)");
printf(" ||exact sol. - computed sol.|| = %g [MACHEPS = %g]\n",
v_norm2(z),MACHEPS);
}
/* Cholesky factorization */
notice("band Choleski LDL' factorization");
m_add(A,B,A); /* symmetric matrix */
for (i=0; i < COL; i++) /* positive definite */
A->me[i][i] += 2*LDIAG;
mat2band(A,LDIAG,LDIAG,bA);
band2mat(bA,A); /* corresponding matrix A */
v_rand(x);
mv_mlt(A,x,w);
z = v_copy(w,z);
bdLDLfactor(bA);
z = bdLDLsolve(bA,z,z);
v_sub(x,z,z);
if (v_norm2(z) > v_norm2(x)*sqrt(MACHEPS)) {
errmesg("incorrect solution (band LDL' factorization)");
printf(" ||exact sol. - computed sol.|| = %g [MACHEPS = %g]\n",
v_norm2(z),MACHEPS);
}
/* new bandwidths */
m_rand(A);
bA = bd_resize(bA,UDIAG,LDIAG,COL);
bB = bd_resize(bB,UDIAG,LDIAG,COL);
mat2band(A,UDIAG,LDIAG,bA);
band2mat(bA,A);
bd_copy(bA,bB);
mv_mlt(A,x,w);
notice("band LU factorization (resized)");
bdLUfactor(bA,pivot);
/* pivot will be changed */
bdLUsolve(bA,pivot,w,z);
v_sub(x,z,z);
if (v_norm2(z) > v_norm2(x)*sqrt(MACHEPS)) {
errmesg("incorrect solution (band LU factorization)");
printf(" ||exact sol. - computed sol.|| = %g [MACHEPS = %g]\n",
v_norm2(z),MACHEPS);
}
/* testing transposition */
notice("band matrix transposition");
m_zero(bA->mat);
bd_copy(bB,bA);
m_zero(bB->mat);
bd_copy(bA,bB);
bd_transp(bB,bB);
bd_transp(bB,bB);
m_zero(bC->mat);
bd_copy(bB,bC);
m_sub(bA->mat,bC->mat,bC->mat);
if (m_norm_inf(bC->mat) > MACHEPS*bC->mat->n) {
errmesg("band transposition");
printf(" difference ||A - (A')'|| = %g\n",m_norm_inf(bC->mat));
}
bd_free(bA);
bd_free(bB);
bd_free(bC);
MEMCHK();
/* now check QRupdate() routine */
notice("update QR routine");
B = m_resize(B,15,7);
A = m_resize(A,B->m,B->n);
m_copy(B,A);
diag = v_resize(diag,A->n);
beta = v_resize(beta,A->n);
QRfactor(A,diag);
Q = m_resize(Q,A->m,A->m);
makeQ(A,diag,Q);
makeR(A,A);
m_resize(C,A->m,A->n);
w = v_resize(w,A->m);
v = v_resize(v,A->n);
u = v_resize(u,A->m);
v_rand(w);
v_rand(v);
vm_mlt(Q,w,u);
QRupdate(Q,A,u,v);
m_mlt(Q,A,C);
for ( i = 0; i < B->m; i++ )
for ( j = 0; j < B->n; j++ )
m_set_val(B,i,j,m_entry(B,i,j)+v_entry(w,i)*v_entry(v,j));
m_sub(B,C,C);
if ( m_norm1(C) >= MACHEPS*m_norm1(A)*m_norm1(Q)*2 )
{
errmesg("QRupdate()");
printf("# Reconstruction error in QR update = %g [cf MACHEPS = %g]\n",
m_norm1(C), MACHEPS);
}
m_resize(D,Q->m,Q->n);
mtrm_mlt(Q,Q,D);
for ( i = 0; i < D->m; i++ )
m_set_val(D,i,i,m_entry(D,i,i)-1.0);
if ( m_norm1(D) >= 10*MACHEPS*m_norm1(Q)*m_norm_inf(Q) )
{
errmesg("QRupdate()");
printf("# QR update orthogonality error = %g [cf MACHEPS = %g]\n",
m_norm1(D), MACHEPS);
}
/* Now check eigenvalue/SVD routines */
notice("eigenvalue and SVD routines");
A = m_resize(A,11,11);
B = m_resize(B,A->m,A->n);
C = m_resize(C,A->m,A->n);
D = m_resize(D,A->m,A->n);
Q = m_resize(Q,A->m,A->n);
m_rand(A);
/* A <- A + A^T for symmetric case */
m_add(A,m_transp(A,C),A);
u = v_resize(u,A->m);
u = symmeig(A,Q,u);
m_zero(B);
for ( i = 0; i < B->m; i++ )
m_set_val(B,i,i,v_entry(u,i));
m_mlt(Q,B,C);
mmtr_mlt(C,Q,D);
m_sub(A,D,D);
if ( m_norm1(D) >= MACHEPS*m_norm1(Q)*m_norm_inf(Q)*v_norm_inf(u)*3 )
{
errmesg("symmeig()");
printf("# Reconstruction error = %g [cf MACHEPS = %g]\n",
m_norm1(D), MACHEPS);
}
mtrm_mlt(Q,Q,D);
for ( i = 0; i < D->m; i++ )
m_set_val(D,i,i,m_entry(D,i,i)-1.0);
if ( m_norm1(D) >= MACHEPS*m_norm1(Q)*m_norm_inf(Q)*3 )
{
errmesg("symmeig()");
printf("# symmeig() orthogonality error = %g [cf MACHEPS = %g]\n",
m_norm1(D), MACHEPS);
}
MEMCHK();
/* now test (real) Schur decomposition */
/* m_copy(A,B); */
M_FREE(A);
A = m_get(11,11);
m_rand(A);
B = m_copy(A,B);
MEMCHK();
B = schur(B,Q);
MEMCHK();
m_mlt(Q,B,C);
mmtr_mlt(C,Q,D);
MEMCHK();
m_sub(A,D,D);
if ( m_norm1(D) >= MACHEPS*m_norm1(Q)*m_norm_inf(Q)*m_norm1(B)*5 )
{
errmesg("schur()");
printf("# Schur reconstruction error = %g [cf MACHEPS = %g]\n",
m_norm1(D), MACHEPS);
}
/* orthogonality check */
mmtr_mlt(Q,Q,D);
for ( i = 0; i < D->m; i++ )
m_set_val(D,i,i,m_entry(D,i,i)-1.0);
if ( m_norm1(D) >= MACHEPS*m_norm1(Q)*m_norm_inf(Q)*10 )
{
errmesg("schur()");
printf("# Schur orthogonality error = %g [cf MACHEPS = %g]\n",
m_norm1(D), MACHEPS);
}
MEMCHK();
/* now test SVD */
A = m_resize(A,11,7);
m_rand(A);
U = m_get(A->n,A->n);
Q = m_resize(Q,A->m,A->m);
u = v_resize(u,max(A->m,A->n));
svd(A,Q,U,u);
/* check reconstruction of A */
D = m_resize(D,A->m,A->n);
C = m_resize(C,A->m,A->n);
m_zero(D);
for ( i = 0; i < min(A->m,A->n); i++ )
m_set_val(D,i,i,v_entry(u,i));
mtrm_mlt(Q,D,C);
m_mlt(C,U,D);
m_sub(A,D,D);
if ( m_norm1(D) >= MACHEPS*m_norm1(U)*m_norm_inf(Q)*m_norm1(A) )
{
errmesg("svd()");
printf("# SVD reconstruction error = %g [cf MACHEPS = %g]\n",
m_norm1(D), MACHEPS);
}
/* check orthogonality of Q and U */
D = m_resize(D,Q->n,Q->n);
mtrm_mlt(Q,Q,D);
for ( i = 0; i < D->m; i++ )
m_set_val(D,i,i,m_entry(D,i,i)-1.0);
if ( m_norm1(D) >= MACHEPS*m_norm1(Q)*m_norm_inf(Q)*5 )
{
errmesg("svd()");
printf("# SVD orthognality error (Q) = %g [cf MACHEPS = %g\n",
m_norm1(D), MACHEPS);
}
D = m_resize(D,U->n,U->n);
mtrm_mlt(U,U,D);
for ( i = 0; i < D->m; i++ )
m_set_val(D,i,i,m_entry(D,i,i)-1.0);
if ( m_norm1(D) >= MACHEPS*m_norm1(U)*m_norm_inf(U)*5 )
{
errmesg("svd()");
printf("# SVD orthognality error (U) = %g [cf MACHEPS = %g\n",
m_norm1(D), MACHEPS);
}
for ( i = 0; i < u->dim; i++ )
if ( v_entry(u,i) < 0 || (i < u->dim-1 &&
v_entry(u,i+1) > v_entry(u,i)) )
break;
if ( i < u->dim )
{
errmesg("svd()");
printf("# SVD sorting error\n");
}
/* test of long vectors */
notice("Long vectors");
x = v_resize(x,100000);
y = v_resize(y,100000);
z = v_resize(z,100000);
v_rand(x);
v_rand(y);
v_mltadd(x,y,3.0,z);
sv_mlt(1.0/3.0,z,z);
v_mltadd(z,x,-1.0/3.0,z);
v_sub(z,y,x);
if (v_norm2(x) >= MACHEPS*(x->dim)) {
errmesg("long vectors");
printf(" norm = %g\n",v_norm2(x));
}
mem_stat_free(1);
MEMCHK();
/**************************************************
VEC *x, *y, *z, *u, *v, *w;
VEC *diag, *beta;
PERM *pi1, *pi2, *pi3, *pivot, *blocks;
MAT *A, *B, *C, *D, *Q, *U;
**************************************************/
V_FREE(x); V_FREE(y); V_FREE(z);
V_FREE(u); V_FREE(v); V_FREE(w);
V_FREE(diag); V_FREE(beta);
PX_FREE(pi1); PX_FREE(pi2); PX_FREE(pi3);
PX_FREE(pivot); PX_FREE(blocks);
M_FREE(A); M_FREE(B); M_FREE(C);
M_FREE(D); M_FREE(Q); M_FREE(U);
MEMCHK();
printf("# Finished torture test\n");
mem_info();
return 0;
}
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?