ztorture.c

C语言版本的矩阵库

/**************************************************************************
**
** Copyright (C) 1993 David E. Steward & 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, complex routines
*/

static char rcsid[] = "$Id: $";

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


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

/* 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()

#define	checkpt()	printf("At line %d in file \"%s\"\n",__LINE__,__FILE__)

/* 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;

void	main(argc, argv)
int	argc;
char	*argv[];
{
    ZVEC 	*x = ZVNULL, *y = ZVNULL, *z = ZVNULL, *u = ZVNULL;
    ZVEC	*diag = ZVNULL;
    PERM	*pi1 = PNULL, *pi2 = PNULL, *pivot = PNULL;
    ZMAT	*A = ZMNULL, *B = ZMNULL, *C = ZMNULL, *D = ZMNULL,
	*Q = ZMNULL;
    complex	ONE;
    complex	z1, z2, z3;
    Real	cond_est, s1, s2, s3;
    int		i, 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);

    /* print out version information */
    m_version();

    printf("# Meschach Complex numbers & vectors torture test\n\n");
    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("\n");

    mem_stat_mark(1);

    notice("complex arithmetic & special functions");

    ONE = zmake(1.0,0.0);
    printf("# ONE = ");	z_output(ONE);
    z1.re = mrand();	z1.im = mrand();
    z2.re = mrand();	z2.im = mrand();
    z3 = zadd(z1,z2);
    if ( fabs(z1.re+z2.re-z3.re) + fabs(z1.im+z2.im-z3.im) > 10*MACHEPS )
	errmesg("zadd");
    z3 = zsub(z1,z2);
    if ( fabs(z1.re-z2.re-z3.re) + fabs(z1.im-z2.im-z3.im) > 10*MACHEPS )
	errmesg("zadd");
    z3 = zmlt(z1,z2);
    if ( fabs(z1.re*z2.re - z1.im*z2.im - z3.re) +
	 fabs(z1.im*z2.re + z1.re*z2.im - z3.im) > 10*MACHEPS )
	errmesg("zmlt");
    s1 = zabs(z1);
    if ( fabs(s1*s1 - (z1.re*z1.re+z1.im*z1.im)) > 10*MACHEPS )
	errmesg("zabs");
    if ( zabs(zsub(z1,zmlt(z2,zdiv(z1,z2)))) > 10*MACHEPS ||
	 zabs(zsub(ONE,zdiv(z1,zmlt(z2,zdiv(z1,z2))))) > 10*MACHEPS )
	errmesg("zdiv");

    z3 = zsqrt(z1);
    if ( zabs(zsub(z1,zmlt(z3,z3))) > 10*MACHEPS )
	errmesg("zsqrt");
    if ( zabs(zsub(z1,zlog(zexp(z1)))) > 10*MACHEPS )
	errmesg("zexp/zlog");
    

    printf("# Check: MACHEPS = %g\n",MACHEPS);
    /* allocate, initialise, copy and resize operations */
    /* ZVEC */
    notice("vector initialise, copy & resize");
    x = zv_get(12);
    y = zv_get(15);
    z = zv_get(12);
    zv_rand(x);
    zv_rand(y);
    z = zv_copy(x,z);
    if ( zv_norm2(zv_sub(x,z,z)) >= MACHEPS )
	errmesg("ZVEC copy");
    zv_copy(x,y);
    x = zv_resize(x,10);
    y = zv_resize(y,10);
    if ( zv_norm2(zv_sub(x,y,z)) >= MACHEPS )
	errmesg("ZVEC copy/resize");
    x = zv_resize(x,15);
    y = zv_resize(y,15);
    if ( zv_norm2(zv_sub(x,y,z)) >= MACHEPS )
	errmesg("VZEC resize");

    /* ZMAT */
    notice("matrix initialise, copy & resize");
    A = zm_get(8,5);
    B = zm_get(3,9);
    C = zm_get(8,5);
    zm_rand(A);
    zm_rand(B);
    C = zm_copy(A,C);
    if ( zm_norm_inf(zm_sub(A,C,C)) >= MACHEPS )
	errmesg("ZMAT copy");
    zm_copy(A,B);
    A = zm_resize(A,3,5);
    B = zm_resize(B,3,5);
    if ( zm_norm_inf(zm_sub(A,B,C)) >= MACHEPS )
	errmesg("ZMAT copy/resize");
    A = zm_resize(A,10,10);
    B = zm_resize(B,10,10);
    if ( zm_norm_inf(zm_sub(A,B,C)) >= MACHEPS )
	errmesg("ZMAT resize");

    MEMCHK();

    /* PERM */
    notice("permutation initialise, inverting & permuting vectors");
    pi1 = px_get(15);
    pi2 = px_get(12);
    px_rand(pi1);
    zv_rand(x);
    px_zvec(pi1,x,z);
    y = zv_resize(y,x->dim);
    pxinv_zvec(pi1,z,y);
    if ( zv_norm2(zv_sub(x,y,z)) >= MACHEPS )
	errmesg("PERMute vector");

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

    /* testing inner products and v_mltadd() etc */
    notice("inner products and linear combinations");
    u = zv_get(x->dim);
    zv_rand(u);
    zv_rand(x);
    zv_resize(y,x->dim);
    zv_rand(y);
    zv_mltadd(y,x,zneg(zdiv(zin_prod(x,y),zin_prod(x,x))),z);
    if ( zabs(zin_prod(x,z)) >= 5*MACHEPS*x->dim )
    {
	errmesg("zv_mltadd()/zin_prod()");
	printf("# error norm = %g\n", zabs(zin_prod(x,z)));
    }

    z1 = zneg(zdiv(zin_prod(x,y),zmake(zv_norm2(x)*zv_norm2(x),0.0)));
    zv_mlt(z1,x,u);
    zv_add(y,u,u);
    if ( zv_norm2(zv_sub(u,z,u)) >= MACHEPS*x->dim )
    {
	errmesg("zv_mlt()/zv_norm2()");
	printf("# error norm = %g\n", zv_norm2(u));
    }

