/[escript]/trunk/paso/src/GMRES2.cpp
ViewVC logotype

Annotation of /trunk/paso/src/GMRES2.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4154 - (hide annotations)
Tue Jan 22 09:30:23 2013 UTC (6 years, 8 months ago) by jfenwick
Original Path: trunk/paso/src/GMRES2.c
File MIME type: text/plain
File size: 5902 byte(s)
Round 1 of copyright fixes
1 gross 1476
2 jfenwick 3981 /*****************************************************************************
3 gross 1476 *
4 jfenwick 4154 * Copyright (c) 2003-2013 by University of Queensland
5 jfenwick 3981 * http://www.uq.edu.au
6 gross 1476 *
7 ksteube 1811 * Primary Business: Queensland, Australia
8     * Licensed under the Open Software License version 3.0
9     * http://www.opensource.org/licenses/osl-3.0.php
10 gross 1476 *
11 jfenwick 3981 * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12     * Development since 2012 by School of Earth Sciences
13     *
14     *****************************************************************************/
15 gross 1476
16 ksteube 1811
17 gross 1476 #include "Common.h"
18     #include "Solver.h"
19 gross 1639 #include "PasoUtil.h"
20 gross 1476 #ifdef _OPENMP
21     #include <omp.h>
22     #endif
23    
24 gross 1639 err_t Paso_Solver_GMRES2(
25 gross 1476 Paso_Function * F,
26     const double* f0,
27     const double* x0,
28 gross 3005 double * dx,
29 gross 1476 dim_t *iter,
30     double* tolerance,
31     Paso_Performance* pp)
32     {
33 caltinay 3489 static double RENORMALIZATION_CONST=0.001;
34 gross 2987 const dim_t l=(*iter)+1, iter_max=*iter;
35     dim_t k=0,i,j;
36 gross 1639 const dim_t n=F->n;
37 gross 2987 const double rel_tol=*tolerance;
38     double abs_tol, normf0, normv, normv2, hh, hr, nu, norm_of_residual=0.;
39 gross 1476 bool_t breakFlag=FALSE, maxIterFlag=FALSE, convergeFlag=FALSE;
40     double *h=NULL, **v=NULL, *c=NULL,*s=NULL,*g=NULL, *work=NULL;
41     err_t Status=SOLVER_NO_ERROR;
42    
43     /* Test the input parameters. */
44     if (n < 0 || iter_max<=0 || l<1 || rel_tol<0) {
45     return SOLVER_INPUT_ERROR;
46     }
47 gross 3005 Paso_zeroes(n,dx);
48 gross 1476 /*
49     * allocate memory:
50     */
51     h=TMPMEMALLOC(l*l,double);
52 gross 1804 v=TMPMEMALLOC(l,double*);
53 gross 1476 c=TMPMEMALLOC(l,double);
54     s=TMPMEMALLOC(l,double);
55     g=TMPMEMALLOC(l,double);
56     work=TMPMEMALLOC(n,double);
57    
58     if (h==NULL || v ==NULL || c== NULL || s == NULL || g== NULL || work==NULL) {
59     Status=SOLVER_MEMORY_ERROR;
60     }
61     if ( Status ==SOLVER_NO_ERROR ) {
62 gross 1639 for (i=0;i<iter_max;i++) v[i]=NULL;
63 gross 1476 /*
64     * the show begins:
65     */
66     normf0=Paso_l2(n,f0,F->mpi_info);
67     k=0;
68     convergeFlag=(ABS(normf0)<=0);
69     if (! convergeFlag) {
70     abs_tol=rel_tol*normf0;
71 gross 2987 printf("GMRES2 initial residual norm %e (rel. tol=%e)\n",normf0,rel_tol);
72 gross 1639 v[0]=TMPMEMALLOC(n,double);
73     if (v[0]==NULL) {
74     Status=SOLVER_MEMORY_ERROR;
75     } else {
76     Paso_zeroes(n,v[0]);
77 gross 3005 Paso_Update(n,1.,v[0],-1./normf0,f0); /* v = -1./normf0*f0 */
78 gross 1639 g[0]=normf0;
79     }
80 gross 1798 while( (! (breakFlag || maxIterFlag || convergeFlag)) && (Status ==SOLVER_NO_ERROR) ) {
81 gross 1476 k++;
82 gross 1639 v[k]=TMPMEMALLOC(n,double);
83     if (v[k]==NULL) {
84     Status=SOLVER_MEMORY_ERROR;
85     } else {
86     /*
87     * call directional derivative function
88     */
89 gross 2987 Paso_FunctionDerivative(v[k],v[k-1],F,f0,x0,work,pp);
90 gross 1639 normv=Paso_l2(n,v[k],F->mpi_info);
91     /*
92     * Modified Gram-Schmidt
93     */
94     for (j=0;j<k;j++){
95     hh=Paso_InnerProduct(n,v[j],v[k],F->mpi_info);
96 gross 3005 Paso_Update(n,1.,v[k],(-hh),v[j]); /* v[k]-hh*v[j] */
97     h[INDEX2(j,k-1,l)]=hh;
98     /* printf("%d : %d = %e\n",k,j,hh); */
99     }
100 gross 1639 normv2=Paso_l2(n,v[k],F->mpi_info);
101     h[INDEX2(k,k-1,l)]=normv2;
102     /*
103     * reorthogonalize
104     */
105     if (! (normv + RENORMALIZATION_CONST*normv2 > normv)) {
106 gross 3005 /* printf("GMRES2: renormalization!"); */
107 gross 1639 for (j=0;j<k;j++){
108     hr=Paso_InnerProduct(n,v[j],v[k],F->mpi_info);
109 gross 3005
110 gross 1639 h[INDEX2(j,k-1,l)]+=hr;
111     Paso_Update(n,1.,v[k],(-hr),v[j]);
112     }
113     normv2=Paso_l2(n,v[k],F->mpi_info);
114     h[INDEX2(k,k-1,l)]=normv2;
115     }
116 gross 3014 /*
117 gross 3005 {int p,q;
118     for (p=0;p<k+1;p++){
119     for (q=0;q<k+1;q++)printf("%e ",h[INDEX2(p,q,l)]);
120     printf("\n");
121    
122     }
123     }
124 gross 3014 */
125 gross 1639 /*
126     * watch out for happy breakdown
127     */
128     if(normv2 > 0.) {
129 gross 3005 Paso_Update(n,1./normv2,v[k],0.,v[k]); /* normalize v[k] */
130 gross 1639 }
131     /*
132     * Form and store the information for the new Givens rotation
133     */
134     ApplyGivensRotations(k,&h[INDEX2(0,k-1,l)],c,s);
135 gross 3005
136 gross 1639 /*
137     * Don't divide by zero if solution has been found
138     */
139     g[k]=0;
140 gross 3005 nu=sqrt(h[INDEX2(k-1,k-1,l)]*h[INDEX2(k-1,k-1,l)]+h[INDEX2(k,k-1,l)]*h[INDEX2(k,k-1,l)]);
141 gross 1639 if (nu>0) {
142 gross 3005 c[k-1]= h[INDEX2(k-1,k-1,l)]/nu;
143 gross 1798 s[k-1]=-h[INDEX2(k,k-1,l)]/nu;
144     h[INDEX2(k-1,k-1,l)]=c[k-1]*h[INDEX2(k-1,k-1,l)]-s[k-1]*h[INDEX2(k,k-1,l)];
145 gross 1639 h[INDEX2(k,k-1,l)]=0;
146     ApplyGivensRotations(2,&(g[k-1]),&(c[k-1]),&(s[k-1]));
147     }
148     norm_of_residual=fabs(g[k]);
149     maxIterFlag = (k>=iter_max);
150 gross 3005 convergeFlag = (norm_of_residual <= abs_tol);
151 gross 2987 printf("GMRES2 step %d: residual %e (abs. tol=%e)\n",k,fabs(g[k]),abs_tol);
152 gross 1639 }
153 gross 1476 }
154     }
155     /*
156     * all done and ready for the forward substitution:
157     */
158 gross 3005
159 gross 1798 for (i=k-1;i>=0;--i) {
160     for (j=i+1;j<k;j++) {
161     g[i]-=h[INDEX2(i,j,l)]*g[j];
162 gross 1476 }
163     g[i]/=h[INDEX2(i,i,l)];
164 gross 3005 Paso_Update(n,1.,dx,g[i],v[i]); /* dx = dx+g[i]*v[i] */
165 gross 1476 }
166     }
167     /*
168     * clean up:
169     */
170     if ( v !=NULL) {
171 gross 1639 for (i=0;i<iter_max;i++) TMPMEMFREE(v[i]);
172 gross 1476 }
173     TMPMEMFREE(h);
174     TMPMEMFREE(v);
175     TMPMEMFREE(c);
176     TMPMEMFREE(s);
177     TMPMEMFREE(g);
178     TMPMEMFREE(work);
179     *iter=k;
180     *tolerance=norm_of_residual;
181     return Status;
182     }
183    

  ViewVC Help
Powered by ViewVC 1.1.26