ztorture.c

C语言版本的矩阵库

	MACHEPS*zm_norm_inf(A)*zm_norm_inf(B)*5 )
	errmesg("zget_row()/zget_col()");
    zv_mlt(zmake(-1.0,0.0),x,x);
    zv_mlt(zmake(-1.0,0.0),y,y);
    zset_row(A,3,x);
    zset_col(B,3,y);
    zm_mlt(A,B,C);
    for ( i = 0; i < C->m; i++ )
	zm_set_val(C,i,i,zsub(zm_entry(C,i,i),ONE));
    if ( zm_norm_inf(C) >= MACHEPS*zm_norm_inf(A)*zm_norm_inf(B)*5 )
	errmesg("zset_row()/zset_col()");

    MEMCHK();

    /* matrix norms */
    notice("matrix norms");
    A = zm_resize(A,11,15);
    zm_rand(A);
    s1 = zm_norm_inf(A);
    B = zm_adjoint(A,B);
    s2 = zm_norm1(B);
    if ( fabs(s1 - s2) >= MACHEPS*A->m )
	errmesg("zm_norm1()/zm_norm_inf()");
    C = zmam_mlt(A,A,C);
    z1.re = z1.im = 0.0;
    for ( i = 0; i < C->m && i < C->n; i++ )
	z1 = zadd(z1,zm_entry(C,i,i));
    if ( fabs(sqrt(z1.re) - zm_norm_frob(A)) >= MACHEPS*A->m*A->n )
	errmesg("zm_norm_frob");

    MEMCHK();
    
    /* permuting rows and columns */
    /******************************
    notice("permuting rows & cols");
    A = zm_resize(A,11,15);
    B = zm_resize(B,11,15);
    pi1 = px_resize(pi1,A->m);
    px_rand(pi1);
    x = zv_resize(x,A->n);
    y = zmv_mlt(A,x,y);
    px_rows(pi1,A,B);
    px_zvec(pi1,y,z);
    zmv_mlt(B,x,u);
    if ( zv_norm2(zv_sub(z,u,u)) >= MACHEPS*A->m )
	errmesg("px_rows()");
    pi1 = px_resize(pi1,A->n);
    px_rand(pi1);
    px_cols(pi1,A,B);
    pxinv_zvec(pi1,x,z);
    zmv_mlt(B,z,u);
    if ( zv_norm2(zv_sub(y,u,u)) >= MACHEPS*A->n )
	errmesg("px_cols()");
    ******************************/

    MEMCHK();

    /* MATLAB save/load */
    notice("MATLAB save/load");
    A = zm_resize(A,12,11);
    if ( (fp=fopen(SAVE_FILE,"w")) == (FILE *)NULL )
	printf("Cannot perform MATLAB save/load test\n");
    else
    {
	zm_rand(A);
	zm_save(fp, A, name);
	fclose(fp);
	if ( (fp=fopen(SAVE_FILE,"r")) == (FILE *)NULL )
	    printf("Cannot open save file \"%s\"\n",SAVE_FILE);
	else
	{
	    ZM_FREE(B);
	    B = zm_load(fp,&cp);
	    if ( strcmp(name,cp) || zm_norm1(zm_sub(A,B,C)) >=
		 MACHEPS*A->m )
	    {
		errmesg("zm_load()/zm_save()");
		printf("# orig. name = %s, restored name = %s\n", name, cp);
		printf("# orig. A =\n");	zm_output(A);
		printf("# restored A =\n");	zm_output(B);
	    }
	}
    }

    MEMCHK();

    /* Now, onto matrix factorisations */
    A = zm_resize(A,10,10);
    B = zm_resize(B,A->m,A->n);
    zm_copy(A,B);
    x = zv_resize(x,A->n);
    y = zv_resize(y,A->m);
    z = zv_resize(z,A->n);
    u = zv_resize(u,A->m);
    zv_rand(x);
    zmv_mlt(B,x,y);
    z = zv_copy(x,z);

    notice("LU factor/solve");
    pivot = px_get(A->m);
    zLUfactor(A,pivot);
    tracecatch(zLUsolve(A,pivot,y,x),"main");
    tracecatch(cond_est = zLUcondest(A,pivot),"main");
    printf("# cond(A) approx= %g\n", cond_est);
    if ( zv_norm2(zv_sub(x,z,u)) >= MACHEPS*zv_norm2(x)*cond_est)
    {
	errmesg("zLUfactor()/zLUsolve()");
	printf("# LU solution error = %g [cf MACHEPS = %g]\n",
	       zv_norm2(zv_sub(x,z,u)), MACHEPS);
    }


    zv_copy(y,x);
    tracecatch(zLUsolve(A,pivot,x,x),"main");
    tracecatch(cond_est = zLUcondest(A,pivot),"main");
    if ( zv_norm2(zv_sub(x,z,u)) >= MACHEPS*zv_norm2(x)*cond_est)
    {
	errmesg("zLUfactor()/zLUsolve()");
	printf("# LU solution error = %g [cf MACHEPS = %g]\n",
	       zv_norm2(zv_sub(x,z,u)), MACHEPS);
    }

    zvm_mlt(B,z,y);
    zv_copy(y,x);
    tracecatch(zLUAsolve(A,pivot,x,x),"main");
    if ( zv_norm2(zv_sub(x,z,u)) >= MACHEPS*zv_norm2(x)*cond_est)
    {
	errmesg("zLUfactor()/zLUAsolve()");
	printf("# LU solution error = %g [cf MACHEPS = %g]\n",
	       zv_norm2(zv_sub(x,z,u)), MACHEPS);
    }

    MEMCHK();

    /* QR factorisation */
    zm_copy(B,A);
    zmv_mlt(B,z,y);
    notice("QR factor/solve:");
    diag = zv_get(A->m);
    zQRfactor(A,diag);
    zQRsolve(A,diag,y,x);
    if ( zv_norm2(zv_sub(x,z,u)) >= MACHEPS*zv_norm2(x)*cond_est )
    {
	errmesg("zQRfactor()/zQRsolve()");
	printf("# QR solution error = %g [cf MACHEPS = %g]\n",
	       zv_norm2(zv_sub(x,z,u)), MACHEPS);
    }
    printf("# QR cond(A) approx= %g\n", zQRcondest(A));
    Q = zm_get(A->m,A->m);
    zmakeQ(A,diag,Q);
    zmakeR(A,A);
    zm_mlt(Q,A,C);
    zm_sub(B,C,C);
    if ( zm_norm1(C) >= MACHEPS*zm_norm1(Q)*zm_norm1(B) )
    {
	errmesg("zQRfactor()/zmakeQ()/zmakeR()");
	printf("# QR reconstruction error = %g [cf MACHEPS = %g]\n",
	       zm_norm1(C), MACHEPS);
    }

    MEMCHK();

    /* now try with a non-square matrix */
    A = zm_resize(A,15,7);
    zm_rand(A);
    B = zm_copy(A,B);
    diag = zv_resize(diag,A->n);
    x = zv_resize(x,A->n);
    y = zv_resize(y,A->m);
    zv_rand(y);
    zQRfactor(A,diag);
    x = zQRsolve(A,diag,y,x);
    /* z is the residual vector */
    zmv_mlt(B,x,z);	zv_sub(z,y,z);
    /* check B*.z = 0 */
    zvm_mlt(B,z,u);
    if ( zv_norm2(u) >= 100*MACHEPS*zm_norm1(B)*zv_norm2(y) )
    {
	errmesg("zQRfactor()/zQRsolve()");
	printf("# QR solution error = %g [cf MACHEPS = %g]\n",
	       zv_norm2(u), MACHEPS);
    }
    Q = zm_resize(Q,A->m,A->m);
    zmakeQ(A,diag,Q);
    zmakeR(A,A);
    zm_mlt(Q,A,C);
    zm_sub(B,C,C);
    if ( zm_norm1(C) >= MACHEPS*zm_norm1(Q)*zm_norm1(B) )
    {
	errmesg("zQRfactor()/zmakeQ()/zmakeR()");
	printf("# QR reconstruction error = %g [cf MACHEPS = %g]\n",
	       zm_norm1(C), MACHEPS);
    }
    D = zm_get(A->m,Q->m);
    zmam_mlt(Q,Q,D);
    for ( i = 0; i < D->m; i++ )
	zm_set_val(D,i,i,zsub(zm_entry(D,i,i),ONE));
    if ( zm_norm1(D) >= MACHEPS*zm_norm1(Q)*zm_norm_inf(Q) )
    {
	errmesg("QRfactor()/makeQ()/makeR()");
	printf("# QR orthogonality error = %g [cf MACHEPS = %g]\n",
	       zm_norm1(D), MACHEPS);
    }

    MEMCHK();

    /* QRCP factorisation */
    zm_copy(B,A);
    notice("QR factor/solve with column pivoting");
    pivot = px_resize(pivot,A->n);
    zQRCPfactor(A,diag,pivot);
    z = zv_resize(z,A->n);
    zQRCPsolve(A,diag,pivot,y,z);
    /* pxinv_zvec(pivot,z,x); */
    /* now compute residual (z) vector */
    zmv_mlt(B,x,z);	zv_sub(z,y,z);
    /* check B^T.z = 0 */
    zvm_mlt(B,z,u);
    if ( zv_norm2(u) >= MACHEPS*zm_norm1(B)*zv_norm2(y) )
    {
	errmesg("QRCPfactor()/QRsolve()");
	printf("# QR solution error = %g [cf MACHEPS = %g]\n",
	       zv_norm2(u), MACHEPS);
    }

    Q = zm_resize(Q,A->m,A->m);
    zmakeQ(A,diag,Q);
    zmakeR(A,A);
    zm_mlt(Q,A,C);
    ZM_FREE(D);
    D = zm_get(B->m,B->n);
    /******************************
    px_cols(pivot,C,D);
    zm_sub(B,D,D);
    if ( zm_norm1(D) >= MACHEPS*zm_norm1(Q)*zm_norm1(B) )
    {
	errmesg("QRCPfactor()/makeQ()/makeR()");
	printf("# QR reconstruction error = %g [cf MACHEPS = %g]\n",
	       zm_norm1(D), MACHEPS);
    }
    ******************************/

    /* Now check eigenvalue/SVD routines */
    notice("complex Schur routines");
    A = zm_resize(A,11,11);
    B = zm_resize(B,A->m,A->n);
    C = zm_resize(C,A->m,A->n);
    D = zm_resize(D,A->m,A->n);
    Q = zm_resize(Q,A->m,A->n);

    MEMCHK();

    /* now test complex Schur decomposition */
    /* zm_copy(A,B); */
    ZM_FREE(A);
    A = zm_get(11,11);
    zm_rand(A);
    B = zm_copy(A,B);
    MEMCHK();

    B = zschur(B,Q);
    checkpt();

    zm_mlt(Q,B,C);
    zmma_mlt(C,Q,D);
    MEMCHK();
    zm_sub(A,D,D);
    if ( zm_norm1(D) >= MACHEPS*zm_norm1(Q)*zm_norm_inf(Q)*zm_norm1(B)*5 )
    {
	errmesg("zschur()");
	printf("# Schur reconstruction error = %g [cf MACHEPS = %g]\n",
	       zm_norm1(D), MACHEPS);
    }

    /* orthogonality check */
    zmma_mlt(Q,Q,D);
    for ( i = 0; i < D->m; i++ )
	zm_set_val(D,i,i,zsub(zm_entry(D,i,i),ONE));
    if ( zm_norm1(D) >= MACHEPS*zm_norm1(Q)*zm_norm_inf(Q)*10 )
    {
	errmesg("zschur()");
	printf("# Schur orthogonality error = %g [cf MACHEPS = %g]\n",
	       zm_norm1(D), MACHEPS);
    }

    MEMCHK();

    /* now test SVD */
    /******************************
    A = zm_resize(A,11,7);
    zm_rand(A);
    U = zm_get(A->n,A->n);
    Q = zm_resize(Q,A->m,A->m);
    u = zv_resize(u,max(A->m,A->n));
    svd(A,Q,U,u);
    ******************************/
    /* check reconstruction of A */
    /******************************
    D = zm_resize(D,A->m,A->n);
    C = zm_resize(C,A->m,A->n);
    zm_zero(D);
    for ( i = 0; i < min(A->m,A->n); i++ )
	zm_set_val(D,i,i,v_entry(u,i));
    zmam_mlt(Q,D,C);
    zm_mlt(C,U,D);
    zm_sub(A,D,D);
    if ( zm_norm1(D) >= MACHEPS*zm_norm1(U)*zm_norm_inf(Q)*zm_norm1(A) )
    {
	errmesg("svd()");
	printf("# SVD reconstruction error = %g [cf MACHEPS = %g]\n",
	       zm_norm1(D), MACHEPS);
    }
    ******************************/
    /* check orthogonality of Q and U */
    /******************************
    D = zm_resize(D,Q->n,Q->n);
    zmam_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 ( zm_norm1(D) >= MACHEPS*zm_norm1(Q)*zm_norm_inf(Q)*5 )
    {
	errmesg("svd()");
	printf("# SVD orthognality error (Q) = %g [cf MACHEPS = %g\n",
	       zm_norm1(D), MACHEPS);
    }
    D = zm_resize(D,U->n,U->n);
    zmam_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 ( zm_norm1(D) >= MACHEPS*zm_norm1(U)*zm_norm_inf(U)*5 )
    {
	errmesg("svd()");
	printf("# SVD orthognality error (U) = %g [cf MACHEPS = %g\n",
	       zm_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");
    }
    ******************************/

    ZV_FREE(x);	ZV_FREE(y);	ZV_FREE(z);
    ZV_FREE(u);	ZV_FREE(diag);
    PX_FREE(pi1);	PX_FREE(pi2);	PX_FREE(pivot);
    ZM_FREE(A);	ZM_FREE(B);	ZM_FREE(C);
    ZM_FREE(D);	ZM_FREE(Q);

    mem_stat_free(1);

    MEMCHK();
    printf("# Finished torture test for complex numbers/vectors/matrices\n");
    mem_info();
}