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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4657 - (hide annotations)
Thu Feb 6 06:12:20 2014 UTC (5 years, 9 months ago) by jfenwick
File size: 5883 byte(s)
I changed some files.
Updated copyright notices, added GeoComp.



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

Properties

Name Value
svn:mergeinfo /branches/amg_from_3530/paso/src/GMRES2.cpp:3531-3826 /branches/lapack2681/paso/src/GMRES2.cpp:2682-2741 /branches/pasowrap/paso/src/GMRES2.cpp:3661-3674 /branches/py3_attempt2/paso/src/GMRES2.cpp:3871-3891 /branches/restext/paso/src/GMRES2.cpp:2610-2624 /branches/ripleygmg_from_3668/paso/src/GMRES2.cpp:3669-3791 /branches/stage3.0/paso/src/GMRES2.cpp:2569-2590 /branches/symbolic_from_3470/paso/src/GMRES2.cpp:3471-3974 /branches/symbolic_from_3470/ripley/test/python/paso/src/GMRES2.cpp:3517-3974 /release/3.0/paso/src/GMRES2.cpp:2591-2601 /trunk/paso/src/GMRES2.cpp:4257-4344 /trunk/ripley/test/python/paso/src/GMRES2.cpp:3480-3515

  ViewVC Help
Powered by ViewVC 1.1.26