00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019 #ifndef _GMRES_H
00020 #define _GMRES_H
00021 #include "vec.h"
00022 #include "smat.h"
00023 #include "mat.h"
00024 #include "sparsemat.h"
00025 #include "nonclass.h"
00026 #include <iomanip>
00027 #include "auxstorage.h"
00028
00029 namespace LA {
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041 template<class T>
00042 void gmres_backsubstitute(const NRMat<T> &R, NRVec<T> &c, const NRVec<T> &d, const int k)
00043 {
00044 c.copyonwrite();
00045 if(R(k,k)==0.) laerror("singular matrix in gmres triangular solution");
00046 c[k] = d[k]/R(k,k);
00047 for (int i=k-1;i>=0;i--) c[i] = (d[i]-xdot(k-i,&R(i,i+1),1,&c[i+1],1)) / R(i,i);
00048 }
00049
00050
00051
00052 template<typename T, typename Matrix>
00053 bool gmres(const Matrix &bigmat, const NRVec<T> &b, NRVec<T> &x, const bool doguess=1, const double eps=1e-7, const int MAXIT=50, const bool verbose=1, bool square=1,const bool precondition=1, int neustart=0, const int incore=1)
00054 {
00055 int zeilen=bigmat.nrows();
00056 int spalten=bigmat.ncols();
00057 if(spalten==1) laerror("gmres does not work for n==1, use conjgrad if you need this trivial case");
00058 if(x.size()!=spalten || b.size() != zeilen) laerror("incompatible vectors and matrix sizes in GMRES");
00059
00060 if(zeilen!=spalten) square=0;
00061 if(!neustart) neustart = zeilen/10;
00062 if (neustart < 10) neustart = 10;
00063 x.copyonwrite();
00064
00065 bool flag;
00066 double beta,beta_0;
00067 double d_alt=0;
00068
00069 AuxStorage<T> *st;
00070 NRVec<T> *v;
00071 NRVec<T> r_k(spalten),z(spalten);
00072 NRVec<T> cci(MAXIT+1),ssi(MAXIT+1),c(MAXIT+1),d(MAXIT+1);
00073 NRMat<T> H(MAXIT+1,MAXIT+1);
00074 T ci,si;
00075 v = new NRVec<T>[incore?MAXIT+1:1];
00076 st = incore?NULL:new AuxStorage<T>;
00077
00078 if(doguess)
00079 {
00080 bigmat.gemv(0,x,'t',-1.,b);
00081 if(precondition) bigmat.diagonalof(x,true);
00082 x.normalize();
00083 }
00084
00085 neustart:
00086 for (int l=0;l<neustart;l++)
00087 {
00088 if(square)
00089 {
00090 bigmat.gemv(0,r_k,'n',1,x);
00091 r_k -= b;
00092 }
00093 else
00094 {
00095 NRVec<T> dum(zeilen);
00096 bigmat.gemv(0,dum,'n',1,x);
00097 bigmat.gemv(0,r_k,'t',1,dum);
00098 bigmat.gemv(0,z,'t',-1.,b);
00099 r_k += z;
00100 }
00101
00102 if(precondition) bigmat.diagonalof(r_k,true);
00103
00104 beta = r_k.norm();
00105 if(l==0) beta_0 = beta;
00106 v[0] = r_k* (1./beta);
00107 if(!incore) st->put(v[0],0);
00108
00109
00110 for (int k=0;k<MAXIT;k++)
00111 {
00112
00113
00114
00115
00116 if(!incore) st->get(v[0],k);
00117 if(square)
00118 {
00119 bigmat.gemv(0,z,'n',1,v[incore?k:0]);
00120 }
00121 else
00122 {
00123 NRVec<T> dum(zeilen);
00124 bigmat.gemv(0,dum,'n',1,v[incore?k:0]);
00125 bigmat.gemv(0,z,'t',1,dum);
00126 }
00127 if(precondition) bigmat.diagonalof(z,true);
00128
00129
00130 for (int i=0;i<=k;i++)
00131 {
00132 if(!incore) st->get(v[0],i);
00133 H(i,k) = z*v[incore?i:0];
00134 z.axpy(-H(i,k),v[incore?i:0]);
00135 }
00136
00137
00138 double tmp;
00139 H(k+1,k) = tmp= z.norm();
00140 if(tmp < 1.e-2*eps )
00141 {
00142 if(verbose) std::cerr <<("gmres restart performed\n");
00143
00144 for (int i=0;i<k;i++)
00145 {
00146 ci = cci[i];si = ssi[i];
00147 for (int j=0;j<k;j++)
00148 {
00149 T a = H(i,j);
00150 H(i,j) = ci*a+si*H(i+1,j);
00151 H(i+1,j) = -si*a+ci*H(i+1,j);
00152 }
00153 }
00154
00155 d *= -1.;
00156 gmres_backsubstitute(H,c,d,k-1);
00157 for (int i=0;i<k-1;i++)
00158 {
00159 if(!incore) st->get(v[0],i);
00160 x.axpy(c[i],v[incore?i:0]);
00161 }
00162 flag=0; goto neustart;
00163 }
00164
00165 v[incore?k+1:0] = z * (1./H(k+1,k));
00166 if(!incore) st->put(v[0],k+1);
00167
00168
00169 for (int j=0;j<k+2;j++) d[j] = H(j,k);
00170 for (int i=0;i<k;i++)
00171 {
00172 ci = cci[i];
00173 si = ssi[i];
00174 T a = d[i];
00175 d[i] = ci*a+si*d[i+1];
00176 d[i+1] = -si*a+ci*d[i+1];
00177 }
00178
00179 ci=hypot(d[k],d[k+1]);
00180 cci[k]=d[k]/ci;
00181 ssi[k]=d[k+1]/ci;
00182
00183
00184 d= 0.;
00185 d[0]=beta;
00186 for (int i=0;i<=k;i++)
00187 {
00188 ci = cci[i];si = ssi[i];
00189 T a = d[i];
00190 d[i] = ci*a+si*d[i+1];
00191 d[i+1] = -si*a+ci*d[i+1];
00192 }
00193
00194
00195 if(verbose)
00196 {
00197 std::cout << "gmres iter "<<l<<" "<<k<<" resid "
00198 <<std::setw(0)<<std::setiosflags(std::ios::scientific)<<std::setprecision(8)
00199 <<std::abs(d[k+1])<< " thr "<<eps*beta_0<< " reduction "
00200 <<std::setw(5)<<std::setprecision(2)<<std::resetiosflags(std::ios::scientific)
00201 <<(d_alt - std::abs(d[k+1]))/d_alt*100<< "\n" <<std::setprecision(12);
00202 std::cout.flush();
00203 }
00204
00205 d_alt = abs(d[k+1]);
00206
00207 if (d_alt < eps*beta_0)
00208 {
00209
00210 for (int i=0;i<k;i++)
00211 {
00212 ci = cci[i];
00213 si = ssi[i];
00214 for (int j=0;j<k;j++)
00215 {
00216 T a = H(i,j);
00217 H(i,j) = ci*a+si*H(i+1,j);
00218 H(i+1,j) = -si*a+ci*H(i+1,j);
00219 }
00220 }
00221
00222
00223 d *= -1.;
00224 gmres_backsubstitute(H,c,d,k-1);
00225 for(int i=0;i<k;i++)
00226 {
00227 if(!incore) st->get(v[0],i);
00228 x.axpy(c[i],v[incore?i:0]);
00229 }
00230 flag=0; goto myreturn;
00231 }
00232 }
00233
00234
00235 for (int i=0;i<MAXIT;i++)
00236 {
00237 ci = cci[i];si = ssi[i];
00238 for (int j=0;j<MAXIT;j++)
00239 {
00240 T a = H(i,j);
00241 H(i,j) = ci*a+si*H(i+1,j);
00242 H(i+1,j) = -si*a+ci*H(i+1,j);
00243 }
00244 }
00245
00246
00247 d *= -1.;
00248 gmres_backsubstitute(H,c,d,MAXIT-1);
00249 for(int i=0;i<MAXIT;i++)
00250 {
00251 if(!incore) st->get(v[0],i);
00252 x.axpy(c[i],v[incore?i:0]);
00253 }
00254
00255 }
00256 flag=1;
00257
00258 myreturn:
00259 delete[] v;
00260 if(!incore) delete st;
00261
00262 return !flag;
00263 }
00264
00265 }
00266
00267 #endif