CPPLapack
 All Classes Files Functions Variables Friends Pages
dsymatrix-lapack.hpp
Go to the documentation of this file.
1 //=============================================================================
2 /*! solve A*X=Y using dsysv\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 dsymatrix::dsysv(dgematrix& 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.m), LWORK(-1), INFO(1);
19  double *WORK( new double[1] );
20  dsysv_(&UPLO, &n, &NRHS, array, &LDA, IPIV, mat.array, &LDB, WORK, &LWORK, &INFO);
21 
22  INFO=1;
23  LWORK = CPPL_INT(WORK[0]);
24  delete [] WORK;
25  WORK = new double[LWORK];
26  dsysv_(&UPLO, &n, &NRHS, array, &LDA, IPIV, mat.array, &LDB, WORK, &LWORK, &INFO);
27  delete [] WORK;
28  delete [] IPIV;
29 
30  if(INFO!=0){
31  WARNING_REPORT;
32  std::cerr << "Serious trouble happend. INFO = " << INFO << "." << std::endl;
33  }
34  return INFO;
35 }
36 
37 //=============================================================================
38 /*! solve A*x=y using dsysv\n
39  The argument is dcovector y. y is overwritten and become the solution x.
40  A is also overwritten.
41 */
42 inline CPPL_INT dsymatrix::dsysv(dcovector& vec)
43 {CPPL_VERBOSE_REPORT;
44 #ifdef CPPL_DEBUG
45  if(n!=vec.l){
46  ERROR_REPORT;
47  std::cerr << "These matrix and vector cannot be solved." << std::endl
48  << "Your input was (" << n << "x" << n << ") and (" << vec.l << ")." << std::endl;
49  exit(1);
50  }
51 #endif//CPPL_DEBUG
52 
53  char UPLO('l');
54  CPPL_INT NRHS(1), LDA(n), *IPIV(new CPPL_INT[n]), LDB(vec.l), LWORK(-1), INFO(1);
55  double *WORK( new double[1] );
56  dsysv_(&UPLO, &n, &NRHS, array, &LDA, IPIV, vec.array, &LDB, WORK, &LWORK, &INFO);
57 
58  INFO=1;
59  LWORK = CPPL_INT(WORK[0]);
60  delete [] WORK;
61  WORK = new double[LWORK];
62  dsysv_(&UPLO, &n, &NRHS, array, &LDA, IPIV, vec.array, &LDB, WORK, &LWORK, &INFO);
63  delete [] WORK;
64  delete [] IPIV;
65 
66  if(INFO!=0){
67  WARNING_REPORT;
68  std::cerr << "Serious trouble happend. INFO = " << INFO << "." << std::endl;
69  }
70  return INFO;
71 }
72 
73 ///////////////////////////////////////////////////////////////////////////////
74 ///////////////////////////////////////////////////////////////////////////////
75 ///////////////////////////////////////////////////////////////////////////////
76 
77 //=============================================================================
78 /*! calculate eigenvalues and eigenvectors.\n
79  All of the arguments need not to be initialized.
80  w is overwitten and become eigenvalues.
81  This matrix is also overwritten.
82  if jobz=1, this matrix becomes eigenvectors.
83 */
84 inline CPPL_INT dsymatrix::dsyev(std::vector<double>& w, const bool& jobz=0)
85 {CPPL_VERBOSE_REPORT;
86  w.resize(n);
87  char JOBZ, UPLO('l');
88  if(jobz==0){ JOBZ='n'; } else{ JOBZ='V'; }
89  CPPL_INT LDA(n), INFO(1), LWORK(-1);
90  double *WORK(new double[1]);
91  dsyev_(&JOBZ, &UPLO, &n, array, &LDA, &w[0], WORK, &LWORK, &INFO);
92 
93  INFO=1;
94  LWORK = CPPL_INT(WORK[0]);
95  delete [] WORK;
96  WORK = new double[LWORK];
97  dsyev_(&JOBZ, &UPLO, &n, array, &LDA, &w[0], WORK, &LWORK, &INFO);
98  delete [] WORK;
99 
100  if(INFO!=0){
101  WARNING_REPORT;
102  std::cerr << "Serious trouble happend. INFO = " << INFO << "." << std::endl;
103  }
104  return INFO;
105 }
106 
107 //=============================================================================
108 /*! calculate eigenvalues and eigenvectors.\n
109  All of the arguments need not to be initialized.
110  w and v are overwitten and become
111  eigenvalues and eigenvectors, respectively.
112  This matrix is also overwritten.
113 */
114 inline CPPL_INT dsymatrix::dsyev(std::vector<double>& w, std::vector<dcovector>& v)
115 {CPPL_VERBOSE_REPORT;
116  w.resize(n);
117  v.resize(n);
118  for(CPPL_INT i=0; i<n; i++){
119  v[i].resize(n);
120  }
121 
122  char JOBZ('V'), UPLO('l');
123  CPPL_INT LDA(n), INFO(1), LWORK(-1);
124  double *WORK(new double[1]);
125  dsyev_(&JOBZ, &UPLO, &n, array, &LDA, &w[0], WORK, &LWORK, &INFO);
126 
127  INFO=1;
128  LWORK = CPPL_INT(WORK[0]);
129  delete [] WORK;
130  WORK = new double[LWORK];
131  dsyev_(&JOBZ, &UPLO, &n, array, &LDA, &w[0], WORK, &LWORK, &INFO);
132  delete [] WORK;
133 
134  //// forming ////
135  for(CPPL_INT j=0; j<n; j++){
136  for(CPPL_INT i=0; i<n; i++){
137  v[j](i) = array[i+n*j];
138  }
139  }
140 
141  if(INFO!=0){
142  WARNING_REPORT;
143  std::cerr << "Serious trouble happend. INFO = " << INFO << "." << std::endl;
144  }
145  return INFO;
146 }
147 
148 //=============================================================================
149 /*! calculate eigenvalues and eigenvectors.\n
150  All of the arguments need not to be initialized.
151  w and v are overwitten and become
152  eigenvalues and eigenvectors, respectively.
153  This matrix is also overwritten.
154 */
155 inline CPPL_INT dsymatrix::dsyev(std::vector<double>& w, std::vector<drovector>& v)
156 {CPPL_VERBOSE_REPORT;
157  w.resize(n);
158  v.resize(n);
159  for(CPPL_INT i=0; i<n; i++){
160  v[i].resize(n);
161  }
162 
163  char JOBZ('V'), UPLO('l');
164  CPPL_INT LDA(n), INFO(1), LWORK(-1);
165  double *WORK(new double[1]);
166  dsyev_(&JOBZ, &UPLO, &n, array, &LDA, &w[0], WORK, &LWORK, &INFO);
167 
168  INFO=1;
169  LWORK = CPPL_INT(WORK[0]);
170  delete [] WORK;
171  WORK = new double[LWORK];
172  dsyev_(&JOBZ, &UPLO, &n, array, &LDA, &w[0], WORK, &LWORK, &INFO);
173  delete [] WORK;
174 
175  //// forming ////
176  for(CPPL_INT j=0; j<n; j++){
177  for(CPPL_INT i=0; i<n; i++){
178  v[j](i) = array[i+n*j];
179  }
180  }
181 
182  if(INFO!=0){
183  WARNING_REPORT;
184  std::cerr << "Serious trouble happend. INFO = " << INFO << "." << std::endl;
185  }
186  return INFO;
187 }
188 
189 ///////////////////////////////////////////////////////////////////////////////
190 ///////////////////////////////////////////////////////////////////////////////
191 ///////////////////////////////////////////////////////////////////////////////
192 
193 //=============================================================================
194 /*! calculate generalized eigenvalues\n
195  w is overwitten and become generalized eigenvalues.
196  This matrix and matB are also overwritten.
197 */
198 inline CPPL_INT dsymatrix::dsygv(dsymatrix& matB, std::vector<double>& w)
199 {CPPL_VERBOSE_REPORT;
200 #ifdef CPPL_DEBUG
201  if(matB.n!=n){
202  ERROR_REPORT;
203  std::cerr << "The matrix B is not a matrix having the same size as \"this\" matrix." << std::endl
204  << "The B matrix is (" << matB.n << "x" << matB.n << ")." << std::endl;
205  exit(1);
206  }
207 #endif//CPPL_DEBUG
208 
209  w.resize(n);
210  char JOBZ('n'), UPLO('l');
211  CPPL_INT ITYPE(1), LDA(n), LDB(n), LWORK(-1), INFO(1);
212  double *WORK(new double[1]);
213  dsygv_(&ITYPE, &JOBZ, &UPLO, &n, array, &LDA, matB.array, &LDB, &w[0], WORK, &LWORK, &INFO);
214 
215  INFO=1;
216  LWORK = CPPL_INT(WORK[0]);
217  delete [] WORK;
218  WORK = new double[LWORK];
219  dsygv_(&ITYPE, &JOBZ, &UPLO, &n, array, &LDA, matB.array, &LDB, &w[0], WORK, &LWORK, &INFO);
220  delete [] WORK;
221 
222  if(INFO!=0){
223  WARNING_REPORT;
224  std::cerr << "Serious trouble happend. INFO = " << INFO << "." << std::endl;
225  }
226  return INFO;
227 }
228 
229 //=============================================================================
230 /*! calculate generalized eigenvalues\n
231  w is overwitten and become generalized eigenvalues.
232  This matrix and matB are also overwritten.
233 */
234 inline CPPL_INT dsymatrix::dsygv(dsymatrix& matB, std::vector<double>& w,
235  std::vector<dcovector>& v)
236 {CPPL_VERBOSE_REPORT;
237 #ifdef CPPL_DEBUG
238  if(matB.n!=n){
239  ERROR_REPORT;
240  std::cerr << "The matrix B is not a matrix having the same size as \"this\" matrix." << std::endl
241  << "The B matrix is (" << matB.n << "x" << matB.n << ")." << std::endl;
242  exit(1);
243  }
244 #endif//CPPL_DEBUG
245 
246  w.resize(n);
247  v.resize(n);
248  char JOBZ('V'), UPLO('l');
249  CPPL_INT ITYPE(1), LDA(n), LDB(n), LWORK(-1), INFO(1);
250  double *WORK(new double[1]);
251  dsygv_(&ITYPE, &JOBZ, &UPLO, &n, array, &LDA, matB.array, &LDB, &w[0], WORK, &LWORK, &INFO);
252 
253  INFO=1;
254  LWORK = CPPL_INT(WORK[0]);
255  delete [] WORK;
256  WORK = new double[LWORK];
257  dsygv_(&ITYPE, &JOBZ, &UPLO, &n, array, &LDA, matB.array, &LDB, &w[0], WORK, &LWORK, &INFO);
258  delete [] WORK;
259 
260  //// reforming ////
261  for(int i=0; i<n; i++){
262  v[i].resize(n);
263  for(int j=0; j<n; j++){
264  v[i](j) =darray[i][j];
265  }
266  }
267 
268  if(INFO!=0){
269  WARNING_REPORT;
270  std::cerr << "Serious trouble happend. INFO = " << INFO << "." << std::endl;
271  }
272  return INFO;
273 }
dgematrix::m
CPPL_INT m
matrix row size
Definition: dgematrix.hpp:9
dgematrix::array
double * array
1D array to store matrix data
Definition: dgematrix.hpp:11
dsymatrix::dsyev
CPPL_INT dsyev(std::vector< double > &, const bool &)
Definition: dsymatrix-lapack.hpp:84
dcovector::l
CPPL_INT l
vector size
Definition: dcovector.hpp:9
dsymatrix::dsygv
CPPL_INT dsygv(dsymatrix &, std::vector< double > &)
Definition: dsymatrix-lapack.hpp:198
dgematrix::n
CPPL_INT n
matrix column size
Definition: dgematrix.hpp:10
dgematrix
Real Double-precision General Dence Matrix Class.
Definition: dgematrix.hpp:3
dcovector::array
double * array
1D array to store vector data
Definition: dcovector.hpp:11
dsymatrix
Real Double-precision Symmetric Matrix Class [l-type (UPLO=l) Strage].
Definition: dsymatrix.hpp:3
dsymatrix::dsysv
CPPL_INT dsysv(dgematrix &)
Definition: dsymatrix-lapack.hpp:6
dsymatrix::i
friend _dsymatrix i(const dsymatrix &)
Definition: dsymatrix-calc.hpp:16
dsymatrix::n
CPPL_INT n
matrix column size
Definition: dsymatrix.hpp:10
dcovector
Real Double-precision Column Vector Class.
Definition: dcovector.hpp:3
dsymatrix::array
double * array
1D array to store matrix data
Definition: dsymatrix.hpp:11
dsymatrix::darray
double ** darray
array of pointers of column head addresses
Definition: dsymatrix.hpp:12
Generated by   doxygen 1.8.6