home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
OS/2 Shareware BBS: 22 gnu
/
22-gnu.zip
/
mesch12a.zip
/
itertort.c
< prev
next >
Wrap
C/C++ Source or Header
|
1994-01-14
|
17KB
|
692 lines
/**************************************************************************
**
** 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 <stdio.h>
#include <math.h>
#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
#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);
m_sub(H,H1,H1);
if (m_norm_inf(H1) > MACHEPS*x->dim)
printf(" (arnoldi) ||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) ||Q*Q^T - I|| = %g [cf. MACHEPS = %g]\n",
m_norm_inf(H1),MACHEPS);
v_rand(u);
mem_stat_mark(3);
iter_sparnoldi(An,u,kk,&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);
m_sub(H,H1,H1);
if (m_norm_inf(H1) > MACHEPS*x->dim)
printf(" (sparnoldi) ||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(" (sparnoldi) ||Q*Q^T - I|| = %g [cf. MACHEPS = %g]\n",
m_norm_inf(H1),MACHEPS);
/****** LANCZOS METHOD ******/
notice("lanczos method");
kk = ipns1->k;
Q = m_resize(Q,kk,x->dim);
Q1 = m_resize(Q1,kk,x->dim);
H = m_resize(H,kk,kk);
ips1->k = kk;
v_rand(u);
v_free(ips1->x);
ips1->x = u;
ips1->shared_x = TRUE;
mem_stat_mark(3);
iter_lanczos(ips1,x,y,&hh,Q);
mem_stat_free(3);
/* check the equality:
Q*A*Q^T = H; */
vt.dim = vt1.dim = Q->n;
vt.max_dim = vt1.max_dim = Q->max_n;
Q1 = m_resize(Q1,Q->m,Q->n);
for (j=0; j < Q->m; j++) {
vt.ve = Q->me[j];
vt1.ve = Q1->me[j];
sp_mv_mlt(A,&vt,&vt1);
}
H1 = m_resize(H1,Q->m,Q->m);
H = m_resize(H,Q->m,Q->m);
mmtr_mlt(Q,Q1,H1);
m_zero(H);
for (j=0; j < Q->m-1; j++) {
H->me[j][j] = x->ve[j];
H->me[j][j+1] = H->me[j+1][j] = y->ve[j];
}
H->me[Q->m-1][Q->m-1] = x->ve[Q->m-1];
m_sub(H,H1,H1);
if (m_norm_inf(H1) > MACHEPS*x->dim)
printf(" (lanczos) ||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 < Q->m; j++)
H1->me[j][j] -= 1.0;
if (m_norm_inf(H1) > MACHEPS*x->dim)
printf(" (lanczos) ||Q*Q^T - I|| = %g [cf. MACHEPS = %g]\n",
m_norm_inf(H1),MACHEPS);
mem_stat_mark(3);
v_rand(u);
iter_splanczos(A,kk,u,x,y,&hh,Q);
mem_stat_free(3);
/* check the equality:
Q*A*Q^T = H; */
vt.dim = vt1.dim = Q->n;
vt.max_dim = vt1.max_dim = Q->max_n;
Q1 = m_resize(Q1,Q->m,Q->n);
for (j=0; j < Q->m; j++) {
vt.ve = Q->me[j];
vt1.ve = Q1->me[j];
sp_mv_mlt(A,&vt,&vt1);
}
H1 = m_resize(H1,Q->m,Q->m);
H = m_resize(H,Q->m,Q->m);
mmtr_mlt(Q,Q1,H1);
for (j=0; j < Q->m-1; j++) {
H->me[j][j] = x->ve[j];
H->me[j][j+1] = H->me[j+1][j] = y->ve[j];
}
H->me[Q->m-1][Q->m-1] = x->ve[Q->m-1];
m_sub(H,H1,H1);
if (m_norm_inf(H1) > MACHEPS*x->dim)
printf(" (splanczos) ||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 < Q->m; j++)
H1->me[j][j] -= 1.0;
if (m_norm_inf(H1) > MACHEPS*x->dim)
printf(" (splanczos) ||Q*Q^T - I|| = %g [cf. MACHEPS = %g]\n",
m_norm_inf(H1),MACHEPS);
/***** LANCZOS2 ****/
notice("lanczos2 method");
kk = 50; /* # of dir. vectors */
ips1->k = kk;
v_rand(u);
ips1->x = u;
ips1->shared_x = TRUE;
for ( i = 0; i < xn->dim; i++ )
xn->ve[i] = i;
iter_Ax(ips1,Dv_mlt,xn);
mem_stat_mark(3);
iter_lanczos2(ips1,y,v);
mem_stat_free(3);
printf("# Number of steps of Lanczos algorithm = %d\n", kk);
printf("# Exact eigenvalues are 0, 1, 2, ..., %d\n",ANON-1);
printf("# Extreme eigenvalues should be accurate; \n");
printf("# interior values usually are not.\n");
printf("# approx e-vals =\n"); v_output(y);
printf("# Error in estimate of bottom e-vec (Lanczos) = %g\n",
fabs(v->ve[0]));
mem_stat_mark(3);
v_rand(u);
iter_splanczos2(A,kk,u,y,v);
mem_stat_free(3);
/***** FINISHING *******/
notice("release ITER variables");
M_FREE(Q);
M_FREE(Q1);
M_FREE(H);
M_FREE(H1);
ITER_FREE(ipns);
ITER_FREE(ips);
ITER_FREE(ipns1);
ITER_FREE(ips1);
SP_FREE(A);
SP_FREE(B);
SP_FREE(An);
SP_FREE(Bn);
V_FREE(x);
V_FREE(y);
V_FREE(u);
V_FREE(v);
V_FREE(xn);
V_FREE(yn);
printf("# Done testing (%s)\n",argv[0]);
mem_info();
}