CPPLapack
 All Classes Files Functions Variables Friends Pages
dssmatrix-pardiso.hpp
Go to the documentation of this file.
1 //=============================================================================
2 /*! solve A*x=b for real and symmetric positive definite matrix using Intel PARDISO.\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 dssmatrix::pardiso_definite(dcovector& b) const
8 {CPPL_VERBOSE_REPORT;
9 #ifndef __INTEL_COMPILER
10  ERROR_REPORT;
11  std::cerr << "dssmatrix::pardiso_definite is only for intel c++ compiler (icpc)." << std::endl;
12  std::cerr << "Recompile your code with icpc to use pardiso." << std::endl;
13  (void)b;
14  exit(1);
15 #else //__INTEL_COMPILER
16 
17 
18  //#ifdef CPPL_DEBUG
19  if(m!=b.l){
20  ERROR_REPORT;
21  std::cerr << "These matrix and vector cannot be solved." << std::endl
22  << "Your input was (" << m << "x" << n << ") and (" << b.l << ")." << std::endl;
23  exit(1);
24  }
25  //#endif//CPPL_DEBUG
26 
27  //////// convert matrix storage into compressed sparse row format ////////
28  //std::cerr << "converting" << std::endl;
29  std::vector<double> a(data.size());
30  std::vector<MKL_INT> ja(data.size());
31  std::vector<MKL_INT> ia(m+1);
32  ia[0] =0;
33  MKL_INT k=0;
34  for(CPPL_INT i=0; i<m; i++){
35  //// make map ////
36  const std::vector<CPPL_INT>::const_iterator line_i_end =line[i].end();
37  std::map<CPPL_INT,CPPL_INT> jc;
38  for(std::vector<CPPL_INT>::const_iterator lit=line[i].begin(); lit!=line_i_end; lit++){
39  if(data[*lit].j==i){//pardiso needs upper triangle part.
40  jc.insert( std::make_pair(data[*lit].i, *lit) );
41  }
42  }
43  //// assign ////
44  const std::map<CPPL_INT,CPPL_INT>::const_iterator jc_end =jc.end();
45  for(std::map<CPPL_INT,CPPL_INT>::const_iterator jcit=jc.begin(); jcit!=jc_end; jcit++){
46  a[k] =data[(*jcit).second].v;
47  ja[k] =MKL_INT((*jcit).first);
48  k++;
49  }
50  ia[i+1] =k;
51  }
52 
53  //////// pardisoinit ////////
54  //std::cerr << "initializing" << std::endl;
55  _MKL_DSS_HANDLE_t pt[64];
56  MKL_INT mtype =2;//real and symmetric positive definite
57  MKL_INT iparm[64];
58  PARDISOINIT(pt, &mtype, iparm);
59  iparm[1] =3;//parallel fill-in reducing ordering
60  iparm[23] =1;//use two-level scheduling factorization algorithm
61  iparm[26] =0;//disable matrix checker
62  iparm[34] =-1;//use zero-base array index
63 
64  //////// pardiso ////////
65  //std::cerr << "solving" << std::endl;
66  MKL_INT maxfct =1;
67  MKL_INT mnum =1;
68  MKL_INT phase =13;
69  MKL_INT MKL_INT_n =MKL_INT(n);
70  std::vector<MKL_INT> perm(n);
71  MKL_INT nrhs =1;
72  MKL_INT msglvl =0;
73  dcovector x(b.l);
74  MKL_INT error =1;
75  PARDISO(pt, &maxfct, &mnum, &mtype, &phase, &MKL_INT_n, &a[0], &ia[0], &ja[0], &perm[0], &nrhs, iparm, &msglvl, b.array, x.array, &error);
76  swap(b,x);//set b as x
77  if(error!=0){
78  WARNING_REPORT;
79  std::cerr << "Serious trouble happend. error = " << error << "." << std::endl;
80  }
81 
82  //////// release memory ////////
83  phase =-1;
84  MKL_INT error2 =1;
85  PARDISO(pt, &maxfct, &mnum, &mtype, &phase, &MKL_INT_n, &a[0], &ia[0], &ja[0], &perm[0], &nrhs, iparm, &msglvl, b.array, x.array, &error2);
86 
87  return error;
88 #endif //__INTEL_COMPILER
89 }
90 
91 ///////////////////////////////////////////////////////////////////////////////
92 ///////////////////////////////////////////////////////////////////////////////
93 ///////////////////////////////////////////////////////////////////////////////
94 
95 //=============================================================================
96 /*! solve A*x=b for real and symmetric indefinite matrix using Intel PARDISO.\n
97  The argument is dcovector b.
98  b is overwritten and become the solution x.
99  A is not overwritten.
100 */
101 inline CPPL_INT dssmatrix::pardiso_indefinite(dcovector& b) const
102 {CPPL_VERBOSE_REPORT;
103 #ifndef __INTEL_COMPILER
104  ERROR_REPORT;
105  std::cerr << "dssmatrix::pardiso_indefinite is only for intel c++ compiler (icpc)." << std::endl;
106  std::cerr << "Recompile your code with icpc to use pardiso." << std::endl;
107  (void)b;
108  exit(1);
109 #else //__INTEL_COMPILER
110 
111 
112  //#ifdef CPPL_DEBUG
113  if(m!=b.l){
114  ERROR_REPORT;
115  std::cerr << "These matrix and vector cannot be solved." << std::endl
116  << "Your input was (" << m << "x" << n << ") and (" << b.l << ")." << std::endl;
117  exit(1);
118  }
119  //#endif//CPPL_DEBUG
120 
121  //////// convert matrix storage into compressed sparse row format ////////
122  //std::cerr << "converting" << std::endl;
123  std::vector<double> a(data.size());
124  std::vector<MKL_INT> ja(data.size());
125  std::vector<MKL_INT> ia(m+1);
126  ia[0] =0;
127  MKL_INT k=0;
128  for(CPPL_INT i=0; i<m; i++){
129  //// make map ////
130  const std::vector<CPPL_INT>::const_iterator line_i_end =line[i].end();
131  std::map<CPPL_INT,CPPL_INT> jc;
132  for(std::vector<CPPL_INT>::const_iterator lit=line[i].begin(); lit!=line_i_end; lit++){
133  if(data[*lit].j==i){//pardiso needs upper triangle part.
134  jc.insert( std::make_pair(data[*lit].i, *lit) );
135  }
136  }
137  //// assign ////
138  const std::map<CPPL_INT,CPPL_INT>::const_iterator jc_end =jc.end();
139  for(std::map<CPPL_INT,CPPL_INT>::const_iterator jcit=jc.begin(); jcit!=jc_end; jcit++){
140  a[k] =data[(*jcit).second].v;
141  ja[k] =MKL_INT((*jcit).first);
142  k++;
143  }
144  ia[i+1] =k;
145  }
146 
147  //////// pardisoinit ////////
148  //std::cerr << "initializing" << std::endl;
149  _MKL_DSS_HANDLE_t pt[64];
150  MKL_INT mtype =-2;//real and symmetric indefinite
151  MKL_INT iparm[64];
152  PARDISOINIT(pt, &mtype, iparm);
153  iparm[1] =3;//parallel fill-in reducing ordering
154  iparm[23] =1;//use two-level scheduling factorization algorithm
155  iparm[26] =0;//disable matrix checker
156  iparm[34] =-1;//use zero-base array index
157 
158  //////// pardiso ////////
159  //std::cerr << "solving" << std::endl;
160  MKL_INT maxfct =1;
161  MKL_INT mnum =1;
162  MKL_INT phase =13;
163  MKL_INT MKL_INT_n =MKL_INT(n);
164  std::vector<MKL_INT> perm(n);
165  MKL_INT nrhs =1;
166  MKL_INT msglvl =0;
167  dcovector x(b.l);
168  MKL_INT error =1;
169  PARDISO(pt, &maxfct, &mnum, &mtype, &phase, &MKL_INT_n, &a[0], &ia[0], &ja[0], &perm[0], &nrhs, iparm, &msglvl, b.array, x.array, &error);
170  swap(b,x);//set b as x
171  if(error!=0){
172  WARNING_REPORT;
173  std::cerr << "Serious trouble happend. error = " << error << "." << std::endl;
174  }
175 
176  //////// release memory ////////
177  phase =-1;
178  MKL_INT error2 =1;
179  PARDISO(pt, &maxfct, &mnum, &mtype, &phase, &MKL_INT_n, &a[0], &ia[0], &ja[0], &perm[0], &nrhs, iparm, &msglvl, b.array, x.array, &error2);
180 
181  return error;
182 #endif //__INTEL_COMPILER
183 }
dssmatrix::data
std::vector< dcomponent > data
matrix data
Definition: dssmatrix.hpp:11
dssmatrix::line
std::vector< std::vector< CPPL_INT > > line
vector of vector to store the entry information of component for each row and column ...
Definition: dssmatrix.hpp:12
dssmatrix::pardiso_indefinite
CPPL_INT pardiso_indefinite(dcovector &) const
Definition: dssmatrix-pardiso.hpp:101
dssmatrix::swap
friend void swap(dssmatrix &, dssmatrix &)
Definition: dssmatrix-misc.hpp:390
dcovector::l
CPPL_INT l
vector size
Definition: dcovector.hpp:9
i
_dgematrix i(const _dgbmatrix &mat)
Definition: _dgbmatrix-calc.hpp:20
dcovector::array
double * array
1D array to store vector data
Definition: dcovector.hpp:11
dssmatrix::pardiso_definite
CPPL_INT pardiso_definite(dcovector &) const
Definition: dssmatrix-pardiso.hpp:7
dcovector
Real Double-precision Column Vector Class.
Definition: dcovector.hpp:3
dssmatrix::n
CPPL_INT n
matrix column size
Definition: dssmatrix.hpp:10
dssmatrix::m
CPPL_INT const & m
matrix row size
Definition: dssmatrix.hpp:9
Generated by   doxygen 1.8.6