CPPLapack
 All Classes Files Functions Variables Friends Pages
Public Member Functions | Public Attributes | Friends | List of all members
dgrmatrix Class Reference

Real Double-precision General Compressed Sparse Row (CSR) Matrix Class. More...

#include <dgrmatrix.hpp>

Public Member Functions

 dgrmatrix ()
 
 dgrmatrix (const CPPL_INT &, const CPPL_INT &)
 
 dgrmatrix (const dgrmatrix &)
 
 dgrmatrix (const char *)
 
 ~dgrmatrix ()
 
_dgematrix to_dgematrix () const
 
double operator() (const CPPL_INT &, const CPPL_INT &) const
 
double & operator() (const CPPL_INT &, const CPPL_INT &)
 
void write (const char *) const
 
void read (const char *)
 
void clear ()
 
dgrmatrixzero ()
 
void copy (const dgrmatrix &)
 
bool isListed (const CPPL_INT &, const CPPL_INT &) const
 
void checkup ()
 
CPPL_INT pardiso (dcovector &) const
 
CPPL_INT dfgmres (dcovector &, const double) const
 
CPPL_INT ilut_dfgmres (dcovector &, const int, const double) const
 
dgrmatrixoperator= (const dgrmatrix &)
 
dgrmatrixoperator*= (const double &)
 
dgrmatrixoperator/= (const double &)
 

Public Attributes

CPPL_INT m
 matrix row size More...
 
CPPL_INT n
 matrix column size More...
 
std::vector< double > a
 matrix component values More...
 
std::vector< CPPL_INT > ia
 rowIndex (NOT zero-based BUT one-based indexing) More...
 
std::vector< CPPL_INT > ja
 columns (NOT zero-based BUT one-based indexing) More...
 

Friends

std::ostream & operator<< (std::ostream &, const dgrmatrix &)
 
void swap (dgrmatrix &, dgrmatrix &)
 
double damax (const dgrmatrix &)
 
dgrmatrix operator* (const dgrmatrix &, const double &)
 
dgrmatrix operator* (const double &, const dgrmatrix &)
 
_dcovector operator* (const dgrmatrix &, const dcovector &)
 
_dcovector operator* (const dgrmatrix &, const _dcovector &)
 
dgrmatrix operator/ (const dgrmatrix &, const double &)
 

Detailed Description

Real Double-precision General Compressed Sparse Row (CSR) Matrix Class.

Definition at line 3 of file dgrmatrix.hpp.

Constructor & Destructor Documentation

dgrmatrix::dgrmatrix ( )
inline

dgrmatrix constructor without arguments

Definition at line 3 of file dgrmatrix-constructor.hpp.

References m, and n.

4 {CPPL_VERBOSE_REPORT;
5  m =0;
6  n =0;
7 }
CPPL_INT n
matrix column size
Definition: dgrmatrix.hpp:10
CPPL_INT m
matrix row size
Definition: dgrmatrix.hpp:9
dgrmatrix::dgrmatrix ( const CPPL_INT &  ,
const CPPL_INT &   
)
inline
dgrmatrix::dgrmatrix ( const dgrmatrix mat)
inline

dgrmatrix copy constructor

Definition at line 15 of file dgrmatrix-constructor.hpp.

References copy().

16 {CPPL_VERBOSE_REPORT;
17  copy(mat);
18 }
void copy(const dgrmatrix &)
dgrmatrix::dgrmatrix ( const char *  filename)
inline

dgrmatrix constructor with filename

Definition at line 22 of file dgrmatrix-constructor.hpp.

References read().

23 {CPPL_VERBOSE_REPORT;
24  read(filename);
25 }
void read(const char *)
dgrmatrix::~dgrmatrix ( )
inline

dgrmatrix destructor

Definition at line 33 of file dgrmatrix-constructor.hpp.

References a, ia, and ja.

34 {CPPL_VERBOSE_REPORT;
35  a.clear();
36  ia.clear();
37  ja.clear();
38 }
std::vector< double > a
matrix component values
Definition: dgrmatrix.hpp:11
std::vector< CPPL_INT > ia
rowIndex (NOT zero-based BUT one-based indexing)
Definition: dgrmatrix.hpp:12
std::vector< CPPL_INT > ja
columns (NOT zero-based BUT one-based indexing)
Definition: dgrmatrix.hpp:13

Member Function Documentation

_dgematrix dgrmatrix::to_dgematrix ( ) const
inline

cast to _zgrmatrix

convert to _dgematrix

Definition at line 28 of file dgrmatrix-cast.hpp.

References _(), a, i(), ia, ja, m, n, and dgematrix::zero().

