My Project
zhematrix-lapack.hpp
1 //=============================================================================
6 inline long zhematrix::zhesv(zgematrix& mat)
7 {
8 #ifdef CPPL_VERBOSE
9  std::cerr << "# [MARK] zhematrix::zhesv(zgematrix&)"
10  << std::endl;
11 #endif//CPPL_VERBOSE
12 
13 #ifdef CPPL_DEBUG
14  if(N!=mat.N){
15  std::cerr << "[ERROR] zhematrix::zhesv(zgematrix&) " << std::endl
16  << "These two matrices cannot be solved." << std::endl
17  << "Your input was (" << N << "x" << N << ") and ("
18  << mat.N << "x" << mat.N << ")." << std::endl;
19  exit(1);
20  }
21 #endif//CPPL_DEBUG
22 
23  char UPLO('L');
24  long NRHS(mat.N), LDA(N), *IPIV(new long[N]), LDB(mat.N), LWORK(-1), INFO(1);
25  std::complex<double> *WORK(new std::complex<double>[1]);
26  zhesv_(UPLO, N, NRHS, Array, LDA, IPIV, mat.Array, LDB, WORK, LWORK, INFO);
27 
28  INFO=1;
29  LWORK = long(std::real(WORK[0]));
30  delete [] WORK; WORK =new std::complex<double>[LWORK];
31  zhesv_(UPLO, N, NRHS, Array, LDA, IPIV, mat.Array, LDB, WORK, LWORK, INFO);
32  delete [] WORK; delete [] IPIV;
33 
34  if(INFO!=0){
35  std::cerr << "[WARNING] zhematrix::zhesv(zhematrix&) "
36  << "Serious trouble happend. INFO = " << INFO << "."
37  << std::endl;
38  }
39  return INFO;
40 }
41 
42 //=============================================================================
47 inline long zhematrix::zhesv(zcovector& vec)
48 {
49 #ifdef CPPL_VERBOSE
50  std::cerr << "# [MARK] zhematrix::zhesv(zcovector&)"
51  << std::endl;
52 #endif//CPPL_VERBOSE
53 
54 #ifdef CPPL_DEBUG
55  if(N!=vec.L){
56  std::cerr << "[ERROR] zhematrix::zhesv(zcovector&) " << std::endl
57  << "These matrix and vector cannot be solved." << std::endl
58  << "Your input was (" << N << "x" << N << ") and ("
59  << vec.L << ")." << std::endl;
60  exit(1);
61  }
62 #endif//CPPL_DEBUG
63 
64  char UPLO('L');
65  long NRHS(1), LDA(N), *IPIV(new long[N]), LDB(vec.L), LWORK(-1), INFO(1);
66  std::complex<double> *WORK( new std::complex<double>[1] );
67  zhesv_(UPLO, N, NRHS, Array, LDA, IPIV, vec.Array, LDB, WORK, LWORK, INFO);
68 
69  INFO=1;
70  LWORK = long(std::real(WORK[0]));
71  delete WORK; WORK = new std::complex<double>[LWORK];
72  zhesv_(UPLO, N, NRHS, Array, LDA, IPIV, vec.Array, LDB, WORK, LWORK, INFO);
73  delete [] WORK; delete [] IPIV;
74 
75  if(INFO!=0){
76  std::cerr << "[WARNING] zhematrix::zhesv(zcovector&) "
77  << "Serious trouble happend. INFO = " << INFO << "."
78  << std::endl;
79  }
80  return INFO;
81 }
82 
86 
87 //=============================================================================
94 inline long zhematrix::zheev(std::vector<double>& w,
95  const bool& jobz=0)
96 {
97 #ifdef CPPL_VERBOSE
98  std::cerr << "# [MARK] zhematrix::zheev(std::vector<double>&, const bool&)"
99  << std::endl;
100 #endif//CPPL_VERBOSE
101 
102  w.resize(N);
103  char JOBZ, UPLO('L');
104  if(jobz==0){ JOBZ='N'; } else{ JOBZ='V'; }
105  long LDA(N), INFO(1), LWORK(-1);
106  double *RWORK(new double[max(1, 3*N-2)]);
107  std::complex<double> *WORK(new std::complex<double>[1]);
108  zheev_(JOBZ, UPLO, N, Array, LDA, &w[0], WORK, LWORK, RWORK, INFO);
109 
110  INFO=1;
111  LWORK = long(std::real(WORK[0]));
112  delete [] WORK; WORK = new std::complex<double>[LWORK];
113  zheev_(JOBZ, UPLO, N, Array, LDA, &w[0], WORK, LWORK, RWORK, INFO);
114  delete [] RWORK; delete [] WORK;
115 
116  if(INFO!=0){
117  std::cerr << "[WARNING] zhematrix::zheev"
118  << "(vector<std::complex<double>>&, const bool&) "
119  << "Serious trouble happend. INFO = " << INFO << "."
120  << std::endl;
121  }
122  return INFO;
123 }
124 
125 //=============================================================================
132 inline long zhematrix::zheev(std::vector<double>& w,
133  std::vector<zcovector>& v)
134 {
135 #ifdef CPPL_VERBOSE
136  std::cerr << "# [MARK] zhematrix::zheev(std::vector<double>&, std::vector<zcovector>&)"
137  << std::endl;
138 #endif//CPPL_VERBOSE
139 
140  w.resize(N); v.resize(N);
141  for(long i=0; i<N; i++){ v[i].resize(N); }
142  char JOBZ('V'), UPLO('L');
143  long LDA(N), INFO(1), LWORK(-1);
144  double *RWORK(new double[max(1, 3*N-2)]);
145  std::complex<double> *WORK(new std::complex<double>[1]);
146  zheev_(JOBZ, UPLO, N, Array, LDA, &w[0], WORK, LWORK, RWORK, INFO);
147 
148  INFO=1;
149  LWORK = long(std::real(WORK[0]));
150  delete [] WORK; WORK = new std::complex<double>[LWORK];
151  zheev_(JOBZ, UPLO, N, Array, LDA, &w[0], WORK, LWORK, RWORK, INFO);
152  delete [] RWORK; delete [] WORK;
153 
155  for(long i=0; i<N; i++){ for(long j=0; j<N; j++){
156  v[j](i) = Array[i+N*j];
157  }}
158 
159  if(INFO!=0){
160  std::cerr << "[WARNING] zhematrix::zheev"
161  << "(vector<std::complex<double>>&, vector<zcovector>&) "
162  << "Serious trouble happend. INFO = " << INFO << "."
163  << std::endl;
164  }
165  return INFO;
166 }
167 
168 //=============================================================================
175 inline long zhematrix::zheev(std::vector<double>& w,
176  std::vector<zrovector>& v)
177 {
178 #ifdef CPPL_VERBOSE
179  std::cerr << "# [MARK] zhematrix::zheev(std::vector<double>&, std::vector<zrovector>&)"
180  << std::endl;
181 #endif//CPPL_VERBOSE
182 
183  w.resize(N); v.resize(N);
184  for(long i=0; i<N; i++){ v[i].resize(N); }
185  char JOBZ('V'), UPLO('L');
186  long LDA(N), INFO(1), LWORK(-1);
187  double *RWORK(new double[max(1, 3*N-2)]);
188  std::complex<double> *WORK(new std::complex<double>[1]);
189  zheev_(JOBZ, UPLO, N, Array, LDA, &w[0], WORK, LWORK, RWORK, INFO);
190 
191  INFO=1;
192  LWORK = long(std::real(WORK[0]));
193  delete [] WORK; WORK = new std::complex<double>[LWORK];
194  zheev_(JOBZ, UPLO, N, Array, LDA, &w[0], WORK, LWORK, RWORK, INFO);
195  delete [] RWORK; delete [] WORK;
196 
198  for(long i=0; i<N; i++){ for(long j=0; j<N; j++){
199  v[j](i) = Array[i+N*j];
200  }}
201 
202  if(INFO!=0){
203  std::cerr << "[WARNING] zhematrix::zheev"
204  << "(vector<std::complex<double>>&, vector<zcovector>&)"
205  << "Serious trouble happend. INFO = " << INFO << "."
206  << std::endl;
207  }
208  return INFO;
209 }
long zheev(std::vector< double > &, const bool &)
Definition: zhematrix-lapack.hpp:94
Complex Double-precision General Dence Matrix Class.
Definition: zgematrix.hpp:3
friend _zgematrix i(const zhematrix &)
Definition: zhematrix-calc.hpp:20
Complex Double-precision Column Vector Class.
Definition: zcovector.hpp:3
long zhesv(zgematrix &)
Definition: zhematrix-lapack.hpp:6
std::complex< double > * Array
1D Array to store vector data
Definition: _zrovector.hpp:8