My Project
dsymatrix-lapack.hpp
1 //=============================================================================
6 inline long dsymatrix::dsysv(dgematrix& mat)
7 {
8 #ifdef CPPL_VERBOSE
9  std::cerr << "# [MARK] dsymatrix::dsysv(dgematrix&)"
10  << std::endl;
11 #endif//CPPL_VERBOSE
12 
13 #ifdef CPPL_DEBUG
14  if(N!=mat.N){
15  std::cerr << "[ERROR] dsymatrix::dsysv(dgematrix&) " << 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  double *WORK( new double[1] );
26  dsysv_(UPLO, N, NRHS, Array, LDA, IPIV, mat.Array, LDB, WORK, LWORK, INFO);
27 
28  INFO=1;
29  LWORK = long(WORK[0]);
30  delete [] WORK; WORK = new double[LWORK];
31  dsysv_(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] dsymatrix::dsysv(dsymatrix&) "
36  << "Serious trouble happend. INFO = " << INFO << "."
37  << std::endl;
38  }
39  return INFO;
40 }
41 
42 //=============================================================================
47 inline long dsymatrix::dsysv(dcovector& vec)
48 {
49 #ifdef CPPL_VERBOSE
50  std::cerr << "# [MARK] dsymatrix::dsysv(dcovector&)"
51  << std::endl;
52 #endif//CPPL_VERBOSE
53 
54 #ifdef CPPL_DEBUG
55  if(N!=vec.L){
56  std::cerr << "[ERROR] dsymatrix::dsysv(dcovector&) " << 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  double *WORK( new double[1] );
67  dsysv_(UPLO, N, NRHS, Array, LDA, IPIV, vec.Array, LDB, WORK, LWORK, INFO);
68 
69  INFO=1;
70  LWORK = long(WORK[0]);
71  delete WORK; WORK = new double[LWORK];
72  dsysv_(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] dsymatrix::dsysv(dcovector&) "
77  << "Serious trouble happend. INFO = " << INFO << "."
78  << std::endl;
79  }
80  return INFO;
81 }
82 
86 
87 //=============================================================================
94 inline long dsymatrix::dsyev(std::vector<double>& w, const bool& jobz=0)
95 {
96 #ifdef CPPL_VERBOSE
97  std::cerr << "# [MARK] dsymatrix::dsyev(std::vector<double>&, const bool&)"
98  << std::endl;
99 #endif//CPPL_VERBOSE
100 
101  w.resize(N);
102  char JOBZ, UPLO('L');
103  if(jobz==0){ JOBZ='N'; } else{ JOBZ='V'; }
104  long LDA(N), INFO(1), LWORK(-1);
105  double *WORK(new double[1]);
106  dsyev_(JOBZ, UPLO, N, Array, LDA, &w[0], WORK, LWORK, INFO);
107 
108  INFO=1;
109  LWORK = long(WORK[0]);
110  delete [] WORK; WORK = new double[LWORK];
111  dsyev_(JOBZ, UPLO, N, Array, LDA, &w[0], WORK, LWORK, INFO);
112  delete [] WORK;
113 
114  if(INFO!=0){
115  std::cerr << "[WARNING] dsymatrix::dsyev(vector<double>&, const bool&)"
116  << std::endl
117  << "Serious trouble happend. INFO = " << INFO << "."
118  << std::endl;
119  }
120  return INFO;
121 }
122 
123 //=============================================================================
130 inline long dsymatrix::dsyev(std::vector<double>& w, std::vector<dcovector>& v)
131 {
132 #ifdef CPPL_VERBOSE
133  std::cerr << "# [MARK] dsymatrix::dsyev(std::vector<double>&, std::vector<dcovector>&)"
134  << std::endl;
135 #endif//CPPL_VERBOSE
136 
137  w.resize(N); v.resize(N);
138  for(long i=0; i<N; i++){ v[i].resize(N); }
139  char JOBZ('V'), UPLO('L');
140  long LDA(N), INFO(1), LWORK(-1);
141  double *WORK(new double[1]);
142  dsyev_(JOBZ, UPLO, N, Array, LDA, &w[0], WORK, LWORK, INFO);
143 
144  INFO=1;
145  LWORK = long(WORK[0]);
146  delete [] WORK; WORK = new double[LWORK];
147  dsyev_(JOBZ, UPLO, N, Array, LDA, &w[0], WORK, LWORK, INFO);
148  delete [] WORK;
149 
151  for(long i=0; i<N; i++){ for(long j=0; j<N; j++){
152  v[j](i) = Array[i+N*j];
153  }}
154 
155  if(INFO!=0){
156  std::cerr << "[WARNING] dsymatrix::dsyev"
157  << "(vector<double>&, vector<dcovector>&)" << std::endl
158  << "Serious trouble happend. INFO = " << INFO << "."
159  << std::endl;
160  }
161  return INFO;
162 }
163 
164 //=============================================================================
171 inline long dsymatrix::dsyev(std::vector<double>& w, std::vector<drovector>& v)
172 {
173 #ifdef CPPL_VERBOSE
174  std::cerr << "# [MARK] dsymatrix::dsyev(std::vector<double>&, std::vector<drovector>&)"
175  << std::endl;
176 #endif//CPPL_VERBOSE
177 
178  w.resize(N); v.resize(N);
179  for(long i=0; i<N; i++){ v[i].resize(N); }
180  char JOBZ('V'), UPLO('L');
181  long LDA(N), INFO(1), LWORK(-1);
182  double *WORK(new double[1]);
183  dsyev_(JOBZ, UPLO, N, Array, LDA, &w[0], WORK, LWORK, INFO);
184 
185  INFO=1;
186  LWORK = long(WORK[0]);
187  delete [] WORK; WORK = new double[LWORK];
188  dsyev_(JOBZ, UPLO, N, Array, LDA, &w[0], WORK, LWORK, INFO);
189  delete [] WORK;
190 
192  for(long i=0; i<N; i++){ for(long j=0; j<N; j++){
193  v[j](i) = Array[i+N*j];
194  }}
195 
196  if(INFO!=0){
197  std::cerr << "[WARNING] dsymatrix::dsyev"
198  << "(vector<double>&, vector<dcovector>&)" << std::endl
199  << "Serious trouble happend. INFO = " << INFO << "."
200  << std::endl;
201  }
202  return INFO;
203 }
204 
208 
209 //=============================================================================
214 inline long dsymatrix::dsygv(dsymatrix& matB, std::vector<double>& w)
215 {
216 #ifdef CPPL_VERBOSE
217  std::cerr << "# [MARK] dsymatrix::dsygv(dsymatrix&, std::vector<double>&)"
218  << std::endl;
219 #endif//CPPL_VERBOSE
220 
221 #ifdef CPPL_DEBUG
222  if(matB.N!=N){
223  std::cerr << "[ERROR] dsymatrix::dsygv"
224  << "(dsymatrix&, vector<double>&) " << std::endl
225  << "The matrix B is not a matrix "
226  << "having the same size as \"this\" matrix." << std::endl
227  << "The B matrix is (" << matB.N << "x" << matB.N << ")."
228  << std::endl;
229  exit(1);
230  }
231 #endif//CPPL_DEBUG
232 
233  w.resize(N);
234  char JOBZ('N'), UPLO('L');
235  long ITYPE(1), LDA(N), LDB(N), LWORK(-1), INFO(1);
236  double *WORK(new double[1]);
237  dsygv_(ITYPE, JOBZ, UPLO, N, Array, LDA, matB.Array, LDB, &w[0],
238  WORK, LWORK, INFO);
239  INFO=1;
240  LWORK = long(WORK[0]);
241  delete [] WORK; WORK = new double[LWORK];
242  dsygv_(ITYPE, JOBZ, UPLO, N, Array, LDA, matB.Array, LDB, &w[0],
243  WORK, LWORK, INFO);
244  delete [] WORK;
245 
246  if(INFO!=0){
247  std::cerr << "[WARNING] dsymatrix::dsygv"
248  << "(dsymatrix&, vector<double>&)"
249  << std::endl
250  << "Serious trouble happend. INFO = " << INFO << "."
251  << std::endl
252  << "Make sure that B matrix is a symmetric-definite."
253  << std::endl;
254  }
255  return INFO;
256 }
257 
258 //=============================================================================
263 inline long dsymatrix::dsygv(dsymatrix& matB, std::vector<double>& w,
264  std::vector<dcovector>& v)
265 {
266 #ifdef CPPL_VERBOSE
267  std::cerr << "# [MARK] dsymatrix::dsygv(dsymatrix&, std::vector<double>&, std::vector<dcovector>&)"
268  << std::endl;
269 #endif//CPPL_VERBOSE
270 
271 #ifdef CPPL_DEBUG
272  if(matB.N!=N){
273  std::cerr << "[ERROR] dsymatrix::dsygv"
274  << "(dsymatrix&, vector<double>&) " << std::endl
275  << "The matrix B is not a matrix "
276  << "having the same size as \"this\" matrix." << std::endl
277  << "The B matrix is (" << matB.N << "x" << matB.N << ")."
278  << std::endl;
279  exit(1);
280  }
281 #endif//CPPL_DEBUG
282 
283  w.resize(N);
284  v.resize(N);
285  char JOBZ('V'), UPLO('L');
286  long ITYPE(1), LDA(N), LDB(N), LWORK(-1), INFO(1);
287  double *WORK(new double[1]);
288  dsygv_(ITYPE, JOBZ, UPLO, N, Array, LDA, matB.Array, LDB, &w[0],
289  WORK, LWORK, INFO);
290  INFO=1;
291  LWORK = long(WORK[0]);
292  std::cout << " LWORK = " << LWORK <<std::endl;
293  delete [] WORK; WORK = new double[LWORK];
294  dsygv_(ITYPE, JOBZ, UPLO, N, Array, LDA, matB.Array, LDB, &w[0],
295  WORK, LWORK, INFO);
296  delete [] WORK;
297 
299  for(int i=0; i<N; i++){
300  v[i].resize(N);
301  for(int j=0; j<N; j++){
302  v[i](j) =Darray[i][j];
303  }
304  }
305 
306  if(INFO!=0){
307  std::cerr << "[WARNING] dsymatrix::dsygv"
308  << "(dsymatrix&, vector<double>&)"
309  << std::endl
310  << "Serious trouble happend. INFO = " << INFO << "."
311  << std::endl
312  << "Make sure that B matrix is a symmetric-definite."
313  << std::endl;
314  }
315  return INFO;
316 }
long dsysv(dgematrix &)
Definition: dsymatrix-lapack.hpp:6
friend _dgematrix i(const dsymatrix &)
Definition: dsymatrix-calc.hpp:22
long dsyev(std::vector< double > &, const bool &)
Definition: dsymatrix-lapack.hpp:94
double * Array
1D Array to store vector data
Definition: _drovector.hpp:8
Real Double-precision General Dence Matrix Class.
Definition: dgematrix.hpp:3
Real Double-precision Symmetric Matrix Class [L-type (UPLO=L) Strage].
Definition: dsymatrix.hpp:3
long dsygv(dsymatrix &, std::vector< double > &)
Definition: dsymatrix-lapack.hpp:214
Real Double-precision Column Vector Class.
Definition: dcovector.hpp:3