29 {CPPL_VERBOSE_REPORT;
30  dgematrix newmat(m,n);
31  newmat.zero();
32 
33  for(CPPL_INT i=0; i<m; i++){
34  int k_beg =ia[i]-1;
35  int k_end =ia[i+1]-1;
36  for(int k=k_beg; k<k_end; k++){
37  int j =ja[k]-1;
38  newmat(i,j) =a[k];
39  }
40  }
41 
42  return _(newmat);
43 }
std::vector< double > a
matrix component values
Definition: dgrmatrix.hpp:11
std::vector< CPPL_INT > ia
rowIndex (NOT zero-based BUT one-based indexing)
Definition: dgrmatrix.hpp:12
_dgematrix i(const _dgbmatrix &mat)
Real Double-precision General Dence Matrix Class.
Definition: dgematrix.hpp:3
CPPL_INT n
matrix column size
Definition: dgrmatrix.hpp:10
CPPL_INT m
matrix row size
Definition: dgrmatrix.hpp:9
std::vector< CPPL_INT > ja
columns (NOT zero-based BUT one-based indexing)
Definition: dgrmatrix.hpp:13
_dcovector _(dcovector &vec)
double dgrmatrix::operator() ( const CPPL_INT &  i,
const CPPL_INT &  j 
) const
inline

operator() for const object

Definition at line 3 of file dgrmatrix-io.hpp.

References a, i(), ia, ja, m, and n.

4 {CPPL_VERBOSE_REPORT;
5 #ifdef CPPL_DEBUG
6  if( i<0 || j<0 || m<=i || n<=j ){
7  ERROR_REPORT;
8  std::cerr << "The required component is out of the matrix size." << std::endl
9  << "Your input is (" << i << "," << j << "), whereas the matrix size is " << m << "x" << n << "." << std::endl;
10  exit(1);
11  }
12 #endif//CPPL_DEBUG
13 
14  //// search (i,j) component ////
15  int k_beg =ia[i]-1;
16  int k_end =ia[i+1]-1;
17  for(int k=k_beg; k<k_end; k++){
18  if(j==ja[k]-1){
19  return a[k];
20  }
21  }
22 
23  //// (i,j) component was not found ////
24  return 0.0;
25 }
std::vector< double > a
matrix component values
Definition: dgrmatrix.hpp:11
std::vector< CPPL_INT > ia
rowIndex (NOT zero-based BUT one-based indexing)
Definition: dgrmatrix.hpp:12
_dgematrix i(const _dgbmatrix &mat)
CPPL_INT n
matrix column size
Definition: dgrmatrix.hpp:10
CPPL_INT m
matrix row size
Definition: dgrmatrix.hpp:9
std::vector< CPPL_INT > ja
columns (NOT zero-based BUT one-based indexing)
Definition: dgrmatrix.hpp:13
double & dgrmatrix::operator() ( const CPPL_INT &  i,
const CPPL_INT &  j 
)
inline

operator() for const object

Definition at line 29 of file dgrmatrix-io.hpp.

References a, i(), ia, ja, m, and n.

30 {CPPL_VERBOSE_REPORT;
31 #ifdef CPPL_DEBUG
32  if( i<0 || j<0 || m<=i || n<=j ){
33  ERROR_REPORT;
34  std::cerr << "The required component is out of the matrix size." << std::endl
35  << "Your input is (" << i << "," << j << "), whereas the matrix size is " << m << "x" << n << "." << std::endl;
36  exit(1);
37  }
38 #endif//CPPL_DEBUG
39 
40  //// search (i,j) component ////
41  int k_beg =ia[i]-1;
42  int k_end =ia[i+1]-1;
43  for(int k=k_beg; k<k_end; k++){
44  if(j==ja[k]-1){
45  return a[k];
46  }
47  }
48 
49  //// (i,j) component was not found ////
50  ERROR_REPORT;
51  std::cerr << "dgrmatrix does not allow component addition with operator()." << std::endl;
52  exit(1);
53 }
std::vector< double > a
matrix component values
Definition: dgrmatrix.hpp:11
std::vector< CPPL_INT > ia
rowIndex (NOT zero-based BUT one-based indexing)
Definition: dgrmatrix.hpp:12
_dgematrix i(const _dgbmatrix &mat)
CPPL_INT n
matrix column size
Definition: dgrmatrix.hpp:10
CPPL_INT m
matrix row size
Definition: dgrmatrix.hpp:9
std::vector< CPPL_INT > ja
columns (NOT zero-based BUT one-based indexing)
Definition: dgrmatrix.hpp:13
void dgrmatrix::write ( const char *  filename) const
inline

