libs/dtgf/pcg_code_C99.h

// $Header: $

#include "ctypes_def.h"

#if TYPE == FLOAT 
   #define PCG spcg
#elif TYPE == DOUBLE 
   #define PCG dpcg
#elif TYPE == COMPLEXFLOAT
   #define PCG cpcg
#elif TYPE == COMPLEXDOUBLE
   #define PCG zpcg
#endif


int PCG( int n, int maxit, RTYPE tol,
         void (*amult)(int n, CTYPE *x, CTYPE *y, void **data),
         void (*msolve)(int n, CTYPE *x, CTYPE *y, void **data),
         CTYPE *b, CTYPE *x, RTYPE *rnorm1,
         void **adata, void **mdata)
{
  int k,i;
  RTYPE rnorm, bnorm;
  CTYPE gamma, alpha, beta, beta_old;
  CTYPE *r, *z, *p, *Ap;

  k = 0;
  /* Compute norm of rhs and return if rhs is zero. */
  bnorm = NRM2(n,b,1);
  if (bnorm ==  RZERO)
  {
    memset(x,0,n*sizeof(CTYPE));
    *rnorm1 = RZERO;
    return k;
  }

  r = (CTYPE *)malloc(n*sizeof(CTYPE));
  Ap = (CTYPE *)malloc(n*sizeof(CTYPE));

  /* z = A x  */
  amult(n,x,Ap,adata);

  /* r = b - A x */
  for (i=0;i<n;i++)
    r[i] = b[i] - Ap[i];
  rnorm = NRM2(n,r,1);

  /* Test if initial guess satisfies tolerance */
  if (rnorm < tol*bnorm)
  {
    *rnorm1 = rnorm/bnorm;
    free(r);
    free(Ap);
    return k;
  }

  z = (CTYPE *)malloc(n*sizeof(CTYPE));
  p = (CTYPE *)malloc(n*sizeof(CTYPE));
  /* solve M z = r  */
  if (msolve)
    msolve(n,r,z,mdata);
  else
    memcpy(z, r, n*sizeof(CTYPE));

  /* p = z  */
  memcpy(p,z,n*sizeof(CTYPE));

  /* < r, z >  */
  beta_old = DOTC(n,z,1,r,1);

  for (k=0; k<maxit && rnorm >= tol*bnorm; k++)
  {
    /* Compute A p  */
    amult(n,p,Ap,adata);

    /* < p, A p >   */
    alpha = DOTC(n,Ap,1,p,1);

    /* x = x + a*p; r = r - a*A p */
    gamma = beta_old / alpha;
    AXPY(n,gamma,p,1,x,1);
    AXPY(n,-gamma,Ap,1,r,1);
    rnorm = NRM2(n,r,1);

     /* solve M z = r  */ 
    if (msolve)
      msolve(n,r,z,mdata);
    else
      memcpy(z, r, n*sizeof(CTYPE));

    /* beta = <r, z>   */ 
    beta = DOTC(n,z,1,r,1);
    gamma = beta / beta_old; 
    beta_old = beta; 
    /* p = z + gamma*p   */ 
    for (i=0;i<n;i++) 
      p[i] = z[i] + gamma*p[i]; 
  } 
 
  free(r); 
  free(z); 
  free(p); 
  free(Ap); 

  *rnorm1 = rnorm/bnorm; 
  return k; 
}


#ifdef NOBLAS

int PCG( int n, int maxit, RTYPE tol,
           void (*amult)(int n, CTYPE *x, CTYPE *y, void **data),
           void (*msolve)(int n, CTYPE *x, CTYPE *y, void **data),
           CTYPE *b, CTYPE *x, RTYPE *rnorm1,
           void **adata, void **mdata)
{
  int k,i;
  RTYPE rnorm, bnorm;
  CTYPE gamma, alpha, beta, beta_old;
  CTYPE *r, *z, *p, *Ap;

  k=0;

  /* Compute norm of rhs and return if rhs is zero. */
  bnorm = RZERO;
  for (i=0;i<n;i++)
    bnorm = CONJ(b[i])*b[i];
  bnorm = sqrt(bnorm);
  if (bnorm ==  RZERO)
  {
    memset(x,0,n*sizeof(CTYPE));
    *rnorm1 = RZERO;
    return k;
  }

  r = (CTYPE *)malloc(n*sizeof(CTYPE));
  Ap = (CTYPE *)malloc(n*sizeof(CTYPE));

  /* z = A x  */
  amult(n,x,Ap,adata);

  /* r = b - A x */
  rnorm = RZERO;
  for (i=0;i<n;i++)
  {
    r[i] = b[i] - Ap[i];
    rnorm += CONJ(r[i])*r[i];
  }
  rnorm = sqrt(rnorm);

  /* Test if initial guess satisfies tolerance */
  if (rnorm < tol*bnorm)
  {
    *rnorm1 = rnorm/bnorm;
    free(r);
    free(Ap);
    return k;
  }

  z = (CTYPE *)malloc(n*sizeof(CTYPE));
  p = (CTYPE *)malloc(n*sizeof(CTYPE));
  /* solve M z = r  */
  if (msolve)
    msolve(n,r,z,mdata);
  else
    memcpy(z, r, n*sizeof(CTYPE));

  /* p = z  */
  memcpy(p,z,n*sizeof(CTYPE));

   /* < r, z >  */
  beta_old = CZERO;
  for (i=0; i<n; i++)
    beta_old += CONJ(z[i]) * r[i];
  
  for (k=0; k<maxit && rnorm >= tol*bnorm; k++)
  {
    /* Compute A p  */
    amult(n,p,Ap,adata);

    /* < p, A p >   */
    alpha = CZERO;
    for (i=0;i<n;i++)
      alpha += CONJ(Ap[i]) * p[i];

    /* x = x + a*p; r = r - a*A p */
    gamma = beta_old / alpha;
    rnorm = RZERO;
    for (i=0;i<n;i++) 
    {
      x[i] += gamma*p[i];
      r[i] -= gamma*Ap[i];
      rnorm += SQR(r[i]);
    } 
    rnorm = sqrt(rnorm);
    
    /* solve M z = r  */ 
    if (msolve)
      msolve(n,r,z,mdata);
    else
      memcpy(z, r, n*sizeof(CTYPE));

    /* beta = <r, z>   */ 
    beta = CZERO;
    for (i=0;i<n;i++) 
      beta += CONJ(z[i]) * r[i]; 
    gamma = beta / beta_old; 
    beta_old = beta; 
    /* p = z + b*p   */ 
    for (i=0;i<n;i++) 
      p[i] = z[i] + gamma*p[i]; 
  } 
 
  free(r); 
  free(z); 
  free(p); 
  free(Ap); 

  *rnorm1 = rnorm/bnorm; 
  return k; 
}
#endif

