torture.c

来自「C语言版本的矩阵库」· C语言 代码 · 共 1,053 行 · 第 1/2 页

/**************************************************************************
**
** Copyright (C) 1993 David E. Stewart & Zbigniew Leyk, all rights reserved.
**
**			     Meschach Library
** 
** This Meschach Library is provided "as is" without any express 
** or implied warranty of any kind with respect to this software. 
** In particular the authors shall not be liable for any direct, 
** indirect, special, incidental or consequential damages arising 
** in any way from use of the software.
** 
** Everyone is granted permission to copy, modify and redistribute this
** Meschach Library, provided:
**  1.  All copies contain this copyright notice.
**  2.  All modified copies shall carry a notice stating who
**      made the last modification and the date of such modification.
**  3.  No charge is made for this software or works derived from it.  
**      This clause shall not be construed as constraining other software
**      distributed on the same medium as this software, nor is a
**      distribution fee considered a charge.
**
***************************************************************************/

/*
	This file contains a series of tests for the Meschach matrix
	library, parts 1 and 2
*/

static char rcsid[] = "$Id: torture.c,v 1.6 1994/08/25 15:22:11 des Exp $";

#include	<stdio.h>
#include	<math.h>
#include	"matrix2.h"
#include        "matlab.h"

#define	errmesg(mesg)	printf("Error: %s error: line %d\n",mesg,__LINE__)
#define notice(mesg)	printf("# Testing %s...\n",mesg);

static char *test_err_list[] = {
   "unknown error",			/* 0 */
   "testing error messages",		/* 1 */
   "unexpected end-of-file"		/* 2 */
};


#define MAX_TEST_ERR   (sizeof(test_err_list)/sizeof(char *))

/* extern	int	malloc_chain_check(); */
/* #define MEMCHK() if ( malloc_chain_check(0) ) \
{ printf("Error in malloc chain: \"%s\", line %d\n", \
	 __FILE__, __LINE__); exit(0); } */
#define	MEMCHK() 

/* cmp_perm -- returns 1 if pi1 == pi2, 0 otherwise */
int	cmp_perm(pi1, pi2)
PERM	*pi1, *pi2;
{
    int		i;

    if ( ! pi1 || ! pi2 )
	error(E_NULL,"cmp_perm");
    if ( pi1->size != pi2->size )
	return 0;
    for ( i = 0; i < pi1->size; i++ )
	if ( pi1->pe[i] != pi2->pe[i] )
	    return 0;
    return 1;
}

/* px_rand -- generates sort-of random permutation */
PERM	*px_rand(pi)
PERM	*pi;
{
    int		i, j, k;

    if ( ! pi )
	error(E_NULL,"px_rand");

    for ( i = 0; i < 3*pi->size; i++ )
    {
	j = (rand() >> 8) % pi->size;
	k = (rand() >> 8) % pi->size;
	px_transp(pi,j,k);
    }

    return pi;
}

#define	SAVE_FILE	"asx5213a.mat"
#define	MATLAB_NAME	"alpha"
char	name[81] = MATLAB_NAME;