Definition at line 89 of file dgrmatrix-io.hpp.

References a, ia, ja, m, and n.

90 {CPPL_VERBOSE_REPORT;
91  std::ofstream ofs(filename, std::ios::trunc);
92  ofs.setf(std::cout.flags());
93  ofs.precision(std::cout.precision());
94  ofs.width(std::cout.width());
95  ofs.fill(std::cout.fill());
96 
97  ofs << "#dgrmatrix" << " " << m << " " << n << " " << a.size() << std::endl;
98 
99  size_t a_size =a.size();
100  for(size_t k=0; k<a_size; k++){
101  ofs << " " << a[k];
102  }
103  ofs << "\n";
104 
105  size_t ia_size =ia.size();
106  for(size_t k=0; k<ia_size; k++){
107  ofs << " " << ia[k];
108  }
109  ofs << "\n";
110 
111  size_t ja_size =ja.size();
112  for(size_t k=0; k<ja_size; k++){
113  ofs << " " << ja[k];
114  }
115  ofs << "\n" << std::flush;
116 
117  ofs.close();
118 }
std::vector< double > a
matrix component values
Definition: dgrmatrix.hpp:11
std::vector< CPPL_INT > ia
rowIndex (NOT zero-based BUT one-based indexing)
Definition: dgrmatrix.hpp:12
CPPL_INT n
matrix column size
Definition: dgrmatrix.hpp:10
CPPL_INT m
matrix row size
Definition: dgrmatrix.hpp:9
std::vector< CPPL_INT > ja
columns (NOT zero-based BUT one-based indexing)
Definition: dgrmatrix.hpp:13
void dgrmatrix::read ( const char *  filename)
inline

Definition at line 121 of file dgrmatrix-io.hpp.

References a, ia, ja, m, and n.

Referenced by dgrmatrix().

122 {CPPL_VERBOSE_REPORT;
123  std::ifstream s( filename );
124  if(!s){
125  ERROR_REPORT;
126  std::cerr << "The file \"" << filename << "\" can not be opened." << std::endl;
127  exit(1);
128  }
129 
130  std::string id;
131  s >> id;
132  if( id != "dgrmatrix" && id != "#dgrmatrix" ){
133  ERROR_REPORT;
134  std::cerr << "The type name of the file \"" << filename << "\" is not dgrmatrix." << std::endl
135  << "Its type name was " << id << " ." << std::endl;
136  exit(1);
137  }
138 
139  //////// read ////////
140  size_t a_size;
141  s >> m >> n >> a_size;
142  a.resize(a_size);
143  ia.resize(m+1);
144  ja.resize(a_size);
145 
146  for(size_t k=0; k<a_size; k++){
147  s >> a[k];
148  if(s.fail()){
149  ERROR_REPORT;
150  exit(1);
151  }
152  }
153  for(CPPL_INT k=0; k<=m; k++){
154  s >> ia[k];
155  if(s.fail()){
156  ERROR_REPORT;
157  exit(1);
158  }
159  }
160  for(size_t k=0; k<a_size; k++){
161  s >> ja[k];
162  if(s.fail()){
163  ERROR_REPORT;
164  exit(1);
165  }
166  }
167 
168  //////// garbage ////////
169  s >> a_size;
170  if(!s.eof()){
171  ERROR_REPORT;
172  std::cerr << "There is something is wrong with the file \"" << filename << " ." << std::endl;
173  exit(1);
174  }
175  s.close();
176 }
std::vector< double > a
matrix component values
Definition: dgrmatrix.hpp:11
std::vector< CPPL_INT > ia
rowIndex (NOT zero-based BUT one-based indexing)
Definition: dgrmatrix.hpp:12
CPPL_INT n
matrix column size
Definition: dgrmatrix.hpp:10
CPPL_INT m
matrix row size
Definition: dgrmatrix.hpp:9
std::vector< CPPL_INT > ja
columns (NOT zero-based BUT one-based indexing)
Definition: dgrmatrix.hpp:13
void dgrmatrix::clear ( )
inline

clear all the matrix data and set the sizes 0

Definition at line 3 of file dgrmatrix-misc.hpp.

References a, ia, ja, m, and n.

4 {CPPL_VERBOSE_REPORT;
5  m =0;
6  n =0;
7  a.clear();
8  ia.clear();
9  ja.clear();
10 }
std::vector< double > a
matrix component values
Definition: dgrmatrix.hpp:11
std::vector< CPPL_INT > ia
rowIndex (NOT zero-based BUT one-based indexing)
Definition: dgrmatrix.hpp:12
CPPL_INT n
matrix column size
Definition: dgrmatrix.hpp:10
CPPL_INT m
matrix row size
Definition: dgrmatrix.hpp:9
std::vector< CPPL_INT > ja
columns (NOT zero-based BUT one-based indexing)
Definition: dgrmatrix.hpp:13
dgrmatrix & dgrmatrix::zero ( )
inline