#undef PCG
#include "ctypes_undef.h"
Revision:	1.1
Committed:	Sun Apr 28 07:46:58 2013 UTC (10 years, 4 months ago) by tplarson
Content type:	text/plain
Branch:	MAIN
CVS Tags:	globalhs_version_5, Ver_8-7, Ver_8-5, globalhs_version_23, globalhs_version_22, globalhs_version_21, globalhs_version_20, Ver_LATEST, globalhs_version_24, Ver_8-3, globalhs_version_8, globalhs_version_9, globalhs_version_0, globalhs_version_1, globalhs_version_2, globalhs_version_3, globalhs_version_4, Ver_9-41, globalhs_version_6, globalhs_version_7, Ver_9-5, Ver_8-8, globalhs_version_19, Ver_8-2, Ver_8-10, Ver_8-1, Ver_8-6, Ver_9-1, Ver_8-4, Ver_9-2, globalhs_version_12, globalhs_version_13, globalhs_version_10, globalhs_version_11, globalhs_version_16, globalhs_version_17, globalhs_version_14, globalhs_version_15, globalhs_version_18, Ver_9-4, Ver_9-3, Ver_8-11, Ver_8-12, Ver_9-0, HEAD
Log Message:	functions needed for detrending and gapfilling
#	Content
1	// $Header: $
2
3	#include "ctypes_def.h"
4
5	#if TYPE == FLOAT
6	#define PCG spcg
7	#elif TYPE == DOUBLE
8	#define PCG dpcg
9	#elif TYPE == COMPLEXFLOAT
10	#define PCG cpcg
11	#elif TYPE == COMPLEXDOUBLE
12	#define PCG zpcg
13	#endif
14
15
16	int PCG( int n, int maxit, RTYPE tol,
17	void (amult)(int n, CTYPE x, CTYPE y, void *data),
18	void (msolve)(int n, CTYPE x, CTYPE y, void *data),
19	CTYPE b, CTYPE x, RTYPE *rnorm1,
20	void adata, void mdata)
21	{
22	int k,i;
23	RTYPE rnorm, bnorm;
24	CTYPE gamma, alpha, beta, beta_old;
25	CTYPE r, z, p, Ap;
26
27	k = 0;
28	/* Compute norm of rhs and return if rhs is zero. */
29	bnorm = NRM2(n,b,1);
30	if (bnorm == RZERO)
31	{
32	memset(x,0,n*sizeof(CTYPE));
33	*rnorm1 = RZERO;
34	return k;
35	}
36
37	r = (CTYPE )malloc(nsizeof(CTYPE));
38	Ap = (CTYPE )malloc(nsizeof(CTYPE));
39
40	/* z = A x */
41	amult(n,x,Ap,adata);
42
43	/* r = b - A x */
44	for (i=0;i<n;i++)
45	r[i] = b[i] - Ap[i];
46	rnorm = NRM2(n,r,1);
47
48	/* Test if initial guess satisfies tolerance */
49	if (rnorm < tol*bnorm)
50	{
51	*rnorm1 = rnorm/bnorm;
52	free(r);
53	free(Ap);
54	return k;
55	}
56
57	z = (CTYPE )malloc(nsizeof(CTYPE));
58	p = (CTYPE )malloc(nsizeof(CTYPE));
59	/* solve M z = r */
60	if (msolve)
61	msolve(n,r,z,mdata);
62	else
63	memcpy(z, r, n*sizeof(CTYPE));
64
65	/* p = z */
66	memcpy(p,z,n*sizeof(CTYPE));
67
68	/* < r, z > */
69	beta_old = DOTC(n,z,1,r,1);
70
71	for (k=0; k<maxit && rnorm >= tol*bnorm; k++)
72	{
73	/* Compute A p */
74	amult(n,p,Ap,adata);
75
76	/* < p, A p > */
77	alpha = DOTC(n,Ap,1,p,1);
78
79	/* x = x + ap; r = r - aA p */
80	gamma = beta_old / alpha;
81	AXPY(n,gamma,p,1,x,1);
82	AXPY(n,-gamma,Ap,1,r,1);
83	rnorm = NRM2(n,r,1);
84
85	/* solve M z = r */
86	if (msolve)
87	msolve(n,r,z,mdata);
88	else
89	memcpy(z, r, n*sizeof(CTYPE));
90
91	/* beta = <r, z> */
92	beta = DOTC(n,z,1,r,1);
93	gamma = beta / beta_old;
94	beta_old = beta;
95	/* p = z + gammap /
96	for (i=0;i<n;i++)
97	p[i] = z[i] + gamma*p[i];
98	}
99
100	free(r);
101	free(z);
102	free(p);
103	free(Ap);
104
105	*rnorm1 = rnorm/bnorm;
106	return k;
107	}
108
109
110	#ifdef NOBLAS
111
112	int PCG( int n, int maxit, RTYPE tol,
113	void (amult)(int n, CTYPE x, CTYPE y, void *data),
114	void (msolve)(int n, CTYPE x, CTYPE y, void *data),
115	CTYPE b, CTYPE x, RTYPE *rnorm1,
116	void adata, void mdata)
117	{
118	int k,i;
119	RTYPE rnorm, bnorm;
120	CTYPE gamma, alpha, beta, beta_old;
121	CTYPE r, z, p, Ap;
122
123	k=0;
124
125	/* Compute norm of rhs and return if rhs is zero. */
126	bnorm = RZERO;
127	for (i=0;i<n;i++)
128	bnorm = CONJ(b[i])*b[i];
129	bnorm = sqrt(bnorm);
130	if (bnorm == RZERO)
131	{
132	memset(x,0,n*sizeof(CTYPE));
133	*rnorm1 = RZERO;
134	return k;
135	}
136
137	r = (CTYPE )malloc(nsizeof(CTYPE));
138	Ap = (CTYPE )malloc(nsizeof(CTYPE));
139
140	/* z = A x */
141	amult(n,x,Ap,adata);
142
143	/* r = b - A x */
144	rnorm = RZERO;
145	for (i=0;i<n;i++)
146	{
147	r[i] = b[i] - Ap[i];
148	rnorm += CONJ(r[i])*r[i];
149	}
150	rnorm = sqrt(rnorm);
151
152	/* Test if initial guess satisfies tolerance */
153	if (rnorm < tol*bnorm)
154	{
155	*rnorm1 = rnorm/bnorm;
156	free(r);
157	free(Ap);
158	return k;
159	}
160
161	z = (CTYPE )malloc(nsizeof(CTYPE));
162	p = (CTYPE )malloc(nsizeof(CTYPE));
163	/* solve M z = r */
164	if (msolve)
165	msolve(n,r,z,mdata);
166	else
167	memcpy(z, r, n*sizeof(CTYPE));
168
169	/* p = z */
170	memcpy(p,z,n*sizeof(CTYPE));
171
172	/* < r, z > */
173	beta_old = CZERO;
174	for (i=0; i<n; i++)
175	beta_old += CONJ(z[i]) * r[i];
176
177	for (k=0; k<maxit && rnorm >= tol*bnorm; k++)
178	{
179	/* Compute A p */
180	amult(n,p,Ap,adata);
181
182	/* < p, A p > */
183	alpha = CZERO;
184	for (i=0;i<n;i++)
185	alpha += CONJ(Ap[i]) * p[i];
186
187	/* x = x + ap; r = r - aA p */
188	gamma = beta_old / alpha;
189	rnorm = RZERO;
190	for (i=0;i<n;i++)
191	{
192	x[i] += gamma*p[i];
193	r[i] -= gamma*Ap[i];
194	rnorm += SQR(r[i]);
195	}
196	rnorm = sqrt(rnorm);
197
198	/* solve M z = r */
199	if (msolve)
200	msolve(n,r,z,mdata);
201	else
202	memcpy(z, r, n*sizeof(CTYPE));
203
204	/* beta = <r, z> */
205	beta = CZERO;
206	for (i=0;i<n;i++)
207	beta += CONJ(z[i]) * r[i];
208	gamma = beta / beta_old;
209	beta_old = beta;
210	/* p = z + bp /
211	for (i=0;i<n;i++)
212	p[i] = z[i] + gamma*p[i];
213	}
214
215	free(r);
216	free(z);
217	free(p);
218	free(Ap);
219
220	*rnorm1 = rnorm/bnorm;
221	return k;
222	}
223	#endif
224
225	#undef PCG
226	#include "ctypes_undef.h"