#!/bin/sh # to unpack, sh this message in an empty directory PATH=/bin:/usr/bin cat > 24048P00 <<'bigmail CUT HERE............' - $(CC) $(CFLAGS) $(DEFS) -o tstmove tstmove.o \ - meschach.a $(LIBS) -tstpxvec: tstpxvec.o meschach.a - $(CC) $(CFLAGS) $(DEFS) -o tstpxvec tstpxvec.o \ - meschach.a $(LIBS) - //GO.SYSIN DD makefile echo FILELIST 1>&2 sed >FILELIST <<'//GO.SYSIN DD FILELIST' 's/^-//' --rw-r--r-- 1 0 30 09:47 FILELIST --rw-r--r-- 1 0 5 1994 README --rw-r--r-- 1 0 12 1994 arnoldi.c --rw-r--r-- 1 0 12 13:50 bdfactor.c --rw-r--r-- 1 0 12 13:44 bkpfacto.c --rw-r--r-- 1 0 12 13:45 chfactor.c --rwxr-xr-x 1 0 7 1994 configure --rw-r--r-- 1 0 7 1994 configure.in --rw-r--r-- 1 0 12 1994 conjgrad.c --rw-r--r-- 1 0 12 1994 copy.c --rw-r--r-- 1 0 12 1994 copyright --rw-r--r-- 1 0 12 1994 dmacheps.c --rw-r--r-- 1 0 30 08:49 err.c --rw-r--r-- 1 0 30 08:49 err.h --rw-r--r-- 1 0 18 1994 extras.c --rw-r--r-- 1 0 12 13:49 fft.c --rw-r--r-- 1 0 12 1994 fmacheps.c --rw-r--r-- 1 0 12 13:46 givens.c --rw-r--r-- 1 0 12 1994 hessen.c --rw-r--r-- 1 0 12 13:47 hsehldr.c --rw-r--r-- 1 0 12 1994 init.c --rw-r--r-- 1 0 13 1994 iotort.c --rw-r--r-- 1 0 7 1994 iter.h --rw-r--r-- 1 0 30 08:51 iter0.c --rw-r--r-- 1 0 30 08:55 iternsym.c --rw-r--r-- 1 0 30 08:57 itersym.c --rw-r--r-- 1 0 12 14:02 itertort.c --rw-r--r-- 1 0 12 1994 ivecop.c --rw-r--r-- 1 0 12 1994 lanczos.c --rw-r--r-- 1 0 12 1994 ls.dat --rw-r--r-- 1 0 12 13:42 lufactor.c --rw-r--r-- 1 0 24 1994 machine.c --rw-r--r-- 1 0 30 09:03 machine.h --rw-r--r-- 1 0 12 13:39 machine.h.in --rw-r--r-- 1 0 30 09:03 makefile --rw-r--r-- 1 0 22 1994 makefile.in --rw-r--r-- 1 0 12 1994 matlab.c --rw-r--r-- 1 0 20 09:39 matlab.h --rw-r--r-- 1 0 12 1994 matop.c --rw-r--r-- 1 0 15 1994 matrix.h --rw-r--r-- 1 0 12 1994 matrix2.h --rw-r--r-- 1 0 12 1994 matrixio.c --rw-r--r-- 1 0 12 1994 maxint.c --rw-r--r-- 1 0 12 1994 meminfo.c --rw-r--r-- 1 0 12 1994 meminfo.h --rw-r--r-- 1 0 4 1994 memory.c --rw-r--r-- 1 0 12 1994 memstat.c --rw-r--r-- 1 0 13 1994 memtort.c --rw-r--r-- 1 0 12 13:50 mfunc.c --rw-r--r-- 1 0 13 1994 mfuntort.c --rw-r--r-- 1 0 12 13:49 norm.c --rw-r--r-- 1 0 12 1994 oldnames.h --rw-r--r-- 1 0 12 1994 otherio.c --rw-r--r-- 1 0 23 1994 pxop.c --rw-r--r-- 1 0 12 13:47 qrfactor.c --rw-r--r-- 1 0 12 1994 rk4.dat --rw-r--r-- 1 0 12 13:45 schur.c --rw-r--r-- 1 0 12 13:48 solve.c --rw-r--r-- 1 0 7 1994 sparse.c --rw-r--r-- 1 0 12 1994 sparse.h --rw-r--r-- 1 0 12 1994 sparse2.h --rw-r--r-- 1 0 12 1994 sparseio.c --rw-r--r-- 1 0 12 13:52 spbkp.c --rw-r--r-- 1 0 12 13:52 spchfctr.c --rw-r--r-- 1 0 12 13:51 splufctr.c --rw-r--r-- 1 0 12 1994 sprow.c --rw-r--r-- 1 0 12 13:53 spswap.c --rw-r--r-- 1 0 28 1994 sptort.c --rw-r--r-- 1 0 12 1994 submat.c --rw-r--r-- 1 0 12 13:46 svd.c --rw-r--r-- 1 0 12 13:49 symmeig.c --rw-r--r-- 1 0 12 14:01 torture.c --rw-r--r-- 1 0 19 1994 tutadv.c --rw-r--r-- 1 0 19 1994 tutorial.c --rw-r--r-- 1 0 12 13:48 update.c --rw-r--r-- 1 0 7 1994 vecop.c --rw-r--r-- 1 0 23 1994 version.c --rw-r--r-- 1 0 12 1994 zcopy.c --rw-r--r-- 1 0 12 13:57 zfunc.c --rw-r--r-- 1 0 12 14:00 zgivens.c --rw-r--r-- 1 0 12 1994 zhessen.c --rw-r--r-- 1 0 12 13:59 zhsehldr.c --rw-r--r-- 1 0 12 13:57 zlufctr.c --rw-r--r-- 1 0 12 13:56 zmachine.c --rw-r--r-- 1 0 12 1994 zmatio.c --rw-r--r-- 1 0 12 1994 zmatlab.c --rw-r--r-- 1 0 12 1994 zmatop.c --rw-r--r-- 1 0 7 1994 zmatrix.h --rw-r--r-- 1 0 12 1994 zmatrix2.h --rw-r--r-- 1 0 22 1994 zmemory.c --rw-r--r-- 1 0 12 13:57 znorm.c --rw-r--r-- 1 0 12 13:57 zqrfctr.c --rw-r--r-- 1 0 12 13:57 zschur.c --rw-r--r-- 1 0 12 13:58 zsolve.c --rw-r--r-- 1 0 12 14:01 ztorture.c --rw-r--r-- 1 0 7 1994 zvecop.c - -DOC: -total 62 --rw------- 1 0 13 1994 fnindex.txt --rw------- 1 0 13 1994 tutorial.txt - -MACHINES: -total 6 -drwx------ 2 0 27 22:19 Cray -drwx------ 2 0 13 1994 GCC -drwx------ 2 0 2 1994 Linux -drwx------ 2 0 13 1994 RS6000 -drwx------ 2 0 27 22:15 SGI -drwx------ 2 0 13 1994 SPARC - -MACHINES/Cray: -total 15 --rw------- 1 0 27 11:18 machine.h --rw------- 1 0 27 11:22 makefile --rw------- 1 0 27 11:18 patch.1 --rw------- 1 0 27 11:18 patch.2 --rw------- 1 0 27 11:18 patch.3 - -MACHINES/GCC: -total 10 --rw------- 1 0 13 1994 machine.h --rw------- 1 0 13 1994 makefile - -MACHINES/Linux: -total 10 --rw------- 1 0 2 1994 machine.h --rw------- 1 0 2 1994 makefile - -MACHINES/RS6000: -total 16 --rw------- 1 0 24 1994 machine.c --rw------- 1 0 13 1994 machine.h --rw------- 1 0 13 1994 makefile - -MACHINES/SGI: -total 11 --rw------- 1 0 27 08:31 machine.h --rw------- 1 0 27 08:55 makefile - -MACHINES/SPARC: -total 10 --rw------- 1 0 13 1994 machine.h --rw------- 1 0 13 1994 makefile //GO.SYSIN DD FILELIST echo torture.c 1>&2 sed >torture.c <<'//GO.SYSIN DD torture.c' 's/^-//' - -/************************************************************************** -** -** 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: m2,v 1.1.1.1 1999/04/14 14:16:22 borland Exp $"; - -#include -#include "matrix2.h" -#include -#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); - 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; -} - - //GO.SYSIN DD torture.c echo sptort.c 1>&2 sed >sptort.c <<'//GO.SYSIN DD sptort.c' 's/^-//' - -/************************************************************************** -** -** 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 tests for the sparse matrix part of Meschach -*/ - -#include -#include -#include "matrix2.h" -#include "sparse2.h" -#include "iter.h" - -#define errmesg(mesg) printf("Error: %s error: line %d\n",mesg,__LINE__) -#define notice(mesg) printf("# Testing %s...\n",mesg); - -/* for iterative methods */ - -#if REAL == DOUBLE -#define EPS 1e-7 -#elif REAL == FLOAT -#define EPS 1e-3 -#endif - -int chk_col_access(A) -SPMAT *A; -{ - int i, j, nxt_idx, nxt_row, scan_cnt, total_cnt; - SPROW *r; - row_elt *e; - - if ( ! A ) - error(E_NULL,"chk_col_access"); - if ( ! A->flag_col ) - return FALSE; - - /* scan down each column, counting the number of entries met */ - scan_cnt = 0; - for ( j = 0; j < A->n; j++ ) - { - i = -1; - nxt_idx = A->start_idx[j]; - nxt_row = A->start_row[j]; - while ( nxt_row >= 0 && nxt_idx >= 0 && nxt_row > i ) - { - i = nxt_row; - r = &(A->row[i]); - e = &(r->elt[nxt_idx]); - nxt_idx = e->nxt_idx; - nxt_row = e->nxt_row; - scan_cnt++; - } - } - - total_cnt = 0; - for ( i = 0; i < A->m; i++ ) - total_cnt += A->row[i].len; - if ( total_cnt != scan_cnt ) - return FALSE; - else - return TRUE; -} - - -void main(argc, argv) -int argc; -char *argv[]; -{ - VEC *x, *y, *z, *u, *v; - Real s1, s2; - PERM *pivot; - SPMAT *A, *B, *C; - SPMAT *B1, *C1; - SPROW *r; - int i, j, k, deg, seed, m, m_old, n, n_old; - - - mem_info_on(TRUE); - - setbuf(stdout, (char *)NULL); - /* get seed if in argument list */ - if ( argc == 1 ) - seed = 1111; - else if ( argc == 2 && sscanf(argv[1],"%d",&seed) == 1 ) - ; - else - { - printf("usage: %s [seed]\n", argv[0]); - exit(0); - } - srand(seed); - - /* set up two random sparse matrices */ - m = 120; - n = 100; - deg = 8; - notice("allocating sparse matrices"); - A = sp_get(m,n,deg); - B = sp_get(m,n,deg); - notice("setting and getting matrix entries"); - for ( k = 0; k < m*deg; k++ ) - { - i = (rand() >> 8) % m; - j = (rand() >> 8) % n; - sp_set_val(A,i,j,rand()/((Real)MAX_RAND)); - i = (rand() >> 8) % m; - j = (rand() >> 8) % n; - sp_set_val(B,i,j,rand()/((Real)MAX_RAND)); - } - for ( k = 0; k < 10; k++ ) - { - s1 = rand()/((Real)MAX_RAND); - i = (rand() >> 8) % m; - j = (rand() >> 8) % n; - sp_set_val(A,i,j,s1); - s2 = sp_get_val(A,i,j); - if ( fabs(s1 - s2) >= MACHEPS ) - break; - } - if ( k < 10 ) - errmesg("sp_set_val()/sp_get_val()"); - - /* test copy routines */ - notice("copy routines"); - x = v_get(n); - y = v_get(m); - z = v_get(m); - /* first copy routine */ - C = sp_copy(A); - for ( k = 0; k < 100; k++ ) - { - v_rand(x); - sp_mv_mlt(A,x,y); - sp_mv_mlt(C,x,z); - if ( v_norm_inf(v_sub(y,z,z)) >= MACHEPS*deg*m ) - break; - } - if ( k < 100 ) - { - errmesg("sp_copy()/sp_mv_mlt()"); - printf("# Error in A.x (inf norm) = %g [cf MACHEPS = %g]\n", - v_norm_inf(z), MACHEPS); - } - /* second copy routine - -- note that A & B have different sparsity patterns */ - - mem_stat_mark(1); - sp_copy2(A,B); - mem_stat_free(1); - for ( k = 0; k < 10; k++ ) - { - v_rand(x); - sp_mv_mlt(A,x,y); - sp_mv_mlt(B,x,z); - if ( v_norm_inf(v_sub(y,z,z)) >= MACHEPS*deg*m ) - break; - } - if ( k < 10 ) - { - errmesg("sp_copy2()/sp_mv_mlt()"); - printf("# Error in A.x (inf norm) = %g [cf MACHEPS = %g]\n", - v_norm_inf(z), MACHEPS); - } - - /* now check compacting routine */ - notice("compacting routine"); - sp_compact(B,0.0); - for ( k = 0; k < 10; k++ ) - { - v_rand(x); - sp_mv_mlt(A,x,y); - sp_mv_mlt(B,x,z); - if ( v_norm_inf(v_sub(y,z,z)) >= MACHEPS*deg*m ) - break; - } - if ( k < 10 ) - { - errmesg("sp_compact()"); - printf("# Error in A.x (inf norm) = %g [cf MACHEPS = %g]\n", - v_norm_inf(z), MACHEPS); - } - for ( i = 0; i < B->m; i++ ) - { - r = &(B->row[i]); - for ( j = 0; j < r->len; j++ ) - if ( r->elt[j].val == 0.0 ) - break; - } - if ( i < B->m ) - { - errmesg("sp_compact()"); - printf("# Zero entry in compacted matrix\n"); - } - - /* check column access paths */ - notice("resizing and access paths"); - m_old = A->m-1; - n_old = A->n-1; - A = sp_resize(A,A->m+10,A->n+10); - for ( k = 0 ; k < 20; k++ ) - { - i = m_old + ((rand() >> 8) % 10); - j = n_old + ((rand() >> 8) % 10); - s1 = rand()/((Real)MAX_RAND); - sp_set_val(A,i,j,s1); - if ( fabs(s1 - sp_get_val(A,i,j)) >= MACHEPS ) - break; - } - if ( k < 20 ) - errmesg("sp_resize()"); - sp_col_access(A); - if ( ! chk_col_access(A) ) - { - errmesg("sp_col_access()"); - } - sp_diag_access(A); - for ( i = 0; i < A->m; i++ ) - { - r = &(A->row[i]); - if ( r->diag != sprow_idx(r,i) ) - break; - } - if ( i < A->m ) - { - errmesg("sp_diag_access()"); - } - - /* test both sp_mv_mlt() and sp_vm_mlt() */ - x = v_resize(x,B->n); - y = v_resize(y,B->m); - u = v_get(B->m); - v = v_get(B->n); - for ( k = 0; k < 10; k++ ) - { - v_rand(x); - v_rand(y); - sp_mv_mlt(B,x,u); - sp_vm_mlt(B,y,v); - if ( fabs(in_prod(x,v) - in_prod(y,u)) >= - MACHEPS*v_norm2(x)*v_norm2(u)*5 ) - break; - } - if ( k < 10 ) - { - errmesg("sp_mv_mlt()/sp_vm_mlt()"); - printf("# Error in inner products = %g [cf MACHEPS = %g]\n", - fabs(in_prod(x,v) - in_prod(y,u)), MACHEPS); - } - - SP_FREE(A); - SP_FREE(B); - SP_FREE(C); - - /* now test Cholesky and LU factorise and solve */ - notice("sparse Cholesky factorise/solve"); - A = iter_gen_sym(120,8); - B = sp_copy(A); - spCHfactor(A); - x = v_resize(x,A->m); - y = v_resize(y,A->m); - v_rand(x); - sp_mv_mlt(B,x,y); - z = v_resize(z,A->m); - spCHsolve(A,y,z); - v = v_resize(v,A->m); - sp_mv_mlt(B,z,v); - /* compute residual */ - v_sub(y,v,v); - if ( v_norm2(v) >= MACHEPS*v_norm2(y)*10 ) - { - errmesg("spCHfactor()/spCHsolve()"); - printf("# Sparse Cholesky residual = %g [cf MACHEPS = %g]\n", - v_norm2(v), MACHEPS); - } - /* compute error in solution */ - v_sub(x,z,z); - if ( v_norm2(z) > MACHEPS*v_norm2(x)*10 ) - { - errmesg("spCHfactor()/spCHsolve()"); - printf("# Solution error = %g [cf MACHEPS = %g]\n", - v_norm2(z), MACHEPS); - } - - /* now test symbolic and incomplete factorisation */ - SP_FREE(A); - A = sp_copy(B); - - mem_stat_mark(2); - spCHsymb(A); - mem_stat_mark(2); - - spICHfactor(A); - spCHsolve(A,y,z); - v = v_resize(v,A->m); - sp_mv_mlt(B,z,v); - /* compute residual */ - v_sub(y,v,v); - if ( v_norm2(v) >= MACHEPS*v_norm2(y)*5 ) - { - errmesg("spCHsymb()/spICHfactor()"); - printf("# Sparse Cholesky residual = %g [cf MACHEPS = %g]\n", - v_norm2(v), MACHEPS); - } - /* compute error in solution */ - v_sub(x,z,z); - if ( v_norm2(z) > MACHEPS*v_norm2(x)*10 ) - { - errmesg("spCHsymb()/spICHfactor()"); - printf("# Solution error = %g [cf MACHEPS = %g]\n", - v_norm2(z), MACHEPS); - } - - /* now test sparse LU factorisation */ - notice("sparse LU factorise/solve"); - SP_FREE(A); - SP_FREE(B); - A = iter_gen_nonsym(100,100,8,1.0); - - B = sp_copy(A); - x = v_resize(x,A->n); - z = v_resize(z,A->n); - y = v_resize(y,A->m); - v = v_resize(v,A->m); - - v_rand(x); - sp_mv_mlt(B,x,y); - pivot = px_get(A->m); - - mem_stat_mark(3); - spLUfactor(A,pivot,0.1); - spLUsolve(A,pivot,y,z); - mem_stat_free(3); - sp_mv_mlt(B,z,v); - - /* compute residual */ - v_sub(y,v,v); - if ( v_norm2(v) >= MACHEPS*v_norm2(y)*A->m ) - { - errmesg("spLUfactor()/spLUsolve()"); - printf("# Sparse LU residual = %g [cf MACHEPS = %g]\n", - v_norm2(v), MACHEPS); - } - /* compute error in solution */ - v_sub(x,z,z); - if ( v_norm2(z) > MACHEPS*v_norm2(x)*100*A->m ) - { - errmesg("spLUfactor()/spLUsolve()"); - printf("# Sparse LU solution error = %g [cf MACHEPS = %g]\n", - v_norm2(z), MACHEPS); - } - - /* now check spLUTsolve */ - mem_stat_mark(4); - sp_vm_mlt(B,x,y); - spLUTsolve(A,pivot,y,z); - sp_vm_mlt(B,z,v); - mem_stat_free(4); - - /* compute residual */ - v_sub(y,v,v); - if ( v_norm2(v) >= MACHEPS*v_norm2(y)*A->m ) - { - errmesg("spLUTsolve()"); - printf("# Sparse LU residual = %g [cf MACHEPS = %g]\n", - v_norm2(v), MACHEPS); - } - /* compute error in solution */ - v_sub(x,z,z); - if ( v_norm2(z) > MACHEPS*v_norm2(x)*100*A->m ) - { - errmesg("spLUTsolve()"); - printf("# Sparse LU solution error = %g [cf MACHEPS = %g]\n", - v_norm2(z), MACHEPS); - } - - /* algebraic operations */ - notice("addition,subtraction and multiplying by a number"); - SP_FREE(A); - SP_FREE(B); - - m = 120; - n = 120; - deg = 5; - A = sp_get(m,n,deg); - B = sp_get(m,n,deg); - C = sp_get(m,n,deg); - C1 = sp_get(m,n,deg); - - for ( k = 0; k < m*deg; k++ ) - { - i = (rand() >> 8) % m; - j = (rand() >> 8) % n; - sp_set_val(A,i,j,rand()/((Real)MAX_RAND)); - i = (rand() >> 8) % m; - j = (rand() >> 8) % n; - sp_set_val(B,i,j,rand()/((Real)MAX_RAND)); - } - - s1 = mrand(); - B1 = sp_copy(B); - - mem_stat_mark(1); - sp_smlt(B,s1,C); - sp_add(A,C,C1); - sp_sub(C1,A,C); - sp_smlt(C,-1.0/s1,C); - sp_add(C,B1,C); - - s2 = 0.0; - for (k=0; k < C->m; k++) { - r = &(C->row[k]); - for (j=0; j < r->len; j++) { - if (s2 < fabs(r->elt[j].val)) - s2 = fabs(r->elt[j].val); - } - } - - if (s2 > MACHEPS*A->m) { - errmesg("add, sub, mlt sparse matrices (args not in situ)\n"); - printf(" difference = %g [MACEPS = %g]\n",s2,MACHEPS); - } - - sp_mltadd(A,B1,s1,C1); - sp_sub(C1,A,A); - sp_smlt(A,1.0/s1,C1); - sp_sub(C1,B1,C1); - mem_stat_free(1); - - s2 = 0.0; - for (k=0; k < C1->m; k++) { - r = &(C1->row[k]); - for (j=0; j < r->len; j++) { - if (s2 < fabs(r->elt[j].val)) - s2 = fabs(r->elt[j].val); - } - } - - if (s2 > MACHEPS*A->m) { - errmesg("add, sub, mlt sparse matrices (args not in situ)\n"); - printf(" difference = %g [MACEPS = %g]\n",s2,MACHEPS); - } - - V_FREE(x); - V_FREE(y); - V_FREE(z); - V_FREE(u); - V_FREE(v); - PX_FREE(pivot); - SP_FREE(A); - SP_FREE(B); - SP_FREE(C); - SP_FREE(B1); - SP_FREE(C1); - - printf("# Done testing (%s)\n",argv[0]); - mem_info(); -} - - - - - //GO.SYSIN DD sptort.c echo ztorture.c 1>&2 sed >ztorture.c <<'//GO.SYSIN DD ztorture.c' 's/^-//' - -/************************************************************************** -** -** 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: m2,v 1.1.1.1 1999/04/14 14:16:22 borland Exp $"; - -#include -#include "zmatrix2.h" -#include -#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() - -/* 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(v_entry(x,i)) >= 0.7 ) - v_set_val(x,i,ONE); - else - v_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++ ) - m_set_val(C,i,i,zsub(m_entry(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++ ) - m_set_val(C,i,i,zsub(m_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++ ) - m_set_val(C,i,i,zsub(m_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++ ) - m_set_val(C,i,i,zsub(m_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)) >= - 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++ ) - m_set_val(C,i,i,zsub(m_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,m_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++ ) - m_set_val(D,i,i,zsub(m_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); - MEMCHK(); - - 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++ ) - m_set_val(D,i,i,zsub(m_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(); -} - //GO.SYSIN DD ztorture.c echo memtort.c 1>&2 sed >memtort.c <<'//GO.SYSIN DD memtort.c' 's/^-//' - -/************************************************************************** -** -** 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. -** -***************************************************************************/ - - -/* - Tests for mem_info.c functions - */ - -static char rcsid[] = "$Id: m2,v 1.1.1.1 1999/04/14 14:16:22 borland Exp $"; - -#include -#include -#include "matrix2.h" -#include "sparse2.h" -#include "zmatrix2.h" - - -#define errmesg(mesg) printf("Error: %s error: line %d\n",mesg,__LINE__) -#define notice(mesg) printf("# Testing %s...\n",mesg) - - -/* new types list */ - -extern MEM_CONNECT mem_connect[MEM_CONNECT_MAX_LISTS]; - -/* the number of a new list */ -#define FOO_LIST 1 - -/* numbers of types */ -#define TYPE_FOO_1 1 -#define TYPE_FOO_2 2 - -typedef struct { - int dim; - int fix_dim; - Real (*a)[10]; -} FOO_1; - -typedef struct { - int dim; - int fix_dim; - Real (*a)[2]; -} FOO_2; - - - -FOO_1 *foo_1_get(dim) -int dim; -{ - FOO_1 *f; - - if ((f = (FOO_1 *)malloc(sizeof(FOO_1))) == NULL) - error(E_MEM,"foo_1_get"); - else if (mem_info_is_on()) { - mem_bytes_list(TYPE_FOO_1,0,sizeof(FOO_1),FOO_LIST); - mem_numvar_list(TYPE_FOO_1,1,FOO_LIST); - } - - f->dim = dim; - f->fix_dim = 10; - if ((f->a = (Real (*)[10])malloc(dim*sizeof(Real [10]))) == NULL) - error(E_MEM,"foo_1_get"); - else if (mem_info_is_on()) - mem_bytes_list(TYPE_FOO_1,0,dim*sizeof(Real [10]),FOO_LIST); - - return f; -} - - -FOO_2 *foo_2_get(dim) -int dim; -{ - FOO_2 *f; - - if ((f = (FOO_2 *)malloc(sizeof(FOO_2))) == NULL) - error(E_MEM,"foo_2_get"); - else if (mem_info_is_on()) { - mem_bytes_list(TYPE_FOO_2,0,sizeof(FOO_2),FOO_LIST); - mem_numvar_list(TYPE_FOO_2,1,FOO_LIST); - } - - f->dim = dim; - f->fix_dim = 2; - if ((f->a = (Real (*)[2])malloc(dim*sizeof(Real [2]))) == NULL) - error(E_MEM,"foo_2_get"); - else if (mem_info_is_on()) - mem_bytes_list(TYPE_FOO_2,0,dim*sizeof(Real [2]),FOO_LIST); - - return f; -} - - - -int foo_1_free(f) -FOO_1 *f; -{ - if ( f != NULL) { - if (mem_info_is_on()) { - mem_bytes_list(TYPE_FOO_1,sizeof(FOO_1)+ - f->dim*sizeof(Real [10]),0,FOO_LIST); - mem_numvar_list(TYPE_FOO_1,-1,FOO_LIST); - } - - free(f->a); - free(f); - } - return 0; -} - -int foo_2_free(f) -FOO_2 *f; -{ - if ( f != NULL) { - if (mem_info_is_on()) { - mem_bytes_list(TYPE_FOO_2,sizeof(FOO_2)+ - f->dim*sizeof(Real [2]),0,FOO_LIST); - mem_numvar_list(TYPE_FOO_2,-1,FOO_LIST); - } - - free(f->a); - free(f); - } - return 0; -} - - - - -char *foo_type_name[] = { - "nothing", - "FOO_1", - "FOO_2" -}; - - -#define FOO_NUM_TYPES (sizeof(foo_type_name)/sizeof(*foo_type_name)) - - -int (*foo_free_func[FOO_NUM_TYPES])() = { - NULL, - foo_1_free, - foo_2_free - }; - - - -static MEM_ARRAY foo_info_sum[FOO_NUM_TYPES]; - - - - /* 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; -} - -#ifdef SPARSE -SPMAT *gen_non_symm(m,n) -int m, n; -{ - SPMAT *A; - static PERM *px = PNULL; - int i, j, k, k_max; - Real s1; - - A = sp_get(m,n,8); - px = px_resize(px,n); - MEM_STAT_REG(px,TYPE_PERM); - for ( i = 0; i < A->m; i++ ) - { - k_max = 1 + ((rand() >> 8) % 10); - for ( k = 0; k < k_max; k++ ) - { - j = (rand() >> 8) % A->n; - s1 = rand()/((double)MAX_RAND); - sp_set_val(A,i,j,s1); - } - } - /* to make it likely that A is nonsingular, use pivot... */ - for ( i = 0; i < 2*A->n; i++ ) - { - j = (rand() >> 8) % A->n; - k = (rand() >> 8) % A->n; - px_transp(px,j,k); - } - for ( i = 0; i < A->n; i++ ) - sp_set_val(A,i,px->pe[i],1.0); - - - return A; -} -#endif - -void stat_test1(par) -int par; -{ - static MAT *AT = MNULL; - static VEC *xt1 = VNULL, *yt1 = VNULL; - static VEC *xt2 = VNULL, *yt2 = VNULL; - static VEC *xt3 = VNULL, *yt3 = VNULL; - static VEC *xt4 = VNULL, *yt4 = VNULL; - - AT = m_resize(AT,10,10); - xt1 = v_resize(xt1,10); - yt1 = v_resize(yt1,10); - xt2 = v_resize(xt2,10); - yt2 = v_resize(yt2,10); - xt3 = v_resize(xt3,10); - yt3 = v_resize(yt3,10); - xt4 = v_resize(xt4,10); - yt4 = v_resize(yt4,10); - - MEM_STAT_REG(AT,TYPE_MAT); - -#ifdef ANSI_C - mem_stat_reg_vars(0,TYPE_VEC,&xt1,&xt2,&xt3,&xt4,&yt1, - &yt2,&yt3,&yt4,NULL); -#else -#ifdef VARARGS - mem_stat_reg_vars(0,TYPE_VEC,&xt1,&xt2,&xt3,&xt4,&yt1, - &yt2,&yt3,&yt4,NULL); -#else - MEM_STAT_REG(xt1,TYPE_VEC); - MEM_STAT_REG(yt1,TYPE_VEC); - MEM_STAT_REG(xt2,TYPE_VEC); - MEM_STAT_REG(yt2,TYPE_VEC); - MEM_STAT_REG(xt3,TYPE_VEC); - MEM_STAT_REG(yt3,TYPE_VEC); - MEM_STAT_REG(xt4,TYPE_VEC); - MEM_STAT_REG(yt4,TYPE_VEC); -#endif -#endif - - v_rand(xt1); - m_rand(AT); - mv_mlt(AT,xt1,yt1); - -} - - -void stat_test2(par) -int par; -{ - static PERM *px = PNULL; - static IVEC *ixt = IVNULL, *iyt = IVNULL; - - px = px_resize(px,10); - ixt = iv_resize(ixt,10); - iyt = iv_resize(iyt,10); - - MEM_STAT_REG(px,TYPE_PERM); - MEM_STAT_REG(ixt,TYPE_IVEC); - MEM_STAT_REG(iyt,TYPE_IVEC); - - px_rand(px); - px_inv(px,px); -} - -#ifdef SPARSE -void stat_test3(par) -int par; -{ - static SPMAT *AT = (SPMAT *)NULL; - static VEC *xt = VNULL, *yt = VNULL; - static SPROW *r = (SPROW *) NULL; - - if (AT == (SPMAT *)NULL) - AT = gen_non_symm(100,100); - else - AT = sp_resize(AT,100,100); - xt = v_resize(xt,100); - yt = v_resize(yt,100); - if (r == NULL) r = sprow_get(100); - - MEM_STAT_REG(AT,TYPE_SPMAT); - MEM_STAT_REG(xt,TYPE_VEC); - MEM_STAT_REG(yt,TYPE_VEC); - MEM_STAT_REG(r,TYPE_SPROW); - - v_rand(xt); - sp_mv_mlt(AT,xt,yt); - -} -#endif - -#ifdef COMPLEX -void stat_test4(par) -int par; -{ - static ZMAT *AT = ZMNULL; - static ZVEC *xt = ZVNULL, *yt = ZVNULL; - - AT = zm_resize(AT,10,10); - xt = zv_resize(xt,10); - yt = zv_resize(yt,10); - - MEM_STAT_REG(AT,TYPE_ZMAT); - MEM_STAT_REG(xt,TYPE_ZVEC); - MEM_STAT_REG(yt,TYPE_ZVEC); - - zv_rand(xt); - zm_rand(AT); - zmv_mlt(AT,xt,yt); - -} -#endif - - -void main(argc, argv) -int argc; -char *argv[]; -{ - VEC *x = VNULL, *y = VNULL, *z = VNULL; - PERM *pi1 = PNULL, *pi2 = PNULL, *pi3 = PNULL; - MAT *A = MNULL, *B = MNULL, *C = MNULL; -#ifdef SPARSE - SPMAT *sA, *sB; - SPROW *r; -#endif - IVEC *ix = IVNULL, *iy = IVNULL, *iz = IVNULL; - int m,n,i,j,deg,k; - Real s1,s2; -#ifdef COMPLEX - ZVEC *zx = ZVNULL, *zy = ZVNULL, *zz = ZVNULL; - ZMAT *zA = ZMNULL, *zB = ZMNULL, *zC = ZMNULL; - complex ONE; -#endif - /* variables for testing attaching new lists of types */ - FOO_1 *foo_1; - FOO_2 *foo_2; - - - mem_info_on(TRUE); - -#if defined(ANSI_C) || defined(VARARGS) - - notice("vector initialize, copy & resize"); - - n = v_get_vars(15,&x,&y,&z,(VEC **)NULL); - if (n != 3) { - errmesg("v_get_vars"); - printf(" n = %d (should be 3)\n",n); - } - - v_rand(x); - v_rand(y); - z = v_copy(x,z); - if ( v_norm2(v_sub(x,z,z)) >= MACHEPS ) - errmesg("v_get_vars"); - v_copy(x,y); - n = v_resize_vars(10,&x,&y,&z,NULL); - if ( n != 3 || v_norm2(v_sub(x,y,z)) >= MACHEPS ) - errmesg("VEC copy/resize"); - - n = v_resize_vars(20,&x,&y,&z,NULL); - if ( n != 3 || v_norm2(v_sub(x,y,z)) >= MACHEPS ) - errmesg("VEC resize"); - - n = v_free_vars(&x,&y,&z,NULL); - if (n != 3) - errmesg("v_free_vars"); - - /* IVEC */ - notice("int vector initialise, copy & resize"); - n = iv_get_vars(15,&ix,&iy,&iz,NULL); - - if (n != 3) { - errmesg("iv_get_vars"); - printf(" n = %d (should be 3)\n",n); - } - for (i=0; i < ix->dim; i++) { - ix->ive[i] = 2*i-1; - iy->ive[i] = 3*i+2; - } - iz = iv_add(ix,iy,iz); - for (i=0; i < ix->dim; i++) - if ( iz->ive[i] != 5*i+1) - errmesg("iv_get_vars"); - - n = iv_resize_vars(10,&ix,&iy,&iz,NULL); - if ( n != 3) errmesg("IVEC copy/resize"); - - iv_add(ix,iy,iz); - for (i=0; i < ix->dim; i++) - if (iz->ive[i] != 5*i+1) - errmesg("IVEC copy/resize"); - - n = iv_resize_vars(20,&ix,&iy,&iz,NULL); - if ( n != 3 ) errmesg("IVEC resize"); - - iv_add(ix,iy,iz); - for (i=0; i < 10; i++) - if (iz->ive[i] != 5*i+1) - errmesg("IVEC copy/resize"); - - n = iv_free_vars(&ix,&iy,&iz,NULL); - if (n != 3) - errmesg("iv_free_vars"); - - /* MAT */ - notice("matrix initialise, copy & resize"); - n = m_get_vars(10,10,&A,&B,&C,NULL); - if (n != 3) { - errmesg("m_get_vars"); - printf(" n = %d (should be 3)\n",n); - } - - 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); - n = m_resize_vars(5,5,&A,&B,&C,NULL); - if ( n != 3 || m_norm_inf(m_sub(A,B,C)) >= MACHEPS ) - errmesg("MAT copy/resize"); - - n = m_resize_vars(20,20,&A,&B,NULL); - if ( m_norm_inf(m_sub(A,B,C)) >= MACHEPS ) - errmesg("MAT resize"); - - k = m_free_vars(&A,&B,&C,NULL); - if ( k != 3 ) - errmesg("MAT free"); - - /* PERM */ - notice("permutation initialise, inverting & permuting vectors"); - n = px_get_vars(15,&pi1,&pi2,&pi3,NULL); - if (n != 3) { - errmesg("px_get_vars"); - printf(" n = %d (should be 3)\n",n); - } - - v_get_vars(15,&x,&y,&z,NULL); - - 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,pi3); - for ( i = 0; i < pi3->size; i++ ) - if ( pi3->pe[i] != i ) - errmesg("PERM inverse/multiply"); - - px_resize_vars(20,&pi1,&pi2,&pi3,NULL); - v_resize_vars(20,&x,&y,&z,NULL); - - px_rand(pi1); - v_rand(x); - px_vec(pi1,x,z); - 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,pi3); - for ( i = 0; i < pi3->size; i++ ) - if ( pi3->pe[i] != i ) - errmesg("PERM inverse/multiply"); - - n = px_free_vars(&pi1,&pi2,&pi3,NULL); - if ( n != 3 ) - errmesg("PERM px_free_vars"); - -#ifdef SPARSE - /* set up two random sparse matrices */ - m = 120; - n = 100; - deg = 5; - notice("allocating sparse matrices"); - k = sp_get_vars(m,n,deg,&sA,&sB,NULL); - if (k != 2) { - errmesg("sp_get_vars"); - printf(" n = %d (should be 2)\n",k); - } - - notice("setting and getting matrix entries"); - for ( k = 0; k < m*deg; k++ ) - { - i = (rand() >> 8) % m; - j = (rand() >> 8) % n; - sp_set_val(sA,i,j,rand()/((Real)MAX_RAND)); - i = (rand() >> 8) % m; - j = (rand() >> 8) % n; - sp_set_val(sB,i,j,rand()/((Real)MAX_RAND)); - } - for ( k = 0; k < 10; k++ ) - { - s1 = rand()/((Real)MAX_RAND); - i = (rand() >> 8) % m; - j = (rand() >> 8) % n; - sp_set_val(sA,i,j,s1); - s2 = sp_get_val(sA,i,j); - if ( fabs(s1 - s2) >= MACHEPS ) { - printf(" s1 = %g, s2 = %g, |s1 - s2| = %g\n", - s1,s2,fabs(s1-s2)); - break; - } - } - if ( k < 10 ) - errmesg("sp_set_val()/sp_get_val()"); - - /* check column access paths */ - notice("resizing and access paths"); - k = sp_resize_vars(sA->m+10,sA->n+10,&sA,&sB,NULL); - if (k != 2) { - errmesg("sp_get_vars"); - printf(" n = %d (should be 2)\n",k); - } - - for ( k = 0 ; k < 20; k++ ) - { - i = sA->m - 1 - ((rand() >> 8) % 10); - j = sA->n - 1 - ((rand() >> 8) % 10); - s1 = rand()/((Real)MAX_RAND); - sp_set_val(sA,i,j,s1); - if ( fabs(s1 - sp_get_val(sA,i,j)) >= MACHEPS ) - break; - } - if ( k < 20 ) - errmesg("sp_resize()"); - sp_col_access(sA); - if ( ! chk_col_access(sA) ) - { - errmesg("sp_col_access()"); - } - sp_diag_access(sA); - for ( i = 0; i < sA->m; i++ ) - { - r = &(sA->row[i]); - if ( r->diag != sprow_idx(r,i) ) - break; - } - if ( i < sA->m ) - { - errmesg("sp_diag_access()"); - } - - k = sp_free_vars(&sA,&sB,NULL); - if (k != 2) - errmesg("sp_free_vars"); -#endif /* SPARSE */ - - -#ifdef COMPLEX - /* complex stuff */ - - ONE = zmake(1.0,0.0); - printf("# ONE = "); z_output(ONE); - printf("# Check: MACHEPS = %g\n",MACHEPS); - /* allocate, initialise, copy and resize operations */ - /* ZVEC */ - notice("vector initialise, copy & resize"); - zv_get_vars(12,&zx,&zy,&zz,NULL); - - zv_rand(zx); - zv_rand(zy); - zz = zv_copy(zx,zz); - if ( zv_norm2(zv_sub(zx,zz,zz)) >= MACHEPS ) - errmesg("ZVEC copy"); - zv_copy(zx,zy); - - zv_resize_vars(10,&zx,&zy,NULL); - if ( zv_norm2(zv_sub(zx,zy,zz)) >= MACHEPS ) - errmesg("ZVEC copy/resize"); - - zv_resize_vars(20,&zx,&zy,NULL); - if ( zv_norm2(zv_sub(zx,zy,zz)) >= MACHEPS ) - errmesg("VZEC resize"); - zv_free_vars(&zx,&zy,&zz,NULL); - - - /* ZMAT */ - notice("matrix initialise, copy & resize"); - zm_get_vars(8,5,&zA,&zB,&zC,NULL); - - zm_rand(zA); - zm_rand(zB); - zC = zm_copy(zA,zC); - if ( zm_norm_inf(zm_sub(zA,zC,zC)) >= MACHEPS ) - errmesg("ZMAT copy"); - - zm_copy(zA,zB); - zm_resize_vars(3,5,&zA,&zB,&zC,NULL); - - if ( zm_norm_inf(zm_sub(zA,zB,zC)) >= MACHEPS ) - errmesg("ZMAT copy/resize"); - zm_resize_vars(20,20,&zA,&zB,&zC,NULL); - - if ( zm_norm_inf(zm_sub(zA,zB,zC)) >= MACHEPS ) - errmesg("ZMAT resize"); - - zm_free_vars(&zA,&zB,&zC,NULL); -#endif /* COMPLEX */ - -#endif /* if defined(ANSI_C) || defined(VARARGS) */ - - printf("# test of mem_info_bytes and mem_info_numvar\n"); - printf(" TYPE VEC: %ld bytes allocated, %d variables allocated\n", - mem_info_bytes(TYPE_VEC,0),mem_info_numvar(TYPE_VEC,0)); - - notice("static memory test"); - mem_info_on(TRUE); - mem_stat_mark(1); - for (i=0; i < 100; i++) - stat_test1(i); - mem_stat_free(1); - - mem_stat_mark(1); - for (i=0; i < 100; i++) { - stat_test1(i); -#ifdef COMPLEX - stat_test4(i); -#endif - } - - mem_stat_mark(2); - for (i=0; i < 100; i++) - stat_test2(i); - - mem_stat_mark(3); -#ifdef SPARSE - for (i=0; i < 100; i++) - stat_test3(i); -#endif - - mem_info(); - mem_dump_list(stdout,0); - - mem_stat_free(1); - mem_stat_free(3); - mem_stat_mark(4); - - for (i=0; i < 100; i++) { - stat_test1(i); -#ifdef COMPLEX - stat_test4(i); -#endif - } - - mem_stat_dump(stdout,0); - if (mem_stat_show_mark() != 4) { - errmesg("not 4 in mem_stat_show_mark()"); - } - - mem_stat_free(2); - mem_stat_free(4); - - if (mem_stat_show_mark() != 0) { - errmesg("not 0 in mem_stat_show_mark()"); - } - - /* add new list of types */ - - mem_attach_list(FOO_LIST,FOO_NUM_TYPES,foo_type_name, - foo_free_func,foo_info_sum); - if (!mem_is_list_attached(FOO_LIST)) - errmesg("list FOO_LIST is not attached"); - - mem_dump_list(stdout,FOO_LIST); - foo_1 = foo_1_get(6); - foo_2 = foo_2_get(3); - for (i=0; i < foo_1->dim; i++) - for (j=0; j < foo_1->fix_dim; j++) - foo_1->a[i][j] = i+j; - for (i=0; i < foo_2->dim; i++) - for (j=0; j < foo_2->fix_dim; j++) - foo_2->a[i][j] = i+j; - printf(" foo_1->a[%d][%d] = %g\n",5,9,foo_1->a[5][9]); - printf(" foo_2->a[%d][%d] = %g\n",2,1,foo_2->a[2][1]); - - mem_stat_mark(5); - mem_stat_reg_list((void **)&foo_1,TYPE_FOO_1,FOO_LIST); - mem_stat_reg_list((void **)&foo_2,TYPE_FOO_2,FOO_LIST); - mem_stat_dump(stdout,FOO_LIST); - mem_info_file(stdout,FOO_LIST); - mem_stat_free_list(5,FOO_LIST); - mem_stat_dump(stdout,FOO_LIST); - if ( foo_1 != NULL ) - errmesg(" foo_1 is not released"); - if ( foo_2 != NULL ) - errmesg(" foo_2 is not released"); - mem_dump_list(stdout,FOO_LIST); - mem_info_file(stdout,FOO_LIST); - - mem_free_vars(FOO_LIST); - if ( mem_is_list_attached(FOO_LIST) ) - errmesg("list FOO_LIST is not detached"); - - mem_info(); - -#if REAL == FLOAT - printf("# SINGLE PRECISION was used\n"); -#elif REAL == DOUBLE - printf("# DOUBLE PRECISION was used\n"); -#endif - -#define ANSI_OR_VAR - -#ifndef ANSI_C -#ifndef VARARGS -#undef ANSI_OR_VAR -#endif -#endif - -#ifdef ANSI_OR_VAR - - printf("# you should get: \n"); -#if (REAL == FLOAT) - printf("# type VEC: 276 bytes allocated, 3 variables allocated\n"); -#elif (REAL == DOUBLE) - printf("# type VEC: 516 bytes allocated, 3 variables allocated\n"); -#endif - printf("# and other types are zeros\n"); - -#endif /*#if defined(ANSI_C) || defined(VARAGS) */ - - printf("# Finished memory torture test\n"); - return; -} //GO.SYSIN DD memtort.c echo itertort.c 1>&2 sed >itertort.c <<'//GO.SYSIN DD itertort.c' 's/^-//' - -/************************************************************************** -** -** 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. -** -***************************************************************************/ - - -/* iter_tort.c 16/09/93 */ - -/* - This file contains tests for the iterative part of Meschach -*/ - -#include -#include "matrix2.h" -#include "sparse2.h" -#include "iter.h" -#include - -#define errmesg(mesg) printf("Error: %s error: line %d\n",mesg,__LINE__) -#define notice(mesg) printf("# Testing %s...\n",mesg); - - /* for iterative methods */ - -#if REAL == DOUBLE -#define EPS 1e-7 -#define KK 20 -#elif REAL == FLOAT -#define EPS 1e-5 -#define KK 8 -#endif - -#define ANON 513 -#define ASYM ANON - - -static VEC *ex_sol = VNULL; - -/* new iter information */ -void iter_mod_info(ip,nres,res,Bres) -ITER *ip; -double nres; -VEC *res, *Bres; -{ - static VEC *tmp; - - if (ip->b == VNULL) return; - tmp = v_resize(tmp,ip->b->dim); - MEM_STAT_REG(tmp,TYPE_VEC); - - if (nres >= 0.0) { - printf(" %d. residual = %g\n",ip->steps,nres); - } - else - printf(" %d. residual = %g (WARNING !!! should be >= 0) \n", - ip->steps,nres); - if (ex_sol != VNULL) - printf(" ||u_ex - u_approx||_2 = %g\n", - v_norm2(v_sub(ip->x,ex_sol,tmp))); -} - - -/* out = A^T*A*x */ -VEC *norm_equ(A,x,out) -SPMAT *A; -VEC *x, *out; -{ - static VEC * tmp; - - tmp = v_resize(tmp,x->dim); - MEM_STAT_REG(tmp,TYPE_VEC); - sp_mv_mlt(A,x,tmp); - sp_vm_mlt(A,tmp,out); - return out; - -} - - -/* - make symmetric preconditioner for nonsymmetric matrix A; - B = 0.5*(A+A^T) and then B is factorized using - incomplete Choleski factorization -*/ - -SPMAT *gen_sym_precond(A) -SPMAT *A; -{ - SPMAT *B; - SPROW *row; - int i,j,k; - Real val; - - B = sp_get(A->m,A->n,A->row[0].maxlen); - for (i=0; i < A->m; i++) { - row = &(A->row[i]); - for (j = 0; j < row->len; j++) { - k = row->elt[j].col; - if (i != k) { - val = 0.5*(sp_get_val(A,i,k) + sp_get_val(A,k,i)); - sp_set_val(B,i,k,val); - sp_set_val(B,k,i,val); - } - else { /* i == k */ - val = sp_get_val(A,i,i); - sp_set_val(B,i,i,val); - } - } - } - - spICHfactor(B); - return B; -} - -/* Dv_mlt -- diagonal by vector multiply; the diagonal matrix is represented - by a vector d */ -VEC *Dv_mlt(d, x, out) -VEC *d, *x, *out; -{ - int i; - - if ( ! d || ! x ) - error(E_NULL,"Dv_mlt"); - if ( d->dim != x->dim ) - error(E_SIZES,"Dv_mlt"); - out = v_resize(out,x->dim); - - for ( i = 0; i < x->dim; i++ ) - out->ve[i] = d->ve[i]*x->ve[i]; - - return out; -} - - - -/************************************************/ -void main(argc, argv) -int argc; -char *argv[]; -{ - VEC *x, *y, *z, *u, *v, *xn, *yn; - SPMAT *A = NULL, *B = NULL; - SPMAT *An = NULL, *Bn = NULL; - int i, k, kk, j; - ITER *ips, *ips1, *ipns, *ipns1; - MAT *Q, *H, *Q1, *H1; - VEC vt, vt1; - Real hh; - - - mem_info_on(TRUE); - notice("allocating sparse matrices"); - - printf(" dim of A = %dx%d\n",ASYM,ASYM); - - A = iter_gen_sym(ASYM,8); - B = sp_copy(A); - spICHfactor(B); - - u = v_get(A->n); - x = v_get(A->n); - y = v_get(A->n); - v = v_get(A->n); - - v_rand(x); - sp_mv_mlt(A,x,y); - ex_sol = x; - - notice(" initialize ITER variables"); - /* ips for symmetric matrices with precondition */ - ips = iter_get(A->m,A->n); - - /* printf(" ips:\n"); - iter_dump(stdout,ips); */ - - ips->limit = 500; - ips->eps = EPS; - - iter_Ax(ips,sp_mv_mlt,A); - iter_Bx(ips,spCHsolve,B); - - ips->b = v_copy(y,ips->b); - v_rand(ips->x); - /* test of iter_resize */ - ips = iter_resize(ips,2*A->m,2*A->n); - ips = iter_resize(ips,A->m,A->n); - - /* printf(" ips:\n"); - iter_dump(stdout,ips); */ - - /* ips1 for symmetric matrices without precondition */ - ips1 = iter_get(0,0); - /* printf(" ips1:\n"); - iter_dump(stdout,ips1); */ - ITER_FREE(ips1); - - ips1 = iter_copy2(ips,ips1); - iter_Bx(ips1,NULL,NULL); - ips1->b = ips->b; - ips1->shared_b = TRUE; - /* printf(" ips1:\n"); - iter_dump(stdout,ips1); */ - - /* ipns for nonsymetric matrices with precondition */ - ipns = iter_copy(ips,INULL); - ipns->k = KK; - ipns->limit = 500; - ipns->info = NULL; - - An = iter_gen_nonsym_posdef(ANON,8); - Bn = gen_sym_precond(An); - xn = v_get(An->n); - yn = v_get(An->n); - v_rand(xn); - sp_mv_mlt(An,xn,yn); - ipns->b = v_copy(yn,ipns->b); - - iter_Ax(ipns, sp_mv_mlt,An); - iter_ATx(ipns, sp_vm_mlt,An); - iter_Bx(ipns, spCHsolve,Bn); - - /* printf(" ipns:\n"); - iter_dump(stdout,ipns); */ - - /* ipns1 for nonsymmetric matrices without precondition */ - ipns1 = iter_copy2(ipns,INULL); - ipns1->b = ipns->b; - ipns1->shared_b = TRUE; - iter_Bx(ipns1,NULL,NULL); - - /* printf(" ipns1:\n"); - iter_dump(stdout,ipns1); */ - - - /******* CG ********/ - - notice(" CG method without preconditioning"); - ips1->info = NULL; - mem_stat_mark(1); - iter_cg(ips1); - - k = ips1->steps; - z = ips1->x; - printf(" cg: no. of iter.steps = %d\n",k); - v_sub(z,x,u); - printf(" (cg:) ||u_ex - u_approx||_2 = %g [EPS = %g]\n", - v_norm2(u),EPS); - - notice(" CG method with ICH preconditioning"); - - ips->info = NULL; - v_zero(ips->x); - iter_cg(ips); - - k = ips->steps; - printf(" cg: no. of iter.steps = %d\n",k); - v_sub(ips->x,x,u); - printf(" (cg:) ||u_ex - u_approx||_2 = %g [EPS = %g]\n", - v_norm2(u),EPS); - - V_FREE(v); - if ((v = iter_spcg(A,B,y,EPS,VNULL,1000,&k)) == VNULL) - errmesg("CG method with precond.: NULL solution"); - - v_sub(ips->x,v,u); - if (v_norm2(u) >= EPS) { - errmesg("CG method with precond.: different solutions"); - printf(" diff. = %g\n",v_norm2(u)); - } - - - mem_stat_free(1); - printf(" spcg: # of iter. steps = %d\n",k); - v_sub(v,x,u); - printf(" (spcg:) ||u_ex - u_approx||_2 = %g [EPS = %g]\n", - v_norm2(u),EPS); - - - /*** CG FOR NORMAL EQUATION *****/ - - notice("CGNE method with ICH preconditioning (nonsymmetric case)"); - - /* ipns->info = iter_std_info; */ - ipns->info = NULL; - v_zero(ipns->x); - - mem_stat_mark(1); - iter_cgne(ipns); - - k = ipns->steps; - z = ipns->x; - printf(" cgne: no. of iter.steps = %d\n",k); - v_sub(z,xn,u); - printf(" (cgne:) ||u_ex - u_approx||_2 = %g [EPS = %g]\n", - v_norm2(u),EPS); - - notice("CGNE method without preconditioning (nonsymmetric case)"); - - v_rand(u); - u = iter_spcgne(An,NULL,yn,EPS,u,1000,&k); - - mem_stat_free(1); - printf(" spcgne: no. of iter.steps = %d\n",k); - v_sub(u,xn,u); - printf(" (spcgne:) ||u_ex - u_approx||_2 = %g [EPS = %g]\n", - v_norm2(u),EPS); - - /*** CGS *****/ - - notice("CGS method with ICH preconditioning (nonsymmetric case)"); - - v_zero(ipns->x); /* new init guess == 0 */ - - mem_stat_mark(1); - ipns->info = NULL; - v_rand(u); - iter_cgs(ipns,u); - - k = ipns->steps; - z = ipns->x; - printf(" cgs: no. of iter.steps = %d\n",k); - v_sub(z,xn,u); - printf(" (cgs:) ||u_ex - u_approx||_2 = %g [EPS = %g]\n", - v_norm2(u),EPS); - - notice("CGS method without preconditioning (nonsymmetric case)"); - - v_rand(u); - v_rand(v); - v = iter_spcgs(An,NULL,yn,u,EPS,v,1000,&k); - - mem_stat_free(1); - printf(" cgs: no. of iter.steps = %d\n",k); - v_sub(v,xn,u); - printf(" (cgs:) ||u_ex - u_approx||_2 = %g [EPS = %g]\n", - v_norm2(u),EPS); - - - - /*** LSQR ***/ - - notice("LSQR method (without preconditioning)"); - - v_rand(u); - v_free(ipns1->x); - ipns1->x = u; - ipns1->shared_x = TRUE; - ipns1->info = NULL; - mem_stat_mark(2); - z = iter_lsqr(ipns1); - - v_sub(xn,z,v); - k = ipns1->steps; - printf(" lsqr: # of iter. steps = %d\n",k); - printf(" (lsqr:) ||u_ex - u_approx||_2 = %g [EPS = %g]\n", - v_norm2(v),EPS); - - v_rand(u); - u = iter_splsqr(An,yn,EPS,u,1000,&k); - mem_stat_free(2); - - v_sub(xn,u,v); - printf(" splsqr: # of iter. steps = %d\n",k); - printf(" (splsqr:) ||u_ex - u_approx||_2 = %g [EPS = %g]\n", - v_norm2(v),EPS); - - - - /***** GMRES ********/ - - notice("GMRES method with ICH preconditioning (nonsymmetric case)"); - - v_zero(ipns->x); -/* ipns->info = iter_std_info; */ - ipns->info = NULL; - - mem_stat_mark(2); - z = iter_gmres(ipns); - v_sub(xn,z,v); - k = ipns->steps; - printf(" gmres: # of iter. steps = %d\n",k); - printf(" (gmres:) ||u_ex - u_approx||_2 = %g [EPS = %g]\n", - v_norm2(v),EPS); - - notice("GMRES method without preconditioning (nonsymmetric case)"); - V_FREE(v); - v = iter_spgmres(An,NULL,yn,EPS,VNULL,10,1004,&k); - mem_stat_free(2); - - v_sub(xn,v,v); - printf(" spgmres: # of iter. steps = %d\n",k); - printf(" (spgmres:) ||u_ex - u_approx||_2 = %g [EPS = %g]\n", - v_norm2(v),EPS); - - - - /**** MGCR *****/ - - notice("MGCR method with ICH preconditioning (nonsymmetric case)"); - - v_zero(ipns->x); - mem_stat_mark(2); - z = iter_mgcr(ipns); - v_sub(xn,z,v); - k = ipns->steps; - printf(" mgcr: # of iter. steps = %d\n",k); - printf(" (mgcr:) ||u_ex - u_approx||_2 = %g [EPS = %g]\n", - v_norm2(v),EPS); - - notice("MGCR method without preconditioning (nonsymmetric case)"); - V_FREE(v); - v = iter_spmgcr(An,NULL,yn,EPS,VNULL,10,1004,&k); - mem_stat_free(2); - - v_sub(xn,v,v); - printf(" spmgcr: # of iter. steps = %d\n",k); - printf(" (spmgcr:) ||u_ex - u_approx||_2 = %g [EPS = %g]\n", - v_norm2(v),EPS); - - - /***** ARNOLDI METHOD ********/ - - - notice("arnoldi method"); - - kk = ipns1->k = KK; - Q = m_get(kk,x->dim); - Q1 = m_get(kk,x->dim); - H = m_get(kk,kk); - v_rand(u); - ipns1->x = u; - ipns1->shared_x = TRUE; - mem_stat_mark(3); - iter_arnoldi_iref(ipns1,&hh,Q,H); - mem_stat_free(3); - - /* check the equality: - Q*A*Q^T = H; */ - - vt.dim = vt.max_dim = x->dim; - vt1.dim = vt1.max_dim = x->dim; - for (j=0; j < kk; j++) { - vt.ve = Q->me[j]; - vt1.ve = Q1->me[j]; - sp_mv_mlt(An,&vt,&vt1); - } - H1 = m_get(kk,kk); - mmtr_mlt(Q,Q1,H1); - m_sub(H,H1,H1); - if (m_norm_inf(H1) > MACHEPS*x->dim) - printf(" (arnoldi_iref) ||Q*A*Q^T - H|| = %g [cf. MACHEPS = %g]\n", - m_norm_inf(H1),MACHEPS); - - /* check Q*Q^T = I */ - - mmtr_mlt(Q,Q,H1); - for (j=0; j < kk; j++) - H1->me[j][j] -= 1.0; - if (m_norm_inf(H1) > MACHEPS*x->dim) - printf(" (arnoldi_iref) ||Q*Q^T - I|| = %g [cf. MACHEPS = %g]\n", - m_norm_inf(H1),MACHEPS); - - ipns1->x = u; - ipns1->shared_x = TRUE; - mem_stat_mark(3); - iter_arnoldi(ipns1,&hh,Q,H); - mem_stat_free(3); - - /* check the equality: - Q*A*Q^T = H; */ - - vt.dim = vt.max_dim = x->dim; - vt1.dim = vt1.max_dim = x->dim; - for (j=0; j < kk; j++) { - vt.ve = Q->me[j]; - vt1.ve = Q1->me[j]; - sp_mv_mlt(An,&vt,&vt1); - } - - mmtr_mlt(Q,Q1,H1); bigmail CUT HERE............ test -w meschach0.shar && test -r 24048P00 && test -r 24048P01 && test -r 24048P02 && test -r 24048P03 && ( cat 24048P00 >> meschach0.shar; rm 24048P00 cat 24048P01 >> meschach0.shar; rm 24048P01 cat 24048P02 >> meschach0.shar; rm 24048P02 cat 24048P03 >> meschach0.shar; rm 24048P03 )