change the matrix into a zero matrix with no change in structure

Definition at line 14 of file dgrmatrix-misc.hpp.

References a.

15 {CPPL_VERBOSE_REPORT;
16  const std::vector<double>::iterator a_end =a.end();
17  for(std::vector<double>::iterator it=a.begin(); it!=a_end; it++){
18  (*it) =0.;
19  }
20  return *this;
21 }
std::vector< double > a
matrix component values
Definition: dgrmatrix.hpp:11
void dgrmatrix::copy ( const dgrmatrix mat)
inline

make a deep copy of the matrix

Definition at line 25 of file dgrmatrix-misc.hpp.

References a, ia, ja, m, and n.

Referenced by dgrmatrix().

26 {CPPL_VERBOSE_REPORT;
27  m =mat.m;
28  n =mat.n;
29  a =mat.a;
30  ia =mat.ia;
31  ja =mat.ja;
32 }
std::vector< double > a
matrix component values
Definition: dgrmatrix.hpp:11
std::vector< CPPL_INT > ia
rowIndex (NOT zero-based BUT one-based indexing)
Definition: dgrmatrix.hpp:12
CPPL_INT n
matrix column size
Definition: dgrmatrix.hpp:10
CPPL_INT m
matrix row size
Definition: dgrmatrix.hpp:9
std::vector< CPPL_INT > ja
columns (NOT zero-based BUT one-based indexing)
Definition: dgrmatrix.hpp:13
bool dgrmatrix::isListed ( const CPPL_INT &  i,
const CPPL_INT &  j 
) const
inline

check if the component is listed

Definition at line 36 of file dgrmatrix-misc.hpp.

References i(), ia, ja, m, and n.

37 {CPPL_VERBOSE_REPORT;
38 #ifdef CPPL_DEBUG
39  if( i<0 || j<0 || m<=i || n<=j ){
40  ERROR_REPORT;
41  std::cerr << "The required component is out of the matrix size." << std::endl
42  << "Your input is (" << i << "," << j << "), whereas the matrix size is " << m << "x" << n << "." << std::endl;
43  exit(1);
44  }
45 #endif//CPPL_DEBUG
46  //////// search ////////
47  int k_beg =ia[i]-1;
48  int k_end =ia[i+1]-1;
49  for(int k=k_beg; k<k_end; k++){
50  if(ja[k]==j+1){
51  return true;//found
52  }
53  }
54  //////// not found ////////
55  return false;
56 }
std::vector< CPPL_INT > ia
rowIndex (NOT zero-based BUT one-based indexing)
Definition: dgrmatrix.hpp:12
_dgematrix i(const _dgbmatrix &mat)
CPPL_INT n
matrix column size
Definition: dgrmatrix.hpp:10
CPPL_INT m
matrix row size
Definition: dgrmatrix.hpp:9
std::vector< CPPL_INT > ja
columns (NOT zero-based BUT one-based indexing)
Definition: dgrmatrix.hpp:13
void dgrmatrix::checkup ( )
inline

health checkup

Definition at line 64 of file dgrmatrix-misc.hpp.

References a, ia, ja, m, and n.

65 {CPPL_VERBOSE_REPORT;
66  //////// size check ////////
67  if(m<0){
68  ERROR_REPORT;
69  std::cerr << "m<0" << std::endl;
70  exit(1);
71  }
72  if(n<0){
73  ERROR_REPORT;
74  std::cerr << "n<0" << std::endl;
75  exit(1);
76  }
77  if(a.size()!=ja.size()){
78  ERROR_REPORT;
79  std::cerr << "a.size()!=ja.size()" << std::endl;
80  exit(1);
81  }
82  if(ia.size()!=size_t(m+1)){
83  ERROR_REPORT;
84  std::cerr << "ia.size()!=m+1" << std::endl;
85  exit(1);
86  }
87  if(a.size()>size_t(m*n)){
88  ERROR_REPORT;
89  std::cerr << "a.size()>m*n" << std::endl;
90  exit(1);
91  }
92 
93  //////// index check ////////
94  //not yet....
95 
96  std::cerr << "# [NOTE]@dgrmatrix::checkup(): This matrix is fine." << std::endl;
97 }
std::vector< double > a
matrix component values
Definition: dgrmatrix.hpp:11
std::vector< CPPL_INT > ia
rowIndex (NOT zero-based BUT one-based indexing)
Definition: dgrmatrix.hpp:12
CPPL_INT n
matrix column size
Definition: dgrmatrix.hpp:10
CPPL_INT m
matrix row size
Definition: dgrmatrix.hpp:9
std::vector< CPPL_INT > ja
columns (NOT zero-based BUT one-based indexing)
Definition: dgrmatrix.hpp:13
CPPL_INT dgrmatrix::pardiso ( dcovector b) const
inline

