www.pudn.com > AdRBF.rar > PBCG.h

#ifndef PBCGSOLVER 
#define PBCGSOLVER 
 
#include  
#include  
#include  
#define NR_END 1 
#define FREE_ARG char* 
#include  
#define NRANSI 
#define EPS 1.0e-14 
 
class PBCG   
{ 
public: 
  unsigned long *ija; 
  double *sa; 
  float diagonal; 
 
public: 
  void dsprstx(double sa[], unsigned long ija[], double x[], double b[], unsigned long n); 
  void dsprsax(double sa[], unsigned long ija[], double x[], double b[], unsigned long n); 
  double snrm(unsigned long n, double sx[], int itol); 
  void atimes(unsigned long n, double x[], double r[], int itrnsp); 
  void linbcg(unsigned long n, double b[], double x[], int itol, double tol, int itmax, int *iter, double *err); 
  PBCG(); 
  virtual ~PBCG(); 
  void asolve(unsigned long n, double b[], double x[], int itrnsp); 
   
   
  static void nrerror(char error_text[]) 
    /* Numerical Recipes standard error handler */ 
    { 
      fprintf(stderr,"Numerical Recipes run-time error...\n"); 
      fprintf(stderr,"%s\n",error_text); 
      fprintf(stderr,"...now exiting to system...\n"); 
      //exit(1); 
    } 
   
  static float *vector(long nl, long nh) 
    /* allocate a float vector with subscript range v[nl..nh] */ 
    { 
      float *v; 
       
      v=(float *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(float))); 
      if (!v) nrerror("allocation failure in vector()"); 
      return v-nl+NR_END; 
    } 
 
int *ivector(long nl, long nh) 
/* allocate an int vector with subscript range v[nl..nh] */ 
{ 
	int *v; 
 
	v=(int *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(int))); 
	if (!v) nrerror("allocation failure in ivector()"); 
	return v-nl+NR_END; 
} 
 
unsigned char *cvector(long nl, long nh) 
/* allocate an unsigned char vector with subscript range v[nl..nh] */ 
{ 
	unsigned char *v; 
 
	v=(unsigned char *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(unsigned char))); 
	if (!v) nrerror("allocation failure in cvector()"); 
	return v-nl+NR_END; 
} 
 
unsigned long *lvector(long nl, long nh) 
/* allocate an unsigned long vector with subscript range v[nl..nh] */ 
{ 
	unsigned long *v; 
 
	v=(unsigned long *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(long))); 
	if (!v) nrerror("allocation failure in lvector()"); 
	return v-nl+NR_END; 
} 
 
double *dvector(long nl, long nh) 
/* allocate a double vector with subscript range v[nl..nh] */ 
{ 
	double *v; 
 
	v=(double *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(double))); 
	if (!v) nrerror("allocation failure in dvector()"); 
	return v-nl+NR_END; 
} 
 
float **matrix(long nrl, long nrh, long ncl, long nch) 
/* allocate a float matrix with subscript range m[nrl..nrh][ncl..nch] */ 
{ 
	long i, nrow=nrh-nrl+1,ncol=nch-ncl+1; 
	float **m; 
 
	/* allocate pointers to rows */ 
	m=(float **) malloc((size_t)((nrow+NR_END)*sizeof(float*))); 
	if (!m) nrerror("allocation failure 1 in matrix()"); 
	m += NR_END; 
	m -= nrl; 
 
	/* allocate rows and set pointers to them */ 
	m[nrl]=(float *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(float))); 
	if (!m[nrl]) nrerror("allocation failure 2 in matrix()"); 
	m[nrl] += NR_END; 
	m[nrl] -= ncl; 
 
	for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol; 
 
	/* return pointer to array of pointers to rows */ 
	return m; 
} 
 
double **dmatrix(long nrl, long nrh, long ncl, long nch) 
/* allocate a double matrix with subscript range m[nrl..nrh][ncl..nch] */ 
{ 
	long i, nrow=nrh-nrl+1,ncol=nch-ncl+1; 
	double **m; 
 
	/* allocate pointers to rows */ 
	m=(double **) malloc((size_t)((nrow+NR_END)*sizeof(double*))); 
	if (!m) nrerror("allocation failure 1 in matrix()"); 
	m += NR_END; 
	m -= nrl; 
 
	/* allocate rows and set pointers to them */ 
	m[nrl]=(double *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(double))); 
	if (!m[nrl]) nrerror("allocation failure 2 in matrix()"); 
	m[nrl] += NR_END; 
	m[nrl] -= ncl; 
 
	for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol; 
 
	/* return pointer to array of pointers to rows */ 
	return m; 
} 
 
int **imatrix(long nrl, long nrh, long ncl, long nch) 
/* allocate a int matrix with subscript range m[nrl..nrh][ncl..nch] */ 
{ 
	long i, nrow=nrh-nrl+1,ncol=nch-ncl+1; 
	int **m; 
 
	/* allocate pointers to rows */ 
	m=(int **) malloc((size_t)((nrow+NR_END)*sizeof(int*))); 
	if (!m) nrerror("allocation failure 1 in matrix()"); 
	m += NR_END; 
	m -= nrl; 
 
 
	/* allocate rows and set pointers to them */ 
	m[nrl]=(int *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(int))); 
	if (!m[nrl]) nrerror("allocation failure 2 in matrix()"); 
	m[nrl] += NR_END; 
	m[nrl] -= ncl; 
 
	for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol; 
 
	/* return pointer to array of pointers to rows */ 
	return m; 
} 
 
float **submatrix(float **a, long oldrl, long oldrh, long oldcl, long oldch, 
	long newrl, long newcl) 
/* point a submatrix [newrl..][newcl..] to a[oldrl..oldrh][oldcl..oldch] */ 
{ 
	long i,j,nrow=oldrh-oldrl+1,ncol=oldcl-newcl; 
	float **m; 
 
	/* allocate array of pointers to rows */ 
	m=(float **) malloc((size_t) ((nrow+NR_END)*sizeof(float*))); 
	if (!m) nrerror("allocation failure in submatrix()"); 
	m += NR_END; 
	m -= newrl; 
 
	/* set pointers to rows */ 
	for(i=oldrl,j=newrl;i<=oldrh;i++,j++) m[j]=a[i]+ncol; 
 
	/* return pointer to array of pointers to rows */ 
	return m; 
} 
 
float **convert_matrix(float *a, long nrl, long nrh, long ncl, long nch) 
/* allocate a float matrix m[nrl..nrh][ncl..nch] that points to the matrix 
declared in the standard C manner as a[nrow][ncol], where nrow=nrh-nrl+1 
and ncol=nch-ncl+1. The routine should be called with the address 
&a[0][0] as the first argument. */ 
{ 
	long i,j,nrow=nrh-nrl+1,ncol=nch-ncl+1; 
	float **m; 
 
	/* allocate pointers to rows */ 
	m=(float **) malloc((size_t) ((nrow+NR_END)*sizeof(float*))); 
	if (!m) nrerror("allocation failure in convert_matrix()"); 
	m += NR_END; 
	m -= nrl; 
 
	/* set pointers to rows */ 
	m[nrl]=a-ncl; 
	for(i=1,j=nrl+1;i