#include <stdio.h>
#include "matrix.h"
#include "matrix2.h"
#include "sparse.h"
#include "iter.h"
#include <math.h>
static
char
rcsid[] =
"$Id: itersym.c,v 1.2 1995/01/30 14:55:54 des Exp $"
;
#ifdef ANSI_C
VEC *spCHsolve(SPMAT *,VEC *,VEC *);
VEC *trieig(VEC *,VEC *,MAT *);
#else
VEC *spCHsolve();
VEC *trieig();
#endif
VEC *iter_spcg(A,LLT,b,eps,x,limit,steps)
SPMAT *A, *LLT;
VEC *b, *x;
double
eps;
int
*steps, limit;
{
ITER *ip;
ip = iter_get(0,0);
ip->Ax = (Fun_Ax) sp_mv_mlt;
ip->A_par = (
void
*)A;
ip->Bx = (Fun_Ax) spCHsolve;
ip->B_par = (
void
*)LLT;
ip->info = (Fun_info) NULL;
ip->b = b;
ip->eps = eps;
ip->limit = limit;
ip->x = x;
iter_cg(ip);
x = ip->x;
if
(steps) *steps = ip->steps;
ip->shared_x = ip->shared_b = TRUE;
iter_free(ip);
return
x;
}
VEC *iter_cg(ip)
ITER *ip;
{
static
VEC *r = VNULL, *p = VNULL, *q = VNULL, *z = VNULL;
Real alpha, beta, inner, old_inner, nres;
VEC *rr;
if
(ip == INULL)
error(E_NULL,
"iter_cg"
);
if
(!ip->Ax || !ip->b)
error(E_NULL,
"iter_cg"
);
if
( ip->x == ip->b )
error(E_INSITU,
"iter_cg"
);
if
(!ip->stop_crit)
error(E_NULL,
"iter_cg"
);
if
( ip->eps <= 0.0 )
ip->eps = MACHEPS;
r = v_resize(r,ip->b->dim);
p = v_resize(p,ip->b->dim);
q = v_resize(q,ip->b->dim);
MEM_STAT_REG(r,TYPE_VEC);
MEM_STAT_REG(p,TYPE_VEC);
MEM_STAT_REG(q,TYPE_VEC);
if
(ip->Bx != (Fun_Ax)NULL) {
z = v_resize(z,ip->b->dim);
MEM_STAT_REG(z,TYPE_VEC);
rr = z;
}
else
rr = r;
if
(ip->x != VNULL) {
if
(ip->x->dim != ip->b->dim)
error(E_SIZES,
"iter_cg"
);
ip->Ax(ip->A_par,ip->x,p);
v_sub(ip->b,p,r);
}
else
{
ip->x = v_get(ip->b->dim);
ip->shared_x = FALSE;
v_copy(ip->b,r);
}
old_inner = 0.0;
for
( ip->steps = 0; ip->steps <= ip->limit; ip->steps++ )
{
if
( ip->Bx )
(ip->Bx)(ip->B_par,r,rr);
inner = in_prod(rr,r);
nres =
sqrt
(
fabs
(inner));
if
(ip->info) ip->info(ip,nres,r,rr);
if
(ip->steps == 0) ip->init_res = nres;
if
( ip->stop_crit(ip,nres,r,rr) )
break
;
if
( ip->steps )
{
beta = inner/old_inner;
p = v_mltadd(rr,p,beta,p);
}
else
{
beta = 0.0;
p = v_copy(rr,p);
old_inner = 0.0;
}
(ip->Ax)(ip->A_par,p,q);
alpha = in_prod(p,q);
if
(
sqrt
(
fabs
(alpha)) <= MACHEPS*ip->init_res)
error(E_BREAKDOWN,
"iter_cg"
);
alpha = inner/alpha;
v_mltadd(ip->x,p,alpha,ip->x);
v_mltadd(r,q,-alpha,r);
old_inner = inner;
}
return
ip->x;
}
void
iter_lanczos(ip,a,b,beta2,Q)
ITER *ip;
VEC *a, *b;
Real *beta2;
MAT *Q;
{
int
j;
static
VEC *v = VNULL, *w = VNULL, *tmp = VNULL;
Real alpha, beta, c;
if
( ! ip )
error(E_NULL,
"iter_lanczos"
);
if
( ! ip->Ax || ! ip->x || ! a || ! b )
error(E_NULL,
"iter_lanczos"
);
if
( ip->k <= 0 )
error(E_BOUNDS,
"iter_lanczos"
);
if
( Q && ( Q->n < ip->x->dim || Q->m < ip->k ) )
error(E_SIZES,
"iter_lanczos"
);
a = v_resize(a,(u_int)ip->k);
b = v_resize(b,(u_int)(ip->k-1));
v = v_resize(v,ip->x->dim);
w = v_resize(w,ip->x->dim);
tmp = v_resize(tmp,ip->x->dim);
MEM_STAT_REG(v,TYPE_VEC);
MEM_STAT_REG(w,TYPE_VEC);
MEM_STAT_REG(tmp,TYPE_VEC);
beta = 1.0;
v_zero(a);
v_zero(b);
if
(Q) m_zero(Q);
c = v_norm2(ip->x);
if
(c <= MACHEPS) {
*beta2 = 0.0;
return
;
}
else
sv_mlt(1.0/c,ip->x,w);
(ip->Ax)(ip->A_par,w,v);
for
( j = 0; j < ip->k; j++ )
{
if
( Q ) set_row(Q,j,w);
alpha = in_prod(w,v);
a->ve[j] = alpha;
v_mltadd(v,w,-alpha,v);
beta = v_norm2(v);
if
( beta == 0.0 )
{
*beta2 = 0.0;
return
;
}
if
( j < ip->k-1 )
b->ve[j] = beta;
v_copy(w,tmp);
sv_mlt(1/beta,v,w);
sv_mlt(-beta,tmp,v);
(ip->Ax)(ip->A_par,w,tmp);
v_add(v,tmp,v);
}
*beta2 = beta;
}
void
iter_splanczos(A,m,x0,a,b,beta2,Q)
SPMAT *A;
int
m;
VEC *x0, *a, *b;
Real *beta2;
MAT *Q;
{
ITER *ip;
ip = iter_get(0,0);
ip->shared_x = ip->shared_b = TRUE;
ip->Ax = (Fun_Ax) sp_mv_mlt;
ip->A_par = (
void
*) A;
ip->x = x0;
ip->k = m;
iter_lanczos(ip,a,b,beta2,Q);
iter_free(ip);
}
extern
double
frexp
(),
ldexp
();
static
double
product(a,offset,expt)
VEC *a;
double
offset;
int
*expt;
{
Real mant, tmp_fctr;
int
i, tmp_expt;
if
( ! a )
error(E_NULL,
"product"
);
mant = 1.0;
*expt = 0;
if
( offset == 0.0 )
for
( i = 0; i < a->dim; i++ )
{
mant *=
frexp
(a->ve[i],&tmp_expt);
*expt += tmp_expt;
if
( ! (i % 10) )
{
mant =
frexp
(mant,&tmp_expt);
*expt += tmp_expt;
}
}
else
for
( i = 0; i < a->dim; i++ )
{
tmp_fctr = a->ve[i] - offset;
tmp_fctr += (tmp_fctr > 0.0 ) ? -MACHEPS*offset :
MACHEPS*offset;
mant *=
frexp
(tmp_fctr,&tmp_expt);
*expt += tmp_expt;
if
( ! (i % 10) )
{
mant =
frexp
(mant,&tmp_expt);
*expt += tmp_expt;
}
}
mant =
frexp
(mant,&tmp_expt);
*expt += tmp_expt;
return
mant;
}
static
double
product2(a,k,expt)
VEC *a;
int
k;
int
*expt;
{
Real mant, mu, tmp_fctr;
int
i, tmp_expt;
if
( ! a )
error(E_NULL,
"product2"
);
if
( k < 0 || k >= a->dim )
error(E_BOUNDS,
"product2"
);
mant = 1.0;
*expt = 0;
mu = a->ve[k];
for
( i = 0; i < a->dim; i++ )
{
if
( i == k )
continue
;
tmp_fctr = a->ve[i] - mu;
tmp_fctr += ( tmp_fctr > 0.0 ) ? -MACHEPS*mu : MACHEPS*mu;
mant *=
frexp
(tmp_fctr,&tmp_expt);
*expt += tmp_expt;
if
( ! (i % 10) )
{
mant =
frexp
(mant,&tmp_expt);
*expt += tmp_expt;
}
}
mant =
frexp
(mant,&tmp_expt);
*expt += tmp_expt;
return
mant;
}
static
int
dbl_cmp(x,y)
Real *x, *y;
{
Real tmp;
tmp = *x - *y;
return
(tmp > 0 ? 1 : tmp < 0 ? -1: 0);
}
VEC *iter_lanczos2(ip,evals,err_est)
ITER *ip;
VEC *evals;
VEC *err_est;
{
VEC *a;
static
VEC *b=VNULL, *a2=VNULL, *b2=VNULL;
Real beta, pb_mant, det_mant, det_mant1, det_mant2;
int
i, pb_expt, det_expt, det_expt1, det_expt2;
if
( ! ip )
error(E_NULL,
"iter_lanczos2"
);
if
( ! ip->Ax || ! ip->x )
error(E_NULL,
"iter_lanczos2"
);
if
( ip->k <= 0 )
error(E_RANGE,
"iter_lanczos2"
);
a = evals;
a = v_resize(a,(u_int)ip->k);
b = v_resize(b,(u_int)(ip->k-1));
MEM_STAT_REG(b,TYPE_VEC);
iter_lanczos(ip,a,b,&beta,MNULL);
pb_mant = 0.0;
if
( err_est )
{
pb_mant = product(b,(
double
)0.0,&pb_expt);
}
a2 = v_resize(a2,a->dim - 1);
b2 = v_resize(b2,b->dim - 1);
MEM_STAT_REG(a2,TYPE_VEC);
MEM_STAT_REG(b2,TYPE_VEC);
for
( i = 0; i < a2->dim - 1; i++ )
{
a2->ve[i] = a->ve[i+1];
b2->ve[i] = b->ve[i+1];
}
a2->ve[a2->dim-1] = a->ve[a2->dim];
trieig(a,b,MNULL);
qsort
((
void
*)(a->ve),(
int
)(a->dim),
sizeof
(Real),(
int
(*)())dbl_cmp);
if
( err_est )
{
err_est = v_resize(err_est,(u_int)ip->k);
trieig(a2,b2,MNULL);
for
( i = 0; i < a->dim; i++ )
{
det_mant1 = product2(a,i,&det_expt1);
det_mant2 = product(a2,(
double
)a->ve[i],&det_expt2);
if
( det_mant1 == 0.0 )
{
err_est->ve[i] = 0.0;
continue
;
}
else
if
( det_mant2 == 0.0 )
{
err_est->ve[i] = HUGE;
continue
;
}
if
( (det_expt1 + det_expt2) % 2 )
det_mant =
sqrt
(2.0*
fabs
(det_mant1*det_mant2));
else
det_mant =
sqrt
(
fabs
(det_mant1*det_mant2));
det_expt = (det_expt1+det_expt2)/2;
err_est->ve[i] =
fabs
(beta*
ldexp
(pb_mant/det_mant,pb_expt-det_expt));
}
}
return
a;
}
VEC *iter_splanczos2(A,m,x0,evals,err_est)
SPMAT *A;
int
m;
VEC *x0;
VEC *evals;
VEC *err_est;
{
ITER *ip;
VEC *a;
ip = iter_get(0,0);
ip->Ax = (Fun_Ax) sp_mv_mlt;
ip->A_par = (
void
*) A;
ip->x = x0;
ip->k = m;
a = iter_lanczos2(ip,evals,err_est);
ip->shared_x = ip->shared_b = TRUE;
iter_free(ip);
return
a;
}
VEC *iter_cg1(ip)
ITER *ip;
{
static
VEC *r = VNULL, *p = VNULL, *q = VNULL, *z = VNULL;
Real alpha;
double
inner,nres;
VEC *rr;
if
(ip == INULL)
error(E_NULL,
"iter_cg"
);
if
(!ip->Ax || !ip->b)
error(E_NULL,
"iter_cg"
);
if
( ip->x == ip->b )
error(E_INSITU,
"iter_cg"
);
if
(!ip->stop_crit)
error(E_NULL,
"iter_cg"
);
if
( ip->eps <= 0.0 )
ip->eps = MACHEPS;
r = v_resize(r,ip->b->dim);
p = v_resize(p,ip->b->dim);
q = v_resize(q,ip->b->dim);
MEM_STAT_REG(r,TYPE_VEC);
MEM_STAT_REG(p,TYPE_VEC);
MEM_STAT_REG(q,TYPE_VEC);
if
(ip->Bx != (Fun_Ax)NULL) {
z = v_resize(z,ip->b->dim);
MEM_STAT_REG(z,TYPE_VEC);
rr = z;
}
else
rr = r;
if
(ip->x != VNULL) {
if
(ip->x->dim != ip->b->dim)
error(E_SIZES,
"iter_cg"
);
ip->Ax(ip->A_par,ip->x,p);
v_sub(ip->b,p,r);
}
else
{
ip->x = v_get(ip->b->dim);
ip->shared_x = FALSE;
v_copy(ip->b,r);
}
if
(ip->Bx) (ip->Bx)(ip->B_par,r,p);
else
v_copy(r,p);
inner = in_prod(p,r);
nres =
sqrt
(
fabs
(inner));
if
(ip->info) ip->info(ip,nres,r,p);
if
( nres == 0.0)
return
ip->x;
for
( ip->steps = 0; ip->steps <= ip->limit; ip->steps++ )
{
ip->Ax(ip->A_par,p,q);
inner = in_prod(q,p);
if
(
sqrt
(
fabs
(inner)) <= MACHEPS*ip->init_res)
error(E_BREAKDOWN,
"iter_cg1"
);
alpha = in_prod(p,r)/inner;
v_mltadd(ip->x,p,alpha,ip->x);
v_mltadd(r,q,-alpha,r);
rr = r;
if
(ip->Bx) {
ip->Bx(ip->B_par,r,z);
rr = z;
}
nres = in_prod(r,rr);
if
(nres < 0.0) {
warning(WARN_RES_LESS_0,
"iter_cg"
);
break
;
}
nres =
sqrt
(
fabs
(nres));
if
(ip->info) ip->info(ip,nres,r,z);
if
(ip->steps == 0) ip->init_res = nres;
if
( ip->stop_crit(ip,nres,r,z) )
break
;
alpha = -in_prod(rr,q)/inner;
v_mltadd(rr,p,alpha,p);
}
return
ip->x;
}