solve A*x=b for real and unsymmetric matrix using Intel PARDISO.
The argument is dcovector b. b is overwritten and become the solution x. A is not overwritten.

Definition at line 7 of file dgrmatrix-pardiso.hpp.

References a, dcovector::array, ia, ja, dcovector::l, m, n, and swap.

8 {CPPL_VERBOSE_REPORT;
9 #ifndef __INTEL_COMPILER
10  ERROR_REPORT;
11  std::cerr << "dgrmatrix::pardiso 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 
16 
17 #else //__INTEL_COMPILER is defined.
18  //#ifdef CPPL_DEBUG
19  if(m!=n || 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  //////// pardisoinit ////////
28  //std::cerr << "initializing" << std::endl;
29  _MKL_DSS_HANDLE_t pt[64];
30  MKL_INT mtype =11;//real unsymmetric
31  MKL_INT iparm[64];
32  PARDISOINIT(pt, &mtype, iparm);
33  iparm[1] =3;//parallel fill-in reducing ordering
34  iparm[23] =1;//use two-level scheduling factorization algorithm
35  iparm[26] =0;//disable matrix checker
36  iparm[34] =0;//use one-base array index
37 
38  //////// pardiso ////////
39  //std::cerr << "solving" << std::endl;
40  MKL_INT maxfct =1;
41  MKL_INT mnum =1;
42  MKL_INT phase =13;
43  MKL_INT MKL_INT_n =MKL_INT(n);
44  double* a0 = const_cast<double*>(&a[0]);
45  MKL_INT* ia0 = const_cast<MKL_INT*>(&ia[0]);
46  MKL_INT* ja0 = const_cast<MKL_INT*>(&ja[0]);
47  std::vector<MKL_INT> perm(n);
48  MKL_INT nrhs =1;
49  MKL_INT msglvl =0;
50  dcovector x(b.l);
51  MKL_INT error =1;
52  PARDISO(pt, &maxfct, &mnum, &mtype, &phase, &MKL_INT_n, a0, ia0, ja0, &perm[0], &nrhs, iparm, &msglvl, b.array, x.array, &error);
53  swap(b,x);//set b as x
54  if(error!=0){
55  WARNING_REPORT;
56  std::cerr << "Serious trouble happend. error = " << error << "." << std::endl;
57  }
58 
59  //////// release memory ////////
60  phase =-1;
61  MKL_INT error2 =1;
62  PARDISO(pt, &maxfct, &mnum, &mtype, &phase, &MKL_INT_n, a0, ia0, ja0, &perm[0], &nrhs, iparm, &msglvl, b.array, x.array, &error2);
63 
64  return error;
65 #endif //__INTEL_COMPILER
66 }
std::vector< double > a
matrix component values
Definition: dgrmatrix.hpp:11
friend void swap(dgrmatrix &, dgrmatrix &)
CPPL_INT l
vector size
Definition: dcovector.hpp:9
std::vector< CPPL_INT > ia
rowIndex (NOT zero-based BUT one-based indexing)
Definition: dgrmatrix.hpp:12
CPPL_INT n
matrix column size
Definition: dgrmatrix.hpp:10
double * array
1D array to store vector data
Definition: dcovector.hpp:11
Real Double-precision Column Vector Class.
Definition: dcovector.hpp:3
CPPL_INT m
matrix row size
Definition: dgrmatrix.hpp:9
std::vector< CPPL_INT > ja
columns (NOT zero-based BUT one-based indexing)
Definition: dgrmatrix.hpp:13
CPPL_INT dgrmatrix::dfgmres ( dcovector b,
const double  eps 
) const
inline

solve A*x=b for real and unsymmetric matrix using DFGMRES of Intel RCI ISS without preconditioning.
The argument is dcovector b. b is overwritten and become the solution x. A is not overwritten.

Definition at line 7 of file dgrmatrix-rci.hpp.

References a, dcovector::array, damax, ia, ja, dcovector::l, m, n, swap, and dcovector::zero().

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 }
std::vector< double > a
matrix component values
Definition: dgrmatrix.hpp:11
friend void swap(dgrmatrix &, dgrmatrix &)
CPPL_INT l
vector size
Definition: dcovector.hpp:9
std::vector< CPPL_INT > ia
rowIndex (NOT zero-based BUT one-based indexing)
Definition: dgrmatrix.hpp:12
dcovector & zero()
CPPL_INT n
matrix column size
Definition: dgrmatrix.hpp:10
double * array
1D array to store vector data
Definition: dcovector.hpp:11
Real Double-precision Column Vector Class.
Definition: dcovector.hpp:3
CPPL_INT m
matrix row size
Definition: dgrmatrix.hpp:9
std::vector< CPPL_INT > ja
columns (NOT zero-based BUT one-based indexing)
Definition: dgrmatrix.hpp:13
friend double damax(const dgrmatrix &)
CPPL_INT dgrmatrix::ilut_dfgmres ( dcovector b,
const int  maxfil,
const double  eps 
) const
inline