#ifdef ANSI_C
    zv_linlist(u,x,z1,y,ONE,VNULL);
    if ( zv_norm2(zv_sub(u,z,u)) >= MACHEPS*x->dim )
	errmesg("zv_linlist()");
#endif
#ifdef VARARGS
    zv_linlist(u,x,z1,y,ONE,VNULL);
    if ( zv_norm2(zv_sub(u,z,u)) >= MACHEPS*x->dim )
	errmesg("zv_linlist()");
#endif

    MEMCHK();

    /* vector norms */
    notice("vector norms");
    x = zv_resize(x,12);
    zv_rand(x);
    for ( i = 0; i < x->dim; i++ )
	if ( zabs(zv_entry(x,i)) >= 0.7 )
	    zv_set_val(x,i,ONE);
        else
	    zv_set_val(x,i,zneg(ONE));
    s1 = zv_norm1(x);
    s2 = zv_norm2(x);	
    s3 = zv_norm_inf(x);
    if ( fabs(s1 - x->dim) >= MACHEPS*x->dim ||
	 fabs(s2 - sqrt((double)(x->dim))) >= MACHEPS*x->dim ||
	 fabs(s3 - 1.0) >= MACHEPS )
	errmesg("zv_norm1/2/_inf()");

    /* test matrix multiply etc */
    notice("matrix multiply and invert");
    A = zm_resize(A,10,10);
    B = zm_resize(B,10,10);
    zm_rand(A);
    zm_inverse(A,B);
    zm_mlt(A,B,C);
    for ( i = 0; i < C->m; i++ )
	zm_sub_val(C,i,i,ONE);
    if ( zm_norm_inf(C) >= MACHEPS*zm_norm_inf(A)*zm_norm_inf(B)*5 )
	errmesg("zm_inverse()/zm_mlt()");

    MEMCHK();

    /* ... and adjoints */
    notice("adjoints and adjoint-multiplies");
    zm_adjoint(A,A);	/* can do square matrices in situ */
    zmam_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("zm_adjoint()/zmam_mlt()");
    zm_adjoint(A,A);
    zm_adjoint(B,B);
    zmma_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("zm_adjoint()/zmma_mlt()");
    zsm_mlt(zmake(3.71,2.753),B,B);
    zmma_mlt(A,B,C);
    for ( i = 0; i < C->m; i++ )
	zm_set_val(C,i,i,zsub(zm_entry(C,i,i),zmake(3.71,-2.753)));
    if ( zm_norm_inf(C) >= MACHEPS*zm_norm_inf(A)*zm_norm_inf(B)*5 )
	errmesg("szm_mlt()/zmma_mlt()");
    zm_adjoint(B,B);
    zsm_mlt(zdiv(ONE,zmake(3.71,-2.753)),B,B);

    MEMCHK();

    /* ... and matrix-vector multiplies */
    notice("matrix-vector multiplies");
    x = zv_resize(x,A->n);
    y = zv_resize(y,A->m);
    z = zv_resize(z,A->m);
    u = zv_resize(u,A->n);
    zv_rand(x);
    zv_rand(y);
    zmv_mlt(A,x,z);
    z1 = zin_prod(y,z);
    zvm_mlt(A,y,u);
    z2 = zin_prod(u,x);
    if ( zabs(zsub(z1,z2)) >= (MACHEPS*x->dim)*x->dim )
    {
	errmesg("zmv_mlt()/zvm_mlt()");
	printf("# difference between inner products is %g\n",
	       zabs(zsub(z1,z2)));
    }
    zmv_mlt(B,z,u);
    if ( zv_norm2(zv_sub(u,x,u)) >= MACHEPS*zm_norm_inf(A)*zm_norm_inf(B)*5 )
	errmesg("zmv_mlt()/zvm_mlt()");

    MEMCHK();

    /* get/set row/col */
    notice("getting and setting rows and cols");
    x = zv_resize(x,A->n);
    y = zv_resize(y,B->m);
    x = zget_row(A,3,x);
    y = zget_col(B,3,y);
    if ( zabs(zsub(_zin_prod(x,y,0,Z_NOCONJ),ONE)) >=