CPPLapack
 All Classes Files Functions Variables Friends Pages
dgrmatrix-rci.hpp
Go to the documentation of this file.
1 //=============================================================================
2 /*! solve A*x=b for real and unsymmetric matrix using DFGMRES of Intel RCI ISS without preconditioning.\n
3  The argument is dcovector b.
4  b is overwritten and become the solution x.
5  A is not overwritten.
6 */
7 inline CPPL_INT dgrmatrix::dfgmres(dcovector& b, const double eps) const
8 {CPPL_VERBOSE_REPORT;
9 #ifndef __INTEL_COMPILER
10  ERROR_REPORT;
11  std::cerr << "dgrmatrix::dfgmres is only for intel c++ compiler (icpc)." << std::endl;
12  std::cerr << "Recompile your code with icpc to use dfgmres." << std::endl;
13  (void)b;
14  (void)eps;
15  exit(1);
16 
17 
18 #else //__INTEL_COMPILER is defined.
19  //#ifdef CPPL_DEBUG
20  if(m!=n || m!=b.l){
21  ERROR_REPORT;
22  std::cerr << "These matrix and vector cannot be solved." << std::endl
23  << "Your input was (" << m << "x" << n << ") and (" << b.l << ")." << std::endl;
24  exit(1);
25  }
26  //#endif//CPPL_DEBUG
27 
28  //////////////// constant ////////////////
29  MKL_INT itmax =1000;
30  MKL_INT n_res =500;//number of iteration before restart
31  MKL_INT N =MKL_INT(n);
32  const dcovector b_orig =b;
33 
34  //////////////// dfgmresinit ////////////////
35  dcovector x =dcovector(N).zero();//initial guess
36  //dcovector x =b;//initial guess
37  MKL_INT RCI_request;
38  MKL_INT ipar[128];
39  double dpar[128];
40  std::vector<double> tmp((2*n_res+1)*N + (n_res*(n_res+9))/2 + 1);
41  DFGMRES_INIT(&N, x.array, b.array, &RCI_request, ipar, dpar, &tmp[0]);
42  if(RCI_request){
43  ERROR_REPORT;
44  std::cerr << "dfgmres_init failed. RCI_request=" << RCI_request << "." << std::endl;
45  exit(1);
46  }
47 
48  //////////////// dfgmres_check ////////////////
49  //////// ipar ////////
50  ipar[7] = 0;//if check iteration count
51  ipar[8] = 0;//if enable residual stopping test
52  ipar[9] = 1;//if enable user defined stopping test
53  ipar[10] = 0;//if enable preconditioner
54  ipar[11] = 1;//if enable check of the norm of the next generated vector automatically
55  ipar[14] = n_res;//number of iteration before restart
56  //////// dpar ////////
57  ////dpar[0] =1e-99;//relative tolerance
58  ////dpar[1] =0.;//absolute tolerance
59  //////// call ////////
60  DFGMRES_CHECK(&N, x.array, b.array, &RCI_request, ipar, dpar, &tmp[0]);
61  if(RCI_request){
62  ERROR_REPORT;
63  std::cerr << "dfgmres_init failed. RCI_request=" << RCI_request << "." << std::endl;
64  exit(1);
65  }
66 
67  //////////////// dfgmres ////////////////
68  char transa ='N';
69  double* A =const_cast<double*>(&a[0]);
70  MKL_INT* IA =const_cast<MKL_INT*>(&ia[0]);
71  MKL_INT* JA =const_cast<MKL_INT*>(&ja[0]);
72  MKL_INT itercount;
73  dcovector x_tmp(N);
74  size_t itc =0;
75  while(1){
76  //// call ////
77  DFGMRES(&N, x.array, b.array, &RCI_request, ipar, dpar, &tmp[0]);
78 
79  //// RCI_request ////
80  if(RCI_request==1){//continue
81  MKL_DCSRGEMV(&transa, &N, A, IA, JA, &tmp[ipar[21]-1], &tmp[ipar[22]-1]);//calc A*v
82  itc++;
83  //std::cerr << "A*v" << std::endl;
84  }
85  else if(RCI_request==2){//convergence check with count check
86  if(itc%100==0){
87  ipar[12] = 1;
88  DFGMRES_GET(&N, x.array, x_tmp.array, &RCI_request, ipar, dpar, &tmp[0], &itercount);
89  double norm_r =fabs(damax(b_orig-(*this)*x_tmp));
90  std::cerr << "itc=" << itc << ": norm_r = " << norm_r << " (eps=" << eps << ")" << std::endl;
91  if(norm_r<eps){
92  std::cerr << "converged (norm_r<eps)" << std::endl;
93  break;
94  }
95  if(itc>=itmax){
96  std::cerr << "failed. (itc>=itmax)" << std::endl;
97  break;
98  }
99  }
100  }
101  else if(RCI_request==0){//converged (ipar[11])
102  std::cerr << "converged (RCI_request=0)" << std::endl;
103  break;
104  }
105  else{//failed
106  WARNING_REPORT;
107  std::cerr << "dfgmres failed. RCI_request=" << RCI_request << "." << std::endl;
108  break;
109  }
110  }
111  //////// info ////////
112  CPPL_INT info =RCI_request;
113 
114  //////////////// dfgmres_get ////////////////
115  ipar[12] =0;
116  DFGMRES_GET(&N, x.array, b.array, &RCI_request, ipar, dpar, &tmp[0], &itercount);
117 
118  //////////////// MKL_Free_Buffers ////////////////
119  MKL_Free_Buffers();
120 
121  //////////////// swap x and b ////////////////
122  swap(x, b);
123 
124  return info;
125 #endif //__INTEL_COMPILER
126 }
127 
128 ///////////////////////////////////////////////////////////////////////////////
129 ///////////////////////////////////////////////////////////////////////////////
130 ///////////////////////////////////////////////////////////////////////////////
131 ///////////////////////////////////////////////////////////////////////////////
132 ///////////////////////////////////////////////////////////////////////////////
133 ///////////////////////////////////////////////////////////////////////////////
134 ///////////////////////////////////////////////////////////////////////////////
135 ///////////////////////////////////////////////////////////////////////////////
136 ///////////////////////////////////////////////////////////////////////////////
137 
138 //=============================================================================
139 /*! solve A*x=b for real and unsymmetric matrix using DFGMRES with ILUT precondition of Intel RCI ISS.\n
140  The argument is dcovector b.
141  b is overwritten and become the solution x.
142  A is not overwritten.
143 */
144 inline CPPL_INT dgrmatrix::ilut_dfgmres(dcovector& b, const int maxfil, const double eps) const
145 {CPPL_VERBOSE_REPORT;
146 #ifndef __INTEL_COMPILER
147  ERROR_REPORT;
148  std::cerr << "dgrmatrix::ilut_dfgmres is only for intel c++ compiler (icpc)." << std::endl;
149  std::cerr << "Recompile your code with icpc to use ilut_dfgmres." << std::endl;
150  (void)b;
151  (void)maxfil;
152  (void)eps;
153  exit(1);
154 
155 
156 #else //__INTEL_COMPILER is defined.
157  //#ifdef CPPL_DEBUG
158  if(m!=n || m!=b.l){
159  ERROR_REPORT;
160  std::cerr << "These matrix and vector cannot be solved." << std::endl
161  << "Your input was (" << m << "x" << n << ") and (" << b.l << ")." << std::endl;
162  exit(1);
163  }
164  //#endif//CPPL_DEBUG
165 
166  //////////////// constants ////////////////
167  MKL_INT itmax =500;
168  MKL_INT n_res =500;//number of iteration before restart
169  //MKL_INT n_res =150;//number of iteration before restart
170  MKL_INT N =MKL_INT(n);
171  MKL_INT MAXFIL =std::min(maxfil, N-1);//limitter for maxfil
172  const dcovector b_orig =b;
173 
174  //////////////// dfgmresinit ////////////////
175  dcovector x =dcovector(b.l).zero();//initial guess
176  //dcovector x =b;//initial guess
177  MKL_INT RCI_request;
178  MKL_INT ipar[128];
179  double dpar[128];
180  std::vector<double> tmp((2*n_res+1)*N + (n_res*(n_res+9))/2 + 1);
181  DFGMRES_INIT(&N, x.array, b.array, &RCI_request, ipar, dpar, &tmp[0]);
182  if(RCI_request){
183  ERROR_REPORT;
184  std::cerr << "dfgmres_init failed. RCI_request=" << RCI_request << "." << std::endl;
185  exit(1);
186  }
187 
188  //////////////// dcsrilut ////////////////
189  double* A =const_cast<double*>(&a[0]);
190  MKL_INT* IA =const_cast<MKL_INT*>(&ia[0]);
191  MKL_INT* JA =const_cast<MKL_INT*>(&ja[0]);
192  std::vector<double> bilut((2*MAXFIL+1)*N - MAXFIL*(MAXFIL+1) + 1);
193  std::vector<MKL_INT> ibilut(N+1);
194  std::vector<MKL_INT> jbilut((2*MAXFIL+1)*N - MAXFIL*(MAXFIL+1) + 1);
195  double tol =1e-3;//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
196  //double tol =1e-4;//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
197  //double tol =1e-6;//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
198  MKL_INT ierr;
199  //ipar[30] =0;//dangerrous
200  ipar[30] =1;//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
201  dpar[30] =tol;//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
202  //dpar[30] =1e-99;//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
203  DCSRILUT(&N, A, IA, JA, &bilut[0], &ibilut[0], &jbilut[0], &tol, &MAXFIL, ipar, dpar, &ierr);
204  if(ierr){
205  ERROR_REPORT;
206  std::cerr << "dcsrilut failed. ierr=" << ierr << "." << std::endl;
207  exit(1);
208  }
209 
210  //////////////// dfgmres_check ////////////////
211  //////// ipar ////////
212  ipar[7] = 0;//if check iteration count
213  ipar[8] = 0;//if enable residual stopping test
214  ipar[9] = 1;//if enable user defined stopping test
215  ipar[10] = 1;//if enable preconditioner
216  ipar[11] = 1;//if enable check of the norm of the next generated vector automatically
217  ipar[14] = n_res;//number of iteration before restart
218  //////// dpar ////////
219  ////dpar[0] =1e-1;//relative tolerance
220  ////dpar[1] =eps;//absolute tolerance
221  //////// call ////////
222  DFGMRES_CHECK(&N, x.array, b.array, &RCI_request, ipar, dpar, &tmp[0]);
223  if(RCI_request){
224  ERROR_REPORT;
225  std::cerr << "dfgmres_init failed. RCI_request=" << RCI_request << "." << std::endl;
226  exit(1);
227  }
228 
229  //////////////// dfgmres ////////////////
230  char transa ='N';
231  char uplo1 ='L';
232  char transa1 ='N';
233  char diag1 ='U';
234  char uplo2 ='U';
235  char transa2 ='N';
236  char diag2 ='N';
237  MKL_INT itercount;
238  dcovector x_tmp(N);
239  size_t itc =0;
240  while(1){
241  //// call ////
242  DFGMRES(&N, x.array, b.array, &RCI_request, ipar, dpar, &tmp[0]);
243 
244  //// RCI_request ////
245  if(RCI_request==1){//continue
246  MKL_DCSRGEMV(&transa, &N, A, IA, JA, &tmp[ipar[21]-1], &tmp[ipar[22]-1]);//calc A*v
247  itc++;
248  }
249  else if(RCI_request==3){//precondition
250  MKL_DCSRTRSV(&uplo1, &transa1, &diag1, &N, &bilut[0], &ibilut[0], &jbilut[0], &tmp[ipar[21]-1], x_tmp.array);
251  MKL_DCSRTRSV(&uplo2, &transa2, &diag2, &N, &bilut[0], &ibilut[0], &jbilut[0], x_tmp.array, &tmp[ipar[22]-1]);
252  //std::cout << "preconditioned" << std::endl;
253  }
254  else if(RCI_request==2){//convergence check with count check
255  if(itc%100==0){
256  ipar[12] = 1;
257  DFGMRES_GET(&N, x.array, x_tmp.array, &RCI_request, ipar, dpar, &tmp[0], &itercount);
258  double norm_r =fabs(damax(b_orig-(*this)*x_tmp));
259  std::cerr << "itc=" << itc << ": norm_r = " << norm_r << " (eps=" << eps << ")" << std::endl;
260  if(norm_r<eps){
261  std::cerr << "converged (norm_r<eps)" << std::endl;
262  break;
263  }
264  if(itc>=itmax){
265  std::cerr << "failed. (itc>=itmax)" << std::endl;
266  break;
267  }
268  }
269  }
270  else if(RCI_request==0){//converged (ipar[11])
271  std::cerr << "converged (RCI_request=0)" << std::endl;
272  break;
273  }
274  else{//failed
275  WARNING_REPORT;
276  std::cerr << "dfgmres failed. RCI_request=" << RCI_request << "." << std::endl;
277  break;
278  }
279  }
280  //////// info ////////
281  CPPL_INT info =RCI_request;
282 
283  //////////////// dfgmres_get ////////////////
284  ipar[12] =0;
285  DFGMRES_GET(&N, x.array, b.array, &RCI_request, ipar, dpar, &tmp[0], &itercount);
286 
287  //////////////// MKL_Free_Buffers ////////////////
288  MKL_Free_Buffers();
289 
290  //////////////// swap x and b ////////////////
291  swap(x, b);
292 
293  return info;
294 #endif //__INTEL_COMPILER
295 }
dgrmatrix::a
std::vector< double > a
matrix component values
Definition: dgrmatrix.hpp:11
dgrmatrix::swap
friend void swap(dgrmatrix &, dgrmatrix &)
Definition: dgrmatrix-misc.hpp:105
dcovector::l
CPPL_INT l
vector size
Definition: dcovector.hpp:9
dgrmatrix::ia
std::vector< CPPL_INT > ia
rowIndex (NOT zero-based BUT one-based indexing)
Definition: dgrmatrix.hpp:12
dcovector::zero
dcovector & zero()
Definition: dcovector-misc.hpp:13
dgrmatrix::n
CPPL_INT n
matrix column size
Definition: dgrmatrix.hpp:10
dcovector::array
double * array
1D array to store vector data
Definition: dcovector.hpp:11
dcovector
Real Double-precision Column Vector Class.
Definition: dcovector.hpp:3
dgrmatrix::m
CPPL_INT m
matrix row size
Definition: dgrmatrix.hpp:9
dgrmatrix::dfgmres
CPPL_INT dfgmres(dcovector &, const double) const
Definition: dgrmatrix-rci.hpp:7
dgrmatrix::ja
std::vector< CPPL_INT > ja
columns (NOT zero-based BUT one-based indexing)
Definition: dgrmatrix.hpp:13
dgrmatrix::ilut_dfgmres
CPPL_INT ilut_dfgmres(dcovector &, const int, const double) const
Definition: dgrmatrix-rci.hpp:144
dgrmatrix::damax
friend double damax(const dgrmatrix &)
Definition: dgrmatrix-calc.hpp:3
Generated by   doxygen 1.8.6