int main(argc, argv)
int	argc;
char	*argv[];
{
   VEC	*x = VNULL, *y = VNULL, *z = VNULL, *u = VNULL, *v = VNULL, 
        *w = VNULL;
   VEC	*diag = VNULL, *beta = VNULL;
   PERM	*pi1 = PNULL, *pi2 = PNULL, *pi3 = PNULL, *pivot = PNULL, 
        *blocks = PNULL;
   MAT	*A = MNULL, *B = MNULL, *C = MNULL, *D = MNULL, *Q = MNULL, 
        *U = MNULL;
   BAND *bA, *bB, *bC;
   Real	cond_est, s1, s2, s3;
   int	i, j, seed;
   FILE	*fp;
   char	*cp;


    mem_info_on(TRUE);

    setbuf(stdout,(char *)NULL);

    seed = 1111;
    if ( argc > 2 )
    {
	printf("usage: %s [seed]\n",argv[0]);
	exit(0);
    }
    else if ( argc == 2 )
	sscanf(argv[1], "%d", &seed);

    /* set seed for rand() */
    smrand(seed);

    mem_stat_mark(1);

    /* print version information */
    m_version();

    printf("# grep \"^Error\" the output for a listing of errors\n");
    printf("# Don't panic if you see \"Error\" appearing; \n");
    printf("# Also check the reported size of error\n");
    printf("# This program uses randomly generated problems and therefore\n");
    printf("# may occasionally produce ill-conditioned problems\n");
    printf("# Therefore check the size of the error compared with MACHEPS\n");
    printf("# If the error is within 1000*MACHEPS then don't worry\n");
    printf("# If you get an error of size 0.1 or larger there is \n");
    printf("# probably a bug in the code or the compilation procedure\n\n");
    printf("# seed = %d\n",seed);

    printf("# Check: MACHEPS = %g\n",MACHEPS);
    /* allocate, initialise, copy and resize operations */
    /* VEC */
    notice("vector initialise, copy & resize");
    x = v_get(12);
    y = v_get(15);
    z = v_get(12);
    v_rand(x);
    v_rand(y);
    z = v_copy(x,z);
    if ( v_norm2(v_sub(x,z,z)) >= MACHEPS )
	errmesg("VEC copy");
    v_copy(x,y);
    x = v_resize(x,10);
    y = v_resize(y,10);
    if ( v_norm2(v_sub(x,y,z)) >= MACHEPS )
	errmesg("VEC copy/resize");
    x = v_resize(x,15);
    y = v_resize(y,15);
    if ( v_norm2(v_sub(x,y,z)) >= MACHEPS )
	errmesg("VEC resize");

    /* MAT */
    notice("matrix initialise, copy & resize");
    A = m_get(8,5);
    B = m_get(3,9);
    C = m_get(8,5);
    m_rand(A);
    m_rand(B);
    C = m_copy(A,C);
    if ( m_norm_inf(m_sub(A,C,C)) >= MACHEPS )
	errmesg("MAT copy");
    m_copy(A,B);
    A = m_resize(A,3,5);
    B = m_resize(B,3,5);
    if ( m_norm_inf(m_sub(A,B,C)) >= MACHEPS )
	errmesg("MAT copy/resize");
    A = m_resize(A,10,10);
    B = m_resize(B,10,10);
    if ( m_norm_inf(m_sub(A,B,C)) >= MACHEPS )
	errmesg("MAT resize");

    MEMCHK();

    /* PERM */
    notice("permutation initialise, inverting & permuting vectors");
    pi1 = px_get(15);
    pi2 = px_get(12);
    px_rand(pi1);
    v_rand(x);
    px_vec(pi1,x,z);
    y = v_resize(y,x->dim);
    pxinv_vec(pi1,z,y);
    if ( v_norm2(v_sub(x,y,z)) >= MACHEPS )
	errmesg("PERMute vector");
    pi2 = px_inv(pi1,pi2);
    pi3 = px_mlt(pi1,pi2,PNULL);
    for ( i = 0; i < pi3->size; i++ )
	if ( pi3->pe[i] != i )
	    errmesg("PERM inverse/multiply");

    /* testing catch() etc */
    notice("error handling routines");
    catch(E_NULL,
	  catchall(v_add(VNULL,VNULL,VNULL);
		     errmesg("tracecatch() failure"),
		     printf("# tracecatch() caught error\n");
		     error(E_NULL,"main"));
	             errmesg("catch() failure"),
	  printf("# catch() caught E_NULL error\n"));

    /* testing attaching a new error list (error list 2) */

    notice("attaching error lists");
    printf("# IT IS NOT A REAL WARNING ... \n");
    err_list_attach(2,MAX_TEST_ERR,test_err_list,TRUE);
    if (!err_is_list_attached(2)) 
       errmesg("attaching the error list 2");
    ev_err(__FILE__,1,__LINE__,"main",2);
    err_list_free(2);
    if (err_is_list_attached(2)) 
       errmesg("detaching the error list 2");

    /* testing inner products and v_mltadd() etc */
    notice("inner products and linear combinations");
    u = v_get(x->dim);
    v_rand(u);
    v_rand(x);
    v_resize(y,x->dim);
    v_rand(y);
    v_mltadd(y,x,-in_prod(x,y)/in_prod(x,x),z);
    if ( fabs(in_prod(x,z)) >= MACHEPS*x->dim )
	errmesg("v_mltadd()/in_prod()");
    s1 = -in_prod(x,y)/(v_norm2(x)*v_norm2(x));
    sv_mlt(s1,x,u);
    v_add(y,u,u);
    if ( v_norm2(v_sub(u,z,u)) >= MACHEPS*x->dim )
	errmesg("sv_mlt()/v_norm2()");

#ifdef ANSI_C 
    v_linlist(u,x,s1,y,1.0,VNULL);
    if ( v_norm2(v_sub(u,z,u)) >= MACHEPS*x->dim )
	errmesg("v_linlist()");
#endif
#ifdef VARARGS
    v_linlist(u,x,s1,y,1.0,VNULL);
    if ( v_norm2(v_sub(u,z,u)) >= MACHEPS*x->dim )
	errmesg("v_linlist()");
#endif


    MEMCHK();

    /* vector norms */
    notice("vector norms");
    x = v_resize(x,12);
    v_rand(x);
    for ( i = 0; i < x->dim; i++ )
	if ( v_entry(x,i) >= 0.5 )
	    v_set_val(x,i,1.0);
        else
	    v_set_val(x,i,-1.0);
    s1 = v_norm1(x);
    s2 = v_norm2(x);	
    s3 = v_norm_inf(x);
    if ( fabs(s1 - x->dim) >= MACHEPS*x->dim ||
	 fabs(s2 - sqrt((Real)(x->dim))) >= MACHEPS*x->dim ||
	 fabs(s3 - 1.0) >= MACHEPS )
	errmesg("v_norm1/2/_inf()");

    /* test matrix multiply etc */
    notice("matrix multiply and invert");
    A = m_resize(A,10,10);
    B = m_resize(B,10,10);
    m_rand(A);
    m_inverse(A,B);
    m_mlt(A,B,C);
    for ( i = 0; i < C->m; i++ )
	m_set_val(C,i,i,m_entry(C,i,i)-1.0);
    if ( m_norm_inf(C) >= MACHEPS*m_norm_inf(A)*m_norm_inf(B)*5 )
	errmesg("m_inverse()/m_mlt()");

    MEMCHK();

    /* ... and transposes */
    notice("transposes and transpose-multiplies");
    m_transp(A,A);	/* can do square matrices in situ */
    mtrm_mlt(A,B,C);
    for ( i = 0; i < C->m; i++ )
	m_set_val(C,i,i,m_entry(C,i,i)-1.0);
    if ( m_norm_inf(C) >= MACHEPS*m_norm_inf(A)*m_norm_inf(B)*5 )
	errmesg("m_transp()/mtrm_mlt()");
    m_transp(A,A);
    m_transp(B,B);
    mmtr_mlt(A,B,C);
    for ( i = 0; i < C->m; i++ )
	m_set_val(C,i,i,m_entry(C,i,i)-1.0);
    if ( m_norm_inf(C) >= MACHEPS*m_norm_inf(A)*m_norm_inf(B)*5 )
	errmesg("m_transp()/mmtr_mlt()");
    sm_mlt(3.71,B,B);
    mmtr_mlt(A,B,C);
    for ( i = 0; i < C->m; i++ )
	m_set_val(C,i,i,m_entry(C,i,i)-3.71);
    if ( m_norm_inf(C) >= MACHEPS*m_norm_inf(A)*m_norm_inf(B)*5 )
	errmesg("sm_mlt()/mmtr_mlt()");
    m_transp(B,B);
    sm_mlt(1.0/3.71,B,B);

    MEMCHK();

    /* ... and matrix-vector multiplies */
    notice("matrix-vector multiplies");
    x = v_resize(x,A->n);
    y = v_resize(y,A->m);
    z = v_resize(z,A->m);
    u = v_resize(u,A->n);
    v_rand(x);
    v_rand(y);
    mv_mlt(A,x,z);
    s1 = in_prod(y,z);
    vm_mlt(A,y,u);
    s2 = in_prod(u,x);
    if ( fabs(s1 - s2) >= (MACHEPS*x->dim)*x->dim )
	errmesg("mv_mlt()/vm_mlt()");
    mv_mlt(B,z,u);
    if ( v_norm2(v_sub(u,x,u)) >= MACHEPS*m_norm_inf(A)*m_norm_inf(B)*5 )
	errmesg("mv_mlt()/m_inverse()");

    MEMCHK();

    /* get/set row/col */
    notice("getting and setting rows and cols");
    x = v_resize(x,A->n);
    y = v_resize(y,B->m);
    x = get_row(A,3,x);
    y = get_col(B,3,y);
    if ( fabs(in_prod(x,y) - 1.0) >= MACHEPS*m_norm_inf(A)*m_norm_inf(B)*5 )
	errmesg("get_row()/get_col()");
    sv_mlt(-1.0,x,x);
    sv_mlt(-1.0,y,y);
    set_row(A,3,x);
    set_col(B,3,y);
    m_mlt(A,B,C);
    for ( i = 0; i < C->m; i++ )
	m_set_val(C,i,i,m_entry(C,i,i)-1.0);
    if ( m_norm_inf(C) >= MACHEPS*m_norm_inf(A)*m_norm_inf(B)*5 )
	errmesg("set_row()/set_col()");

    MEMCHK();

    /* matrix norms */
    notice("matrix norms");
    A = m_resize(A,11,15);
    m_rand(A);
    s1 = m_norm_inf(A);
    B = m_transp(A,B);
    s2 = m_norm1(B);
    if ( fabs(s1 - s2) >= MACHEPS*A->m )
	errmesg("m_norm1()/m_norm_inf()");
    C = mtrm_mlt(A,A,C);
    s1 = 0.0;
    for ( i = 0; i < C->m && i < C->n; i++ )
	s1 += m_entry(C,i,i);
    if ( fabs(sqrt(s1) - m_norm_frob(A)) >= MACHEPS*A->m*A->n )
	errmesg("m_norm_frob");

    MEMCHK();
    
    /* permuting rows and columns */
    notice("permuting rows & cols");
    A = m_resize(A,11,15);
    B = m_resize(B,11,15);
    pi1 = px_resize(pi1,A->m);
    px_rand(pi1);
    x = v_resize(x,A->n);
    y = mv_mlt(A,x,y);
    px_rows(pi1,A,B);
    px_vec(pi1,y,z);
    mv_mlt(B,x,u);
    if ( v_norm2(v_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_vec(pi1,x,z);
    mv_mlt(B,z,u);
    if ( v_norm2(v_sub(y,u,u)) >= MACHEPS*A->n )
	errmesg("px_cols()");

    MEMCHK();

    /* MATLAB save/load */
    notice("MATLAB save/load");
    A = m_resize(A,12,11);
    if ( (fp=fopen(SAVE_FILE,"w")) == (FILE *)NULL )
	printf("Cannot perform MATLAB save/load test\n");
    else
    {
	m_rand(A);
	m_save(fp, A, name);
	fclose(fp);
	if ( (fp=fopen(SAVE_FILE,"r")) == (FILE *)NULL )
	    printf("Cannot open save file \"%s\"\n",SAVE_FILE);
	else
	{
	    M_FREE(B);
	    B = m_load(fp,&cp);
	    if ( strcmp(name,cp) || m_norm1(m_sub(A,B,B)) >= MACHEPS*A->m )
		errmesg("mload()/m_save()");
	}
    }

    MEMCHK();

    /* Now, onto matrix factorisations */
    A = m_resize(A,10,10);
    B = m_resize(B,A->m,A->n);
    m_copy(A,B);
    x = v_resize(x,A->n);
    y = v_resize(y,A->m);
    z = v_resize(z,A->n);
    u = v_resize(u,A->m);
    v_rand(x);
    mv_mlt(B,x,y);
    z = v_copy(x,z);

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

    v_copy(y,x);
    tracecatch(LUsolve(A,pivot,x,x),"main");
    tracecatch(cond_est = LUcondest(A,pivot),"main");
    if ( v_norm2(v_sub(x,z,u)) >= MACHEPS*v_norm2(x)*cond_est)
    {
	errmesg("LUfactor()/LUsolve()");
	printf("# LU solution error = %g [cf MACHEPS = %g]\n",
	       v_norm2(v_sub(x,z,u)), MACHEPS);
    }

    vm_mlt(B,z,y);
    v_copy(y,x);
    tracecatch(LUTsolve(A,pivot,x,x),"main");
    if ( v_norm2(v_sub(x,z,u)) >= MACHEPS*v_norm2(x)*cond_est)
    {
	errmesg("LUfactor()/LUTsolve()");
	printf("# LU solution error = %g [cf MACHEPS = %g]\n",
	       v_norm2(v_sub(x,z,u)), MACHEPS);
    }
 
    MEMCHK();

    /* QR factorisation */
    m_copy(B,A);
    mv_mlt(B,z,y);
    notice("QR factor/solve:");
    diag = v_get(A->m);
    beta = v_get(A->m);
    QRfactor(A,diag);
    QRsolve(A,diag,y,x);
    if ( v_norm2(v_sub(x,z,u)) >= MACHEPS*v_norm2(x)*cond_est )
    {
	errmesg("QRfactor()/QRsolve()");
	printf("# QR solution error = %g [cf MACHEPS = %g]\n",
	       v_norm2(v_sub(x,z,u)), MACHEPS);
    }
    Q = m_get(A->m,A->m);
    makeQ(A,diag,Q);
    makeR(A,A);
    m_mlt(Q,A,C);
    m_sub(B,C,C);
    if ( m_norm1(C) >= MACHEPS*m_norm1(Q)*m_norm1(B) )
    {
	errmesg("QRfactor()/makeQ()/makeR()");
	printf("# QR reconstruction error = %g [cf MACHEPS = %g]\n",
	       m_norm1(C), MACHEPS);
    }

    MEMCHK();

    /* now try with a non-square matrix */
    A = m_resize(A,15,7);
    m_rand(A);
    B = m_copy(A,B);
    diag = v_resize(diag,A->n);
    beta = v_resize(beta,A->n);
    x = v_resize(x,A->n);
    y = v_resize(y,A->m);
    v_rand(y);
    QRfactor(A,diag);
    x = QRsolve(A,diag,y,x);
    /* z is the residual 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("QRfactor()/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_sub(B,C,C);
    if ( m_norm1(C) >= MACHEPS*m_norm1(Q)*m_norm1(B) )
    {
	errmesg("QRfactor()/makeQ()/makeR()");
	printf("# QR reconstruction error = %g [cf MACHEPS = %g]\n",
	       m_norm1(C), MACHEPS);
    }
    D = m_get(A->m,Q->m);
    mtrm_mlt(Q,Q,D);