CPPLapack
 All Classes Files Functions Variables Friends Pages
zhematrix-lapack.hpp
Go to the documentation of this file.
1 //=============================================================================
2 /*! solve A*X=Y using zhesv\n
3  The argument is dmatrix Y. Y is overwritten and become the solution X.
4  A is also overwritten.
5 */
6 inline CPPL_INT zhematrix::zhesv(zgematrix& mat)
7 {CPPL_VERBOSE_REPORT;
8 #ifdef CPPL_DEBUG
9  if(n!=mat.n){
10  ERROR_REPORT;
11  std::cerr << "These two matrices cannot be solved." << std::endl
12  << "Your input was (" << n << "x" << n << ") and (" << mat.n << "x" << mat.n << ")." << std::endl;
13  exit(1);
14  }
15 #endif//CPPL_DEBUG
16 
17  char UPLO('l');
18  CPPL_INT NRHS(mat.n), LDA(n), *IPIV(new CPPL_INT[n]), LDB(mat.n), LWORK(-1), INFO(1);
19  comple *WORK(new comple[1]);
20  zhesv_(&UPLO, &n, &NRHS, array, &LDA, IPIV, mat.array, &LDB, WORK, &LWORK, &INFO);
21 
22  INFO=1;
23  LWORK = CPPL_INT(std::real(WORK[0]));
24  delete [] WORK;
25  WORK =new comple[LWORK];
26  zhesv_(&UPLO, &n, &NRHS, array, &LDA, IPIV, mat.array, &LDB, WORK, &LWORK, &INFO);
27 
28  delete [] WORK;
29  delete [] IPIV;
30 
31  if(INFO!=0){
32  WARNING_REPORT;
33  std::cerr << "Serious trouble happend. INFO = " << INFO << "." << std::endl;
34  }
35  return INFO;
36 }
37 
38 //=============================================================================
39 /*! solve A*x=y using zhesv\n
40  The argument is zcovector y. y is overwritten and become the solution x.
41  A is also overwritten.
42 */
43 inline CPPL_INT zhematrix::zhesv(zcovector& vec)
44 {CPPL_VERBOSE_REPORT;
45 #ifdef CPPL_DEBUG
46  if(n!=vec.l){
47  ERROR_REPORT;
48  std::cerr << "These matrix and vector cannot be solved." << std::endl
49  << "Your input was (" << n << "x" << n << ") and (" << vec.l << ")." << std::endl;
50  exit(1);
51  }
52 #endif//CPPL_DEBUG
53 
54  char UPLO('l');
55  CPPL_INT NRHS(1), LDA(n), *IPIV(new CPPL_INT[n]), LDB(vec.l), LWORK(-1), INFO(1);
56  comple *WORK( new comple[1] );
57  zhesv_(&UPLO, &n, &NRHS, array, &LDA, IPIV, vec.array, &LDB, WORK, &LWORK, &INFO);
58 
59  INFO=1;
60  LWORK = CPPL_INT(std::real(WORK[0]));
61  delete [] WORK;
62  WORK = new comple[LWORK];
63  zhesv_(&UPLO, &n, &NRHS, array, &LDA, IPIV, vec.array, &LDB, WORK, &LWORK, &INFO);
64 
65  delete [] WORK;
66  delete [] IPIV;
67 
68  if(INFO!=0){
69  WARNING_REPORT;
70  std::cerr << "Serious trouble happend. INFO = " << INFO << "." << std::endl;
71  }
72  return INFO;
73 }
74 
75 ///////////////////////////////////////////////////////////////////////////////
76 ///////////////////////////////////////////////////////////////////////////////
77 ///////////////////////////////////////////////////////////////////////////////
78 
79 //=============================================================================
80 /*! calculate eigenvalues and eigenvectors.\n
81  All of the arguments need not to be initialized.
82  w is overwitten and become eigenvalues.
83  This matrix is also overwritten.
84  if jobz=1, this matrix becomes eigenvectors.
85 */
86 inline CPPL_INT zhematrix::zheev(std::vector<double>& w,
87  const bool& jobz=0)
88 {CPPL_VERBOSE_REPORT;
89  w.resize(n);
90 
91  char JOBZ, UPLO('l');
92  if(jobz==0){ JOBZ='n'; } else{ JOBZ='V'; }
93  CPPL_INT LDA(n), INFO(1), LWORK(-1);
94  double *RWORK(new double[std::max(CPPL_INT(1), 3*n-2)]);
95  comple *WORK(new comple[1]);
96  zheev_(&JOBZ, &UPLO, &n, array, &LDA, &w[0], WORK, &LWORK, RWORK, &INFO);
97 
98  INFO=1;
99  LWORK = CPPL_INT(std::real(WORK[0]));
100  delete [] WORK;
101  WORK = new comple[LWORK];
102  zheev_(&JOBZ, &UPLO, &n, array, &LDA, &w[0], WORK, &LWORK, RWORK, &INFO);
103 
104  delete [] RWORK;
105  delete [] WORK;
106 
107  if(INFO!=0){
108  WARNING_REPORT;
109  std::cerr << "Serious trouble happend. INFO = " << INFO << "." << std::endl;
110  }
111  return INFO;
112 }
113 
114 //=============================================================================
115 /*! calculate eigenvalues and eigenvectors.\n
116  All of the arguments need not to be initialized.
117  w and v are overwitten and become
118  eigenvalues and eigenvectors, respectively.
119  This matrix is also overwritten.
120 */
121 inline CPPL_INT zhematrix::zheev(std::vector<double>& w,
122  std::vector<zcovector>& v)
123 {CPPL_VERBOSE_REPORT;
124  w.resize(n);
125  v.resize(n);
126  for(CPPL_INT i=0; i<n; i++){
127  v[i].resize(n);
128  }
129 
130  char JOBZ('V'), UPLO('l');
131  CPPL_INT LDA(n), INFO(1), LWORK(-1);
132  double *RWORK(new double[std::max(CPPL_INT(1), 3*n-2)]);
133  comple *WORK(new comple[1]);
134  zheev_(&JOBZ, &UPLO, &n, array, &LDA, &w[0], WORK, &LWORK, RWORK, &INFO);
135 
136  INFO=1;
137  LWORK = CPPL_INT(std::real(WORK[0]));
138  delete [] WORK;
139  WORK = new comple[LWORK];
140  zheev_(&JOBZ, &UPLO, &n, array, &LDA, &w[0], WORK, &LWORK, RWORK, &INFO);
141 
142  delete [] RWORK;
143  delete [] WORK;
144 
145  //// forming ////
146  for(CPPL_INT i=0; i<n; i++){
147  for(CPPL_INT j=0; j<n; j++){
148  v[j](i) = array[i+n*j];
149  }
150  }
151 
152  if(INFO!=0){
153  WARNING_REPORT;
154  std::cerr << "Serious trouble happend. INFO = " << INFO << "." << std::endl;
155  }
156  return INFO;
157 }
158 
159 //=============================================================================
160 /*! calculate eigenvalues and eigenvectors.\n
161  All of the arguments need not to be initialized.
162  w and v are overwitten and become
163  eigenvalues and eigenvectors, respectively.
164  This matrix is also overwritten.
165 */
166 inline CPPL_INT zhematrix::zheev(std::vector<double>& w,
167  std::vector<zrovector>& v)
168 {CPPL_VERBOSE_REPORT;
169  w.resize(n);
170  v.resize(n);
171  for(CPPL_INT i=0; i<n; i++){
172  v[i].resize(n);
173  }
174 
175  char JOBZ('V'), UPLO('l');
176  CPPL_INT LDA(n), INFO(1), LWORK(-1);
177  double *RWORK(new double[std::max(CPPL_INT(1), 3*n-2)]);
178  comple *WORK(new comple[1]);
179  zheev_(&JOBZ, &UPLO, &n, array, &LDA, &w[0], WORK, &LWORK, RWORK, &INFO);
180 
181  INFO=1;
182  LWORK = CPPL_INT(std::real(WORK[0]));
183  delete [] WORK;
184  WORK = new comple[LWORK];
185  zheev_(&JOBZ, &UPLO, &n, array, &LDA, &w[0], WORK, &LWORK, RWORK, &INFO);
186 
187  delete [] RWORK;
188  delete [] WORK;
189 
190  //// forming ////
191  for(CPPL_INT i=0; i<n; i++){
192  for(CPPL_INT j=0; j<n; j++){
193  v[j](i) = array[i+n*j];
194  }
195  }
196 
197  if(INFO!=0){
198  WARNING_REPORT;
199  std::cerr << "Serious trouble happend. INFO = " << INFO << "." << std::endl;
200  }
201  return INFO;
202 }
comple * array
1D array to store vector data
Definition: zcovector.hpp:10
CPPL_INT l
vector size
Definition: zcovector.hpp:9
CPPL_INT zhesv(zgematrix &)
CPPL_INT n
matrix column size
Definition: zgematrix.hpp:10
Complex Double-precision General Dence Matrix Class.
Definition: zgematrix.hpp:3
CPPL_INT n
matrix column size
Definition: zhematrix.hpp:11
comple * array
1D array to store matrix data
Definition: zgematrix.hpp:11
friend _zgematrix i(const zhematrix &)
comple * array
1D array to store matrix data
Definition: zhematrix.hpp:12
CPPL_INT zheev(std::vector< double > &, const bool &)
Complex Double-precision Column Vector Class.
Definition: zcovector.hpp:3