#define NRANSI#include "nrutil.h"#define NPRE 1#define NPOST 1#define ALPHA 0.33#define NGMAX 15void mgfas(double **u, int n, int maxcyc){	double anorm2(double **a, int n);	void copy(double **aout, double **ain, int n);	void interp(double **uf, double **uc, int nf);	void lop(double **out, double **u, int n);	void matadd(double **a, double **b, double **c, int n);	void matsub(double **a, double **b, double **c, int n);	void relax2(double **u, double **rhs, int n);	void rstrct(double **uc, double **uf, int nc);	void slvsm2(double **u, double **rhs);	unsigned int j,jcycle,jj,jm1,jpost,jpre,nf,ng=0,ngrid,nn;	double **irho[NGMAX+1],**irhs[NGMAX+1],**itau[NGMAX+1],		**itemp[NGMAX+1],**iu[NGMAX+1];	double res,trerr;	nn=n;	while (nn >>= 1) ng++;	if (n != 1+(1L << ng)) nrerror("n-1 must be a power of 2 in mgfas.");	if (ng > NGMAX) nrerror("increase NGMAX in mglin.");	nn=n/2+1;	ngrid=ng-1;	irho[ngrid]=dmatrix(1,nn,1,nn);	rstrct(irho[ngrid],u,nn);	while (nn > 3) {		nn=nn/2+1;		irho[--ngrid]=dmatrix(1,nn,1,nn);		rstrct(irho[ngrid],irho[ngrid+1],nn);	}	nn=3;	iu[1]=dmatrix(1,nn,1,nn);	irhs[1]=dmatrix(1,nn,1,nn);	itau[1]=dmatrix(1,nn,1,nn);	itemp[1]=dmatrix(1,nn,1,nn);	slvsm2(iu[1],irho[1]);	free_dmatrix(irho[1],1,nn,1,nn);	ngrid=ng;	for (j=2;j<=ngrid;j++) {		nn=2*nn-1;		iu[j]=dmatrix(1,nn,1,nn);		irhs[j]=dmatrix(1,nn,1,nn);		itau[j]=dmatrix(1,nn,1,nn);		itemp[j]=dmatrix(1,nn,1,nn);		interp(iu[j],iu[j-1],nn);		copy(irhs[j],(j != ngrid ? irho[j] : u),nn);		for (jcycle=1;jcycle<=maxcyc;jcycle++) {		nf=nn;			for (jj=j;jj>=2;jj--) {				for (jpre=1;jpre<=NPRE;jpre++)					relax2(iu[jj],irhs[jj],nf);				lop(itemp[jj],iu[jj],nf);				nf=nf/2+1;				jm1=jj-1;				rstrct(itemp[jm1],itemp[jj],nf);				rstrct(iu[jm1],iu[jj],nf);				lop(itau[jm1],iu[jm1],nf);				matsub(itau[jm1],itemp[jm1],itau[jm1],nf);				if (jj == j)					trerr=ALPHA*anorm2(itau[jm1],nf);				rstrct(irhs[jm1],irhs[jj],nf);				matadd(irhs[jm1],itau[jm1],irhs[jm1],nf);			}			slvsm2(iu[1],irhs[1]);			nf=3;			for (jj=2;jj<=j;jj++) {			jm1=jj-1;			rstrct(itemp[jm1],iu[jj],nf);			matsub(iu[jm1],itemp[jm1],itemp[jm1],nf);			nf=2*nf-1;			interp(itau[jj],itemp[jm1],nf);			matadd(iu[jj],itau[jj],iu[jj],nf);			for (jpost=1;jpost<=NPOST;jpost++)				relax2(iu[jj],irhs[jj],nf);			}			lop(itemp[j],iu[j],nf);			matsub(itemp[j],irhs[j],itemp[j],nf);			res=anorm2(itemp[j],nf);			if (res < trerr) break;		}	}	copy(u,iu[ngrid],n);	for (nn=n,j=ng;j>=1;j--,nn=nn/2+1) {		free_dmatrix(itemp[j],1,nn,1,nn);		free_dmatrix(itau[j],1,nn,1,nn);		free_dmatrix(irhs[j],1,nn,1,nn);		free_dmatrix(iu[j],1,nn,1,nn);		if (j != ng && j != 1) free_dmatrix(irho[j],1,nn,1,nn);	}}#undef NGMAX#undef NPRE#undef NPOST#undef ALPHA#undef NRANSI