solve A*x=b for real and unsymmetric matrix using DFGMRES with ILUT precondition of Intel RCI ISS.
The argument is dcovector b. b is overwritten and become the solution x. A is not overwritten.

Definition at line 144 of file dgrmatrix-rci.hpp.

References a, dcovector::array, damax, ia, ja, dcovector::l, m, n, swap, and dcovector::zero().

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 }
std::vector< double > a
matrix component values
Definition: dgrmatrix.hpp:11
friend void swap(dgrmatrix &, dgrmatrix &)
CPPL_INT l
vector size
Definition: dcovector.hpp:9
std::vector< CPPL_INT > ia
rowIndex (NOT zero-based BUT one-based indexing)
Definition: dgrmatrix.hpp:12
dcovector & zero()
CPPL_INT n
matrix column size
Definition: dgrmatrix.hpp:10
double * array
1D array to store vector data
Definition: dcovector.hpp:11
Real Double-precision Column Vector Class.
Definition: dcovector.hpp:3
CPPL_INT m
matrix row size
Definition: dgrmatrix.hpp:9
std::vector< CPPL_INT > ja
columns (NOT zero-based BUT one-based indexing)
Definition: dgrmatrix.hpp:13
friend double damax(const dgrmatrix &)
dgrmatrix& dgrmatrix::operator= ( const dgrmatrix )
inline
dgrmatrix & dgrmatrix::operator*= ( const double &  d)
inline

dgrmatrix*=double operator

Definition at line 3 of file dgrmatrix-double.hpp.

References a.

4 {CPPL_VERBOSE_REPORT;
5  const std::vector<double>::iterator a_end =a.end();
6  for(std::vector<double>::iterator it=a.begin(); it!=a_end; it++){
7  (*it) *=d;
8  }
9 
10  return *this;
11 }
std::vector< double > a
matrix component values
Definition: dgrmatrix.hpp:11
dgrmatrix & dgrmatrix::operator/= ( const double &  d)
inline

dgrmatrix/=double operator

Definition at line 15 of file dgrmatrix-double.hpp.

References a.

16 {CPPL_VERBOSE_REPORT;
17  const std::vector<double>::iterator a_end =a.end();
18  for(std::vector<double>::iterator it=a.begin(); it!=a_end; it++){
19  (*it) /=d;
20  }
21 
22  return *this;
23 }
std::vector< double > a
matrix component values
Definition: dgrmatrix.hpp:11

Friends And Related Function Documentation

std::ostream& operator<< ( std::ostream &  s,
const dgrmatrix mat 
)
friend

Definition at line 60 of file dgrmatrix-io.hpp.

61 {CPPL_VERBOSE_REPORT;
62  for(CPPL_INT i=0; i<mat.m; i++){
63  int k_beg =mat.ia[i]-1;
64  int k_end =mat.ia[i+1]-1;
65  int j =0;
66  for(int k=k_beg; k<k_end; k++){
67  if(j<mat.ja[k]-1){
68  for(; j<mat.ja[k]-1; j++){
69  s << " x";
70  }
71  }
72  s << " " << mat.a[k];
73  j++;
74  }
75  for(; j<mat.n; j++){
76  s << " x";
77  }
78  s << std::endl;
79  }
80 
81  return s;
82 }
std::vector< double > a
matrix component values
Definition: dgrmatrix.hpp:11
std::vector< CPPL_INT > ia
rowIndex (NOT zero-based BUT one-based indexing)
Definition: dgrmatrix.hpp:12
_dgematrix i(const _dgbmatrix &mat)
CPPL_INT n
matrix column size
Definition: dgrmatrix.hpp:10
CPPL_INT m
matrix row size
Definition: dgrmatrix.hpp:9
std::vector< CPPL_INT > ja
columns (NOT zero-based BUT one-based indexing)
Definition: dgrmatrix.hpp:13
void swap ( dgrmatrix A,
dgrmatrix B 
)
friend

swap two matrices

Definition at line 105 of file dgrmatrix-misc.hpp.

Referenced by dfgmres(), ilut_dfgmres(), and pardiso().

106 {CPPL_VERBOSE_REPORT;
107  std::swap(A.n,B.n);
108  std::swap(A.m,B.m);
109  std::swap(A.a,B.a);
110  std::swap(A.ia,B.ia);
111  std::swap(A.ja,B.ja);
112 }
std::vector< double > a
matrix component values
Definition: dgrmatrix.hpp:11
void swap(dgrmatrix &A, dgrmatrix &B)
std::vector< CPPL_INT > ia
rowIndex (NOT zero-based BUT one-based indexing)
Definition: dgrmatrix.hpp:12
CPPL_INT n
matrix column size
Definition: dgrmatrix.hpp:10
CPPL_INT m
matrix row size
Definition: dgrmatrix.hpp:9
std::vector< CPPL_INT > ja
columns (NOT zero-based BUT one-based indexing)
Definition: dgrmatrix.hpp:13
double damax ( const dgrmatrix mat)
friend

return its largest absolute value

Definition at line 3 of file dgrmatrix-calc.hpp.

Referenced by dfgmres(), and ilut_dfgmres().

4 {CPPL_VERBOSE_REPORT;
5  //////// exception ////////
6  if(mat.a.size()==0){
7  return 0.;
8  }
9 
10  //////// find ////////
11  const size_t mat_a_size =mat.a.size();
12  double amax =0.;
13  double vmax;
14  for(size_t k=0; k<mat_a_size; k++){
15  if( amax < fabs(mat.a[k]) ){
16  amax =fabs(mat.a[k]);
17  vmax =mat.a[k];
18  }
19  }
20 
21  return vmax;
22 }
std::vector< double > a
matrix component values
Definition: dgrmatrix.hpp:11
dgrmatrix operator* ( const dgrmatrix mat,
const double &  d 
)
friend

dgrmatrix*double operator

Definition at line 31 of file dgrmatrix-double.hpp.

32 {CPPL_VERBOSE_REPORT;
33  dgrmatrix newmat =mat;
34 
35  const std::vector<double>::iterator newmat_a_end =newmat.a.end();
36  for(std::vector<double>::iterator it=newmat.a.begin(); it!=newmat_a_end; it++){
37  (*it) *=d;
38  }
39 
40  return newmat;
41 }
std::vector< double > a
matrix component values
Definition: dgrmatrix.hpp:11
Real Double-precision General Compressed Sparse Row (CSR) Matrix Class.
Definition: dgrmatrix.hpp:3
dgrmatrix operator* ( const double &  d,
const dgrmatrix mat 
)
friend

double*dgrmatrix operator

Definition at line 3 of file double-dgrmatrix.hpp.

4 {CPPL_VERBOSE_REPORT;
5  dgrmatrix newmat =mat;
6 
7  const size_t newmat_a_size =newmat.a.size();
8  for(size_t k=0; k<newmat_a_size; k++){
9  newmat.a[k] *= d;
10  }
11 
12  return newmat;
13 }
std::vector< double > a
matrix component values
Definition: dgrmatrix.hpp:11
Real Double-precision General Compressed Sparse Row (CSR) Matrix Class.
Definition: dgrmatrix.hpp:3
_dcovector operator* ( const dgrmatrix mat,
const dcovector vec 
)
friend

dgrmatrix*dcovector operator

Definition at line 3 of file dgrmatrix-dcovector.hpp.

4 {CPPL_VERBOSE_REPORT;
5 #ifdef CPPL_DEBUG
6  if(mat.n!=vec.l){
7  ERROR_REPORT;
8  std::cerr << "These matrix and vector can not make a product." << std::endl
9  << "Your input was (" << mat.m << "x" << mat.n << ") * (" << vec.l << ")." << std::endl;
10  exit(1);
11  }
12 #endif//CPPL_DEBUG
13 
14 #ifdef __INTEL_COMPILER
15  dcovector newvec(mat.m);
16  char transa ='N';
17  MKL_INT m =MKL_INT(mat.m);
18  double* a =const_cast<double*>(&mat.a[0]);
19  MKL_INT* ia =const_cast<MKL_INT*>(&mat.ia[0]);
20  MKL_INT* ja =const_cast<MKL_INT*>(&mat.ja[0]);
21  MKL_DCSRGEMV(&transa, &m, a, ia, ja, vec.array, newvec.array);
22  return _(newvec);
23 
24 
25 #else //__INTEL_COMPILER is not defined
26  dcovector newvec(mat.m);
27 #pragma omp parallel for
28  for(CPPL_INT i=0; i<mat.m; i++){
29  double sum =0.;
30  int k_beg =mat.ia[i]-1;
31  int k_end =mat.ia[i+1]-1;
32  for(int k=k_beg; k<k_end; k++){
33  int j =mat.ja[k]-1;
34  sum += mat.a[k] * vec(j);
35  }
36  newvec(i) =sum;
37  }
38  return _(newvec);
39 #endif//__INTEL_COMPILER
40 }
std::vector< double > a
matrix component values
Definition: dgrmatrix.hpp:11
CPPL_INT l
vector size
Definition: dcovector.hpp:9
std::vector< CPPL_INT > ia
rowIndex (NOT zero-based BUT one-based indexing)
Definition: dgrmatrix.hpp:12
_dgematrix i(const _dgbmatrix &mat)
CPPL_INT n
matrix column size
Definition: dgrmatrix.hpp:10
double * array
1D array to store vector data
Definition: dcovector.hpp:11
Real Double-precision Column Vector Class.
Definition: dcovector.hpp:3
CPPL_INT m
matrix row size
Definition: dgrmatrix.hpp:9
std::vector< CPPL_INT > ja
columns (NOT zero-based BUT one-based indexing)
Definition: dgrmatrix.hpp:13
_dcovector _(dcovector &vec)
_dcovector operator* ( const dgrmatrix mat,
const _dcovector vec 
)
friend

dgrmatrix*_dcovector operator

Definition at line 3 of file dgrmatrix-_dcovector.hpp.

4 {CPPL_VERBOSE_REPORT;
5 #ifdef CPPL_DEBUG
6  if(mat.n!=vec.l){
7  ERROR_REPORT;
8  std::cerr << "These matrix and vector can not make a product." << std::endl
9  << "Your input was (" << mat.m << "x" << mat.n << ") * (" << vec.l << ")." << std::endl;
10  exit(1);
11  }
12 #endif//CPPL_DEBUG
13 
14  dcovector VEC =vec;
15  dcovector newvec =mat*VEC;
16  return _(newvec);
17 }
CPPL_INT l
vector size
Definition: _dcovector.hpp:9
CPPL_INT n
matrix column size
Definition: dgrmatrix.hpp:10
Real Double-precision Column Vector Class.
Definition: dcovector.hpp:3
CPPL_INT m
matrix row size
Definition: dgrmatrix.hpp:9
_dcovector _(dcovector &vec)
dgrmatrix operator/ ( const dgrmatrix mat,
const double &  d 
)
friend

dgrmatrix/double operator

Definition at line 45 of file dgrmatrix-double.hpp.

46 {CPPL_VERBOSE_REPORT;
47  dgrmatrix newmat =mat;
48 
49  const std::vector<double>::iterator newmat_a_end =newmat.a.end();
50  for(std::vector<double>::iterator it=newmat.a.begin(); it!=newmat_a_end; it++){
51  (*it) /=d;
52  }
53 
54  return newmat;
55 }
std::vector< double > a
matrix component values
Definition: dgrmatrix.hpp:11
Real Double-precision General Compressed Sparse Row (CSR) Matrix Class.
Definition: dgrmatrix.hpp:3

Member Data Documentation

CPPL_INT dgrmatrix::m
CPPL_INT dgrmatrix::n
std::vector<double> dgrmatrix::a
std::vector<CPPL_INT> dgrmatrix::ia

rowIndex (NOT zero-based BUT one-based indexing)

Definition at line 12 of file dgrmatrix.hpp.

Referenced by checkup(), clear(), copy(), dfgmres(), ilut_dfgmres(), isListed(), operator()(), operator*(), operator<<(), pardiso(), read(), swap(), to_dgematrix(), dgsmatrix::to_dgrmatrix(), write(), and ~dgrmatrix().

std::vector<CPPL_INT> dgrmatrix::ja

columns (NOT zero-based BUT one-based indexing)

Definition at line 13 of file dgrmatrix.hpp.

Referenced by checkup(), clear(), copy(), dfgmres(), ilut_dfgmres(), isListed(), operator()(), operator*(), operator<<(), pardiso(), read(), swap(), to_dgematrix(), dgsmatrix::to_dgrmatrix(), write(), and ~dgrmatrix().


The documentation for this class was generated from the following files: