file: Hlib/src/Solve.c
/**************************************************************************/
// //
// Author: T.Warburton //
// Design: T.Warburton && S.Sherwin //
// Date : 12/4/96 //
// //
// Copyright notice: This code shall not be replicated or used without //
// the permission of the author. //
// //
/**************************************************************************/
#include <math.h>
#include <veclib.h>
#include "hotel.h"
/* general utilities routines */
void ScatrBndry(double *u, Element_List *U, Bsystem *B);
static void setupRHS (Element_List *, Element_List *, double *, double *, //OVERLOAD CALL: setupRHS: Solve.c(?), Solve_Stokes.c(?)
Bndry *, Bsystem *, SolveType);
static void saveinit (Element_List *, double *, Bsystem *); //OVERLOAD CALL: saveinit: Solve.c(?), Solve_Stokes.c(?)
void solve_boundary
(Element_List *, Element_List *, double *, double *,Bsystem *);
void Solve(Element_List *U, Element_List *Uf, Bndry *Ubc, Bsystem *Ubsys,
SolveType Stype){
static double *RHS = (double*)0;
static double *U0 = (double*)0;
static int nglob = 0;
if(Ubsys->nglobal > nglob){
if(nglob){
free(U0);
free(RHS);
}
nglob = Ubsys->nglobal;
RHS = dvector(0,nglob-1);
U0 = dvector(0,nglob-1);
}
dzero(Ubsys->nglobal, RHS, 1);
dzero(Ubsys->nglobal, U0, 1);
/* ----------------------------*/
/* Compute the Right-Hand Side */ setupRHS(U,Uf,RHS,U0,Ubc,Ubsys,Stype); //OVERLOAD CALL: setupRHS: Solve.c(?), Solve_Stokes.c(?)
/*-----------------------------*/
/*-----------------------------*/
/* Solve the boundary system */ solve_boundary(U,Uf,RHS,U0,Ubsys);
/*-----------------------------*/
/*-----------------------------*/
/* solve the Interior system */ solve_interior(U,Ubsys);
/*-----------------------------*/
// free(RHS); free(U0);
}
static void saveinit(Element_List *U, double *u0, Bsystem *B){
int **bmap = B->bmap;
Element *E;
SignChange(U,B);
for(E=U->fhead;E;E=E->next)
dscatr(E->Nbmodes,E->vert->hj,bmap[E->id],u0);
SignChange(U,B);
}
double tol;
static double tolv[3];
static int toli=0;
static void Set_tolerance(Element_List *E_L, Element_List *Ef_L){
Element *E = E_L->fhead;
Element *Ef = Ef_L->fhead;
for(tol = 0.0;Ef;Ef = Ef->next)
tol += ddot(Ef->Nbmodes ,Ef->vert->hj,1,Ef->vert->hj,1);
#ifdef PARALLEL
DO_PARALLEL{ /* gather tolerances together to make global tolerance */
double tmp;
gdsum(&tol,1,&tmp);
}
#endif
tol = sqrt(tol);
if(E->type != 'p'){
tolv[toli] = tol;
if(E->dim() == 2){ //OVERLOAD CALL: dim: Pyr.h(?), Quad.h(?), Tet.h(?), Tri.h(?), Hex.h(?), Prism.h(?)
tol = sqrt(tolv[0]*tolv[0] + tolv[1]*tolv[1])/2.;
toli = (++toli)%2;
}
else{
tol = sqrt(tolv[0]*tolv[0] + tolv[1]*tolv[1] + tolv[2]*tolv[2])
/3.;
toli = (++toli)%3;
}
}
/* for those wierd cases make sure tol is not less than 10e-3 */
tol = (tol < 10e-3)? 1.0: tol;
}
/* compute right hand side and put leave in U */
static void setupRHS (Element_List *U, Element_List *Uf,
double *rhs, double *u0,
Bndry *Ubc, Bsystem *B, SolveType Stype){
register int k;
int nel = B->nel;
int N,nbl;
double **binvc = B->Gmat->binvc;
Element *E;
Bndry *Ebc;
/* save initial condition */
saveinit(U,u0,B); //OVERLOAD CALL: saveinit: Solve.c(?), Solve_Stokes.c(?)
if(Uf->fhead->state == 'p') /* take inner product if in physical space */
Uf->Iprod(Uf); //OVERLOAD CALL: Iprod: InnerProd.c(Element), Element_List.c(Element_List), InnerProd.c(Tri), InnerProd.c(Quad), InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex)
// fixed for contiguous memory model
if(Stype == Helm) /* negate if helmholtz solve */
dneg(Uf->hjtot, Uf->base_hj, 1);
if(B->nsolve){ /* subtract off interior coupling */
for(E=Uf->fhead;E;E=E->next){
nbl = E->Nbmodes;
N = E->Nmodes - nbl;
if(N) dgemv('T', N, nbl, -1., binvc[E->geom->id], N,
E->vert->hj+nbl, 1, 1., E->vert->hj,1);
}
}
// needs to be fixed
/* add flux terms */
for(Ebc = Ubc; Ebc; Ebc = Ebc->next)
if(Ebc->type == 'F' || Ebc->type == 'R')
Uf->flist[Ebc->elmt->id]->Add_flux_terms(Ebc); //OVERLOAD CALL: Add_flux_terms: Boundary.c(Tri), Boundary.c(Quad), Boundary.c(Tet), Boundary.c(Pyr), Boundary.c(Prism), Boundary.c(Hex), Boundary.c(Element)
/* subtract off Boundary initial condition */
if(B->smeth == iterative&&!B->rslv){
double **a = B->Gmat->a;
Set_tolerance(U,Uf); /* set solver tolerances */ //OVERLOAD CALL: Set_tolerance: Solve.c(?), Solve_cg.c(?)
for(k = 0; k < nel; ++k){
dspmv('L',Uf->flist[k]->Nbmodes,-1.0,a[U->flist[k]->geom->id],
U->flist[k]->vert->hj,1,1.,Uf->flist[k]->vert->hj,1);
dcopy(Uf->flist[k]->Nmodes,Uf->flist[k]->vert->hj,1, U->flist[k]->vert->hj,1);
U->flist[k]->state = 't';
}
GathrBndry(U,rhs,B);
dzero(B->nglobal-B->nsolve, rhs + B->nsolve, 1);
}
else{
if(Ubc&&Ubc->DirRHS){
/* for non pressure solve have preset vector */ //OVERLOAD CALL: vector: Coords.c(?), Curvi.c(?)
/* note that this implies that the bc's are time independent */
/* or can be expressed as B(x)F(t) where F(t) is given */
// fixed for contiguous memory model
dcopy(Uf->hjtot, Uf->base_hj, 1, U->base_hj, 1);
U->Set_state('t');
GathrBndry(U,rhs,B);
dzero(B->nglobal-B->nsolve, rhs + B->nsolve, 1);
/* subtract of bcs */
dvsub(B->nsolve,rhs,1,Ubc->DirRHS,1,rhs,1);
/* zero ic vector */ //OVERLOAD CALL: vector: Coords.c(?), Curvi.c(?)
dzero(B->nsolve,u0,1);
}
else{
for(k = 0; k < nel; ++k)
dzero(U->flist[k]->Nmodes-U->flist[k]->Nbmodes,
U->flist[k]->vert->hj + U->flist[k]->Nbmodes, 1);
if(Stype == Helm){
U->Trans(U, J_to_Q); //OVERLOAD CALL: Trans: Element_List.c(Element_List), Transforms.c(Element), Transforms.c(Tri), Transforms.c(Quad), Transforms.c(Tet), Transforms.c(Pyr), Transforms.c(Prism), Transforms.c(Hex)
U->HelmHoltz(B->lambda); //OVERLOAD CALL: HelmHoltz: InnerProd.c(Quad), InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex), InnerProd.c(Element), Element_List.c(Element_List), InnerProd.c(Tri)
}
else{
U->Trans(U, J_to_Q); //OVERLOAD CALL: Trans: Element_List.c(Element_List), Transforms.c(Element), Transforms.c(Tri), Transforms.c(Quad), Transforms.c(Tet), Transforms.c(Pyr), Transforms.c(Prism), Transforms.c(Hex)
U->Iprod(U); //OVERLOAD CALL: Iprod: InnerProd.c(Element), Element_List.c(Element_List), InnerProd.c(Tri), InnerProd.c(Quad), InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex)
}
for(k = 0; k < nel; ++k){
nbl = U->flist[k]->Nbmodes;
N = U->flist[k]->Nmodes - nbl;
if(N) dgemv('T',N, nbl, -1., binvc[U->flist[k]->geom->id], N,
U->flist[k]->vert->hj+nbl, 1, 1.0,
U->flist[k]->vert->hj,1);
dvsub(nbl,Uf->flist[k]->vert->hj,1, U->flist[k]->vert->hj, 1,
U->flist[k]->vert->hj, 1);
dcopy(N ,Uf->flist[k]->vert->hj+nbl,1, U->flist[k]->vert->hj+nbl, 1);
}
GathrBndry(U,rhs,B);
dzero(B->nglobal-B->nsolve, rhs + B->nsolve, 1);
}
}
}
static void Bsolve_CG (Element_List *, Element_List *, Bsystem *, double *); //OVERLOAD CALL: Bsolve_CG: Solve.c(?), Solve.c(?)
void solve_boundary
(Element_List *U, Element_List *Uf, double *rhs, double *u0, Bsystem *B){
const int nsolve = B->nsolve;
int info;
if(nsolve){
if(B->rslv){ /* recursive Static condensation solver */
Rsolver *R = B->rslv;
int nrecur = R->nrecur;
int aslv = R->rdata[nrecur-1].cstart, bw = R->Ainfo.bwidth_a;
Recur_setrhs(R,rhs);
if(B->smeth == direct){
if(B->singular) rhs[B->singular-1] = 0.0;
if(aslv)
if(2*bw < aslv) /* banded matrix */
dpbtrs('L', aslv, bw-1, 1, R->A.inva, bw, rhs, aslv, info);
else /* symmatric matrix */
dpptrs('L', aslv, 1, R->A.inva, rhs, aslv, info);
}
else
Recur_Bsolve_CG(B,rhs,U->flist[0]->type);
#if 0
Recur_backslv(R,rhs); //OVERLOAD CALL: Recur_backslv: SolveR.c(?), SolveR.c(?)
#endif
Recur_backslv(R,rhs,'p'); //OVERLOAD CALL: Recur_backslv: SolveR.c(?), SolveR.c(?)
}
else{
const int bwidth = B->Gmat->bwidth_a;
if(B->smeth == iterative)
Bsolve_CG(U, Uf, B, rhs); //OVERLOAD CALL: Bsolve_CG: Solve.c(?), Solve.c(?)
else{
if(B->singular) rhs[B->singular-1] = 0.0;
if(2*bwidth < nsolve) /* banded matrix */
dpbtrs('L', nsolve, bwidth-1, 1, *B->Gmat->inva, bwidth,
rhs, nsolve, info);
else /* symmatric matrix */
dpptrs('L', nsolve, 1, *B->Gmat->inva, rhs, nsolve, info);
}
}
}
/* add intial conditions */
dvadd(B->nglobal,u0,1,rhs,1,rhs,1);
ScatrBndry(rhs,U,B);
}
void solve_interior(Element_List *U, Bsystem *Ubsys){
int N,nbl,id,info;
int *bw = Ubsys->Gmat->bwidth_c;
double *hj;
double **invc = Ubsys->Gmat-> invc;
double **binvc = Ubsys->Gmat->binvc;
Element *E;
for(E=U->fhead;E;E=E->next)
{
N = E->Nmodes - E->Nbmodes;
if(!N) continue;
id = E->geom->id;
nbl = E->Nbmodes;
hj = E->vert->hj + nbl;
if(N > 2*bw[id])
dpbtrs('L', N, bw[id]-1, 1, invc[id], bw[id], hj, N, info);
else
dpptrs('L', N, 1, invc[id], hj, N, info);
dgemv('N', N, nbl,-1., binvc[id], N, E->vert->hj, 1, 1.,hj,1);
}
}
#if 0
void SetBCs(Element_List *U_L, Bndry *Ubc, Bsystem *Ubsys){
register int i;
double *bc = dvector(Ubsys->nsolve,Ubsys->nglobal);
Element *U;
dzero((Ubsys->nglobal-Ubsys->nsolve),bc+Ubsys->nsolve,1);
for(;Ubc;Ubc= Ubc->next)
U_L->flist[Ubc->elmt->id]->setbcs(Ubc,bc); //OVERLOAD CALL: setbcs: Boundary.c(Tri), Boundary.c(Quad), Boundary.c(Tet), Boundary.c(Pyr), Boundary.c(Prism), Boundary.c(Hex), Boundary.c(Element)
/* need to do this so that all vertices get boundary condition */
for(U=U_L->fhead;U;U = U->next)
for(i = 0; i < U->Nverts; ++i)
if(!U->vert[i].solve)
U->vert[i].hj[0] = bc[U->vert[i].gid];
free(bc+Ubsys->nsolve);
}
#else
void SetBCs(Element_List *EL, Bndry *Ubc, Bsystem *Ubsys){
if(!Ubsys->signchange)
setup_signchange(EL,Ubsys);
register int i;
int eid,face,nsolve,nglobal,skip,gid,eDIM,ll;
Vert *v;
Edge *e;
double *bc,scal;
register int j;
int nfv;
Face *f;
Element *E;
eDIM = EL->fhead->dim(); //OVERLOAD CALL: dim: Pyr.h(?), Quad.h(?), Tet.h(?), Tri.h(?), Hex.h(?), Prism.h(?)
nsolve = Ubsys->nsolve;
nglobal = Ubsys->nglobal;
#ifdef PARALLEL
DO_PARALLEL{
if(Ubsys->pll->nsolve == Ubsys->pll->nglobal)
return;
}
else
#endif
if(!(nglobal-nsolve)) return;
bc = dvector(0,nglobal-nsolve);
dzero(nglobal-nsolve, bc, 1);
#if 0
double *mult;
if(option("E")){
mult = dvector(0,nglobal-nsolve);
dzero(nglobal-nsolve, mult, 1);
for(;Ubc;Ubc= Ubc->next){
switch(Ubc->type){
case 'W': case 'V': case 'v': case 'o': case 'm': case 'M':
eid = Ubc->elmt->id;
face = Ubc->face;
E = EL->flist[eid];
nfv = E->Nfverts(face); //OVERLOAD CALL: Nfverts: Nums.c(Tri), Nums.c(Quad), Nums.c(Tet), Nums.c(Pyr), Nums.c(Prism), Nums.c(Hex), Nums.c(Element)
for(i=0;i<nfv;++i)
mult[E->vert[E->vnum(face,i)].gid-nsolve] += 1.; //OVERLOAD CALL: vnum: Nums.c(Quad), Nums.c(Tet), Nums.c(Pyr), Nums.c(Prism), Nums.c(Hex), Nums.c(Element), Nums.c(Tri)
}
}
}
#endif
/* copy all Dirichlet values from Ubc into bc */
for(;Ubc;Ubc= Ubc->next){
switch(Ubc->type){
case 'W': case 'V': case 'v': case 'o': case 'm': case 'M':
scal = 1.0;
eid = Ubc->elmt->id;
face = Ubc->face;
E = EL->flist[eid];
nfv = E->Nfverts(face); //OVERLOAD CALL: Nfverts: Nums.c(Tri), Nums.c(Quad), Nums.c(Tet), Nums.c(Pyr), Nums.c(Prism), Nums.c(Hex), Nums.c(Element)
for(i = 0; i < nfv; ++i){
v = E->vert + E->vnum(face,i); //OVERLOAD CALL: vnum: Nums.c(Quad), Nums.c(Tet), Nums.c(Pyr), Nums.c(Prism), Nums.c(Hex), Nums.c(Element), Nums.c(Tri)
#if 0
if(!v->solve)
if(!option("E"))
bc[v->gid-nsolve] = scal*Ubc->bvert[i];
else
bc[v->gid-nsolve] += scal*Ubc->bvert[i];
#else
bc[v->gid-nsolve] = scal*Ubc->bvert[i];
#endif
}
if(eDIM == 2){
e = E->edge + face;
skip = Ubsys->edge[e->gid] - nsolve;
if(e->l)
dsmul(e->l, scal, Ubc->bedge[0], 1, bc + skip,1);
}
else{
for(i = 0; i < nfv; ++i){
e = E->edge+E->ednum(face,i); //OVERLOAD CALL: ednum: Nums.c(Tri), Nums.c(Quad), Nums.c(Tet), Nums.c(Pyr), Nums.c(Prism), Nums.c(Hex), Nums.c(Element)
skip = Ubsys->edge[e->gid] - nsolve;
dsmul(e->l,scal,Ubc->bedge[i],1,bc+skip,1);
if(e->con) /* keep in universal format */
for(j =1;j < e->l; j+=2) bc[skip+j] *=-1.;
}
f = E->face + face;
skip = Ubsys->face[f->gid] - nsolve;
ll = (nfv == 3) ? f->l*(f->l+1)/2 : f->l*f->l;
if(ll)
dsmul(ll, scal, *Ubc->bface, 1, bc + skip,1);
}
}
}
#ifdef PARALLEL
DO_PARALLEL BCreduce(bc,Ubsys);
#endif
/* copy all values from bc to U*/
for(E=EL->fhead;E;E = E->next){
for(i = 0; i < E->Nverts; ++i)
if(!E->vert[i].solve){
#if 0
if(!option("E"))
E->vert[i].hj[0] = bc[E->vert[i].gid-nsolve];
else
E->vert[i].hj[0] = bc[E->vert[i].gid-nsolve]/mult[E->vert[i].gid-nsolve];
#else
E->vert[i].hj[0] = bc[E->vert[i].gid-nsolve];
#endif
}
for(i = 0; i < E->Nedges; ++i)
if((gid=E->edge[i].gid) >= Ubsys->ne_solve){
e = E->edge + i;
skip = Ubsys->edge[gid] - nsolve;
dcopy(e->l,bc+skip,1,e->hj,1);
if(e->con && eDIM == 3)
dscal(e->l/2,-1.0,e->hj+1,2);
}
if(eDIM == 3)
for(i = 0; i < E->Nfaces; ++i)
if((gid=E->face[i].gid) >= Ubsys->nf_solve){
f = E->face + i;
skip = Ubsys->face[gid] - nsolve;
ll = (E->Nfverts(i) == 3) ? f->l*(f->l+1)/2 : f->l*f->l; //OVERLOAD CALL: Nfverts: Nums.c(Tri), Nums.c(Quad), Nums.c(Tet), Nums.c(Pyr), Nums.c(Prism), Nums.c(Hex), Nums.c(Element)
if(ll)
dcopy(ll,bc+skip,1,*f->hj,1);
}
}
#if 0
if(option("E"))
free(mult);
#endif
free(bc);
}
#endif
void A(Element_List *, Element_List *, Bsystem *, double *, double *);
static void Precon (Element_List *, Element_List *,
Bsystem *, double *, double *);
#ifndef PARALLEL
#define MAX_ITERATIONS nsolve
#ifdef Lanczos
#define MAX_lanc_iter 990
#endif
static void
Bsolve_CG(Element_List *U, Element_List *Uf, Bsystem *B, double *p)
{
const int nsolve = B->nsolve, nvs = B->nv_solve;
int iter = 0;
double tolcg, alpha, beta, eps, rtz, rtz_old, epsfac;
static double *u = (double*)0;
static double *r = (double*)0;
static double *w = (double*)0;
static double *z = (double*)0;
static int nsol = 0, nglob = 0;
// Multi_RHS *mrhs; //OVERLOAD CALL: Multi_RHS: nekstruct.h(?), nekstruct.h(?)
#ifdef Lanczos
double *beta_l,*alpha_l;/* define the vectors to be used in the Lanczos
Algorithm */
double *D_l, *E_l; /* define additional vectors for the Lanczos Alg. */
int info =0;
if(nsolve > nsol){
alpha_l = dvector(0,MAX_ITERATIONS-1); /* Lanczos Alg */ //OVERLOAD CALL: MAX_ITERATIONS: Solve.c(?), SolveR.c(?), Solve_Stokes.c(?)
beta_l = dvector(0,MAX_ITERATIONS-1); /* Lanczos Alg */ //OVERLOAD CALL: MAX_ITERATIONS: Solve.c(?), SolveR.c(?), Solve_Stokes.c(?)
dzero (MAX_ITERATIONS,alpha_l,1); /* zero the matrices */ //OVERLOAD CALL: MAX_ITERATIONS: Solve.c(?), SolveR.c(?), Solve_Stokes.c(?)
dzero (MAX_ITERATIONS,beta_l,1); /* zero the matrices */ //OVERLOAD CALL: MAX_ITERATIONS: Solve.c(?), SolveR.c(?), Solve_Stokes.c(?)
}
#endif
if(nsolve > nsol){
if(nsol){
free(u); free(r); free(z);
}
/* Temporary arrays */
u = dvector(0,nsolve-1); /* Solution */
r = dvector(0,nsolve-1); /* residual */
z = dvector(0,nsolve-1); /* precondition solution */
nsol = nsolve;
}
if(B->nglobal > nglob){
if(nglob)
free(w);
w = dvector(0,B->nglobal-1); /* A*Search direction */
nglob = B->nglobal;
}
dzero (B->nglobal, w, 1);
dzero (nsolve, u, 1);
dcopy (nsolve, p, 1, r, 1);
if(B->singular){ /* take off mean of vertex modes */
epsfac = dsum(nvs,r,1);
epsfac /= (double)nvs;
dsadd(nvs, -epsfac, r, 1, r, 1);
}
// mrhs = Get_Multi_RHS(B,option("MULTIRHS"),nsolve,U->fhead->type);
// Mrhs_rhs(mrhs,r,r);
tolcg = (U->fhead->type == 'p')? dparam("TOLCGP"):dparam("TOLCG");
epsfac = (tol)? 1.0/tol : 1.0;
eps = sqrt(ddot(nsolve,r,1,r,1))*epsfac;
if (option("verbose") > 1)
printf("\tField %c: %3d iterations, error = %#14.6g %lg %lg\n",
U->fhead->type, iter, eps/epsfac, epsfac, tolcg);
/* =========================================================== *
* ---- Conjugate Gradient Iteration ---- *
* =========================================================== */
// Inset MRHS part (I) here
while (eps > tolcg && iter++ < MAX_ITERATIONS ) //OVERLOAD CALL: MAX_ITERATIONS: Solve.c(?), SolveR.c(?), Solve_Stokes.c(?)
{
Precon(U,Uf,B,r,z);
rtz = ddot (nsolve, r, 1, z, 1);
if (iter > 1) { /* Update search direction */
beta = rtz / rtz_old;
#ifdef Lanczos
if (iter <= MAX_lanc_iter) beta_l[iter] = beta; //OVERLOAD CALL: MAX_lanc_iter: Solve.c(?), Solve_cg.c(?)
#endif
dsvtvp(nsolve, beta, p, 1, z, 1, p, 1);
}
else
dcopy(nsolve, z, 1, p, 1);
A(U,Uf,B,p,w);
alpha = rtz/ddot(nsolve, p, 1, w, 1);
#ifdef Lanczos
if (iter <= MAX_lanc_iter) alpha_l[iter] = alpha; //OVERLOAD CALL: MAX_lanc_iter: Solve.c(?), Solve_cg.c(?)
#endif
daxpy(nsolve, alpha, p, 1, u, 1); /* Update solution... */
daxpy(nsolve,-alpha, w, 1, r, 1); /* ...and residual */
rtz_old = rtz;
eps = sqrt(ddot(nsolve, r, 1, r, 1))*epsfac; /* Compute new L2-error */
if (option("verbose") > 1)
printf("\tField %c: %3d iterations, error = %#14.6g %lg %lg\n",
U->fhead->type, iter, eps/epsfac, epsfac, tolcg);
}
/* =========================================================== *
* End of Loop *
* =========================================================== */
// printf("User time of bsystem solver (seconds): %lf \n",
// (clock()-tm)/(double)CLOCKS_PER_SEC);
#ifdef Lanczos
/* Start to set up the Lanczos matrix, T, which is a tria-diagonal
matrix and call the Lapack routine "dsterf" to compute all the
eigenvalues of a tria-diagonal matrix */
if (iter > MAX_lanc_iter){ //OVERLOAD CALL: MAX_lanc_iter: Solve.c(?), Solve_cg.c(?)
iter = MAX_lanc_iter; //OVERLOAD CALL: MAX_lanc_iter: Solve.c(?), Solve_cg.c(?)
printf("Lanczos: Conj.Gradient not fully converged, stop at 990 iter.\n");
};
D_l = dvector(0, iter-1); /* allocate memory */
E_l = dvector(0, iter-1);
dzero (iter, D_l, 1);
dzero (iter, E_l, 1);
D_l[0] = 1.0/alpha_l[1];
int i;
for (i=1; i < iter; ++i)
D_l[i] = beta_l[i+1]/alpha_l[i] + 1.0/alpha_l[i+1];
for (i=0;i<iter-1;++i)
E_l[i] = -1.0*sqrt(beta_l[i+2])/(alpha_l[i+1]);
dsterf(iter,D_l,E_l,info);
if (info) error_msg("Lanczos Method, dsterf -- info is not zero!!"); //OVERLOAD CALL: error_msg: gen_utils.h(?), hotel.h(?)
/* print maximum and minimum eigenvalues */
printf("Minimum and Maximum Eigenvalues are:%.8f,%.8f\n",D_l[0],D_l[iter-1]);
printf("Condition number(MAXeig/MINeig) is:%.8f\n",(D_l[iter-1])/D_l[0]);
#endif
/* Save solution and clean up */
dcopy(nsolve,u,1,p,1);
if (iter > MAX_ITERATIONS){ //OVERLOAD CALL: MAX_ITERATIONS: Solve.c(?), SolveR.c(?), Solve_Stokes.c(?)
error_msg (Bsolve_CG failed to converge); //OVERLOAD CALL: error_msg: gen_utils.h(?), hotel.h(?)
}
else if (option("verbose") > 1)
printf("\tField %c: %3d iterations, error = %#14.6g %lg %lg\n",
U->fhead->type, iter, eps/epsfac, epsfac, tolcg);
/* Free temporary vectors */
// free(u); free(r); free(w); free(z);
return;
}
#else
void gddot(double *alpha, double *r, double *s, double *mult,
int nsolve){
int i;
double tmp = 0.;
DO_PARALLEL{
for(i = 0, *alpha = 0.0; i < nsolve; ++i)
*alpha += mult[i]*r[i]*s[i];
gdsum(alpha,1,&tmp);
}
else
*alpha = ddot(nsolve, r, 1, s, 1);
}
static void Bsolve_CG(Element_List *U, Element_List *Uf, Bsystem *B, double *p){
register int i;
const int nsolve = B->nsolve, nvs = B->nv_solve;
int iter = 0;
double tolcg = 0.0, beta = 0.0, eps = 0.0;
double rtz_old = 0.0, epsfac = 0.0, rtz = 0.0;
double epsloc = 0.0, alpha = 0.0, tmp = 0.0;
static double *u = (double*)0;
static double *r = (double*)0;
static double *w = (double*)0;
static double *z = (double*)0;
static int nsol = 0, nglob = 0;
// Multi_RHS *mrhs; //OVERLOAD CALL: Multi_RHS: nekstruct.h(?), nekstruct.h(?)
if(nsolve > nsol){
if(nsol){
free(u); free(r); free(z);
}
/* Temporary arrays */
u = dvector(0,nsolve-1); /* Solution */
r = dvector(0,nsolve-1); /* residual */
z = dvector(0,nsolve-1); /* precondition solution */
nsol = nsolve;
}
if(B->nglobal > nglob){
if(nglob)
free(w);
w = dvector(0,B->nglobal-1); /* A*Search direction */
nglob = B->nglobal;
}
dzero (B->nglobal, w, 1);
dzero (nsolve, u, 1);
dcopy (nsolve, p, 1, r, 1);
tolcg = (U->fhead->type == 'p')? dparam("TOLCGP"):dparam("TOLCG");
epsfac = (tol)? 1.0/tol : 1.0;
// mrhs = Get_Multi_RHS(B,option("MULTIRHS"),nsolve,U->fhead->type);
/* =========================================================== *
* ---- Conjugate Gradient Iteration ---- *
* =========================================================== */
DO_PARALLEL{
int MAX_ITERATIONS = max(B->pll->nsolve, 1000); //OVERLOAD CALL: MAX_ITERATIONS: Solve.c(?), SolveR.c(?), Solve_Stokes.c(?)
parallel_gather(r,B);
if(B->singular){ /* take off mean of vertex modes */
double mean[2], work[2];
mean[0] = 0.0; mean[1] = 0.0;
for(i = 0; i < nvs; ++i) {
mean[0] += B->pll->mult[i]*r[i];
mean[1] += B->pll->mult[i];
}
gdsum(mean,2,work);
mean[0] = -1.0*mean[0]/mean[1];
dsadd(nvs, mean[0], r, 1, r, 1);
}
// Mrhs_rhs(mrhs,r,B->pll->mult);
/* sort out initial eps so that exact solution is treated properly */
gddot(&epsloc, r, r, B->pll->mult, nsolve);
eps = sqrt(epsloc)*epsfac;
if (option("verbose") > 1)
ROOTONLY printf("\tInitial eps : %lg\n",eps);
while (eps > tolcg && iter++ < MAX_ITERATIONS){ //OVERLOAD CALL: MAX_ITERATIONS: Solve.c(?), SolveR.c(?), Solve_Stokes.c(?)
Precon(U,Uf,B,r,z);
gddot(&rtz, r, z, B->pll->mult, nsolve);
if (iter > 1) { /* Update search direction */
beta = rtz / rtz_old;
dsvtvp(nsolve, beta, p, 1, z, 1, p, 1);
}
else
dcopy(nsolve, z, 1, p, 1);
/* calculate local residual */
A_fast(U,Uf,B,p,w);
alpha = ddot(B->nsolve,p,1,w,1);
gdsum(&alpha, 1, &tmp);
alpha = rtz/alpha;
parallel_gather(w,B);
daxpy(nsolve, alpha, p, 1, u, 1); /* Update solution... */
daxpy(nsolve,-alpha, w, 1, r, 1); /* ...and residual */
rtz_old = rtz;
gddot(&epsloc, r, r, B->pll->mult, nsolve);
eps = sqrt(epsloc)*epsfac; /* local L2 error */
}
// Update_Mrhs(U,Uf,B,mrhs,u);
/* =========================================================== *
* End of Loop *
* =========================================================== */
/* Save solution and clean up */
dcopy(nsolve,u,1,p,1);
if (iter > MAX_ITERATIONS){ //OVERLOAD CALL: MAX_ITERATIONS: Solve.c(?), SolveR.c(?), Solve_Stokes.c(?)
fprintf(stderr, "Iter = %d, Nprocs = %d\n",
iter, pllinfo.nprocs);
error_msg (Bsolve_CG failed to converge); //OVERLOAD CALL: error_msg: gen_utils.h(?), hotel.h(?)
}
else if (option("verbose") > 1)
ROOTONLY printf("\tField %c: %3d iterations, error = %#14.6g %lg %lg\n",
U->fhead->type, iter, eps/epsfac, epsfac, tolcg);
}
else{
int MAX_ITERATIONS = max(nsolve, 1000); //OVERLOAD CALL: MAX_ITERATIONS: Solve.c(?), SolveR.c(?), Solve_Stokes.c(?)
if(B->singular){ /* take off mean of vertex modes */
epsfac = dsum(nvs, r, 1);
epsfac /= (double)nvs;
dsadd(nvs, -epsfac, r, 1, r, 1);
}
// Mrhs_rhs(mrhs,r,r);
tolcg = (U->fhead->type == 'p')? dparam("TOLCGP"):dparam("TOLCG");
epsfac = (tol)? 1.0/tol : 1.0;
eps = sqrt(ddot(nsolve,r,1,r,1))*epsfac;
if (option("verbose") > 1)
printf("\tField %c: %3d iterations, error = %#14.6g %lg %lg\n",
U->fhead->type, iter, eps/epsfac, epsfac, tolcg);
/* =========================================================== *
* ---- Conjugate Gradient Iteration ---- *
* =========================================================== */
while (eps > tolcg && iter++ < MAX_ITERATIONS ) //OVERLOAD CALL: MAX_ITERATIONS: Solve.c(?), SolveR.c(?), Solve_Stokes.c(?)
{
Precon(U,Uf,B, r,z);
rtz = ddot (nsolve, r, 1, z, 1);
if (iter > 1) { /* Update search direction */
beta = rtz / rtz_old;
dsvtvp(nsolve, beta, p, 1, z, 1, p, 1);
}
else
dcopy(nsolve, z, 1, p, 1);
A(U,Uf,B,p,w);
alpha = rtz/ddot(nsolve, p, 1, w, 1);
daxpy(nsolve, alpha, p, 1, u, 1); /* Update solution... */
daxpy(nsolve,-alpha, w, 1, r, 1); /* ...and residual */
rtz_old = rtz;
eps = sqrt(ddot(nsolve, r, 1, r, 1))*epsfac; /* Compute new L2-error */
}
/* =========================================================== *
* End of Loop *
* =========================================================== */
// Update_Mrhs(U,Uf,B,mrhs,u);
/* Save solution and clean up */
dcopy(nsolve,u,1,p,1);
if (iter > MAX_ITERATIONS){ //OVERLOAD CALL: MAX_ITERATIONS: Solve.c(?), SolveR.c(?), Solve_Stokes.c(?)
fprintf(stderr, "Iter = %d, Nprocs = %d\n",
iter, pllinfo.nprocs);
error_msg (Bsolve_CG failed to converge); //OVERLOAD CALL: error_msg: gen_utils.h(?), hotel.h(?)
}
else if (option("verbose") > 1)
ROOTONLY printf("\tField %c: %3d iterations, error = %#14.6g %lg %lg\n",
U->fhead->type, iter, eps/epsfac, epsfac, tolcg);
}
return;
}
void parallel_gather(double *w, Bsystem *B){
gs_gop(B->pll->solve,w,"+");
}
#endif
void PreconOverlap(Element_List *U, Element_List *Uf,
Bsystem *B, double *pin, double *zout);
void PreconDirectOverlap(Bsystem *B, double *pin, double *zout);
static void Precon(Element_List *U, Element_List *Uf, Bsystem *B,
double *r, double *z){
static int counter = 0;
switch(B->Precon){
case Pre_Diag:
dvmul(B->Pmat->info.diag.ndiag,B->Pmat->info.diag.idiag,1,r,1,z,1);
break;
case Pre_Block:
break;
case Pre_None:
break;
case Pre_Overlap:
// dvmul(B->Pmat->info.overlap.ndiag,B->Pmat->info.overlap.idiag,1,r,1,z,1);
// PreconOverlap(U, Uf, B, r, z);
PreconDirectOverlap(B, r, z);
break;
}
}
void setup_signchange(Element_List *U, Bsystem *B){
double *tmp = dvector(0, U->hjtot-1);
dcopy(U->hjtot, U->base_hj, 1, tmp, 1);
Element *E;
int cnt = 0;
double *sc;
for(E=U->fhead;E;E=E->next)
cnt += E->Nbmodes;
sc = B->signchange = dvector(0, cnt-1);
dfill(U->hjtot, 1., U->base_hj, 1);
for(E=U->fhead;E;E=E->next){
E->Sign_Change(); //OVERLOAD CALL: Sign_Change: Boundary.c(Tri), Boundary.c(Quad), Boundary.c(Tet), Boundary.c(Pyr), Boundary.c(Prism), Boundary.c(Hex), Boundary.c(Element)
dcopy(E->Nbmodes, E->vert[0].hj, 1, sc, 1);
sc += E->Nbmodes;
}
dcopy(U->hjtot, tmp, 1, U->base_hj, 1);
free(tmp);
}
/* multiply out connectivity sign changes */
void SignChange(Element_List *U, Bsystem *B)
{
double *sc = B->signchange;
Element *E;
for(E=U->fhead;E;sc+=E->Nbmodes, E=E->next)
dvmul(E->Nbmodes, sc, 1, E->vert[0].hj, 1, E->vert[0].hj, 1);
}
/* Gather points from U into 'u' using the map in B->bmap */
void GathrBndry(Element_List *U, double *u, Bsystem *B){
register int i,k;
int *bmap, Nbmodes;
double *s;
dzero(B->nsolve,u,1);
double *sc = B->signchange;
for(k = 0; k < B->nel; ++k){
Nbmodes = U->flist[k]->Nbmodes;
bmap = B->bmap[k];
s = U->flist[k]->vert[0].hj;
for(i = 0; i < Nbmodes; ++i)
u[bmap[i]] += sc[i]*s[i];
sc += Nbmodes;
}
}
/* Scatter the points from u to U using the map in B->bmap */
void ScatrBndry(double *u, Element_List *U, Bsystem *B){
register int k;
int nel = B->nel, Nbmodes;
double *hj;
double *sc = B->signchange;
for(k = 0; k < nel; ++k) {
Nbmodes = U->flist[k]->Nbmodes;
hj = U->flist[k]->vert->hj;
dgathr(Nbmodes, u, B->bmap[k], hj);
dvmul (Nbmodes, sc, 1, hj, 1, hj,1);
sc += Nbmodes;
}
}
void
A(Element_List *U, Element_List *Uf, Bsystem *B, double *p, double *w){
int nel = B->nel;
double **a = B->Gmat->a;
register int k;
/* put p boundary modes into U and impose continuity */
ScatrBndry (p,Uf,B);
for(k = 0; k < nel; ++k)
dspmv('L',U->flist[k]->Nbmodes,1.0,a[U->flist[k]->geom->id],
Uf->flist[k]->vert->hj,1,0.0,U->flist[k]->vert->hj,1);
/* do sign change back and put into w */
GathrBndry (U,w,B);
}
void A_fast(Element_List *U, Element_List *Uf,
Bsystem *B, double *p, double *w){
int nel = B->nel, Nbmodes, geomid, *bmap;
double **a = B->Gmat->a, *Uhj, *Ufhj;
double *sc = B->signchange;
register int i,k;
dzero(B->nsolve, w, 1);
/* put p boundary modes into U and impose continuity */
for(k = 0; k < nel; ++k){
Uhj = U->flist[k]->vert[0].hj;
Ufhj = Uf->flist[k]->vert[0].hj;
Nbmodes = U->flist[k]->Nbmodes;
geomid = U->flist[k]->geom->id;
bmap = B->bmap[k];
// gather modes
dgathr(Nbmodes, p, bmap, Ufhj);
// fix connectivity
dvmul (Nbmodes, sc, 1, Ufhj, 1, Ufhj, 1);
// A.p
dspmv('L',Nbmodes, 1.0, a[geomid], Ufhj,1, 0.0, Uhj,1);
// scatter modes
for(i = 0; i < Nbmodes; ++i)
w[bmap[i]] += sc[i]*Uhj[i];
sc += Nbmodes;
}
}
C++ to HTML Conversion by ctoohtml