Warning: this documentation for the development version is under construction.
/home/vivien/public_html/.src_seldon/computation/interfaces/Lapack_Eigenvalues.cxx
00001 // Copyright (C) 2003-2009 Marc Duruflé 00002 // Copyright (C) 2001-2009 Vivien Mallet 00003 // 00004 // This file is part of the linear-algebra library Seldon, 00005 // http://seldon.sourceforge.net/. 00006 // 00007 // Seldon is free software; you can redistribute it and/or modify it under the 00008 // terms of the GNU Lesser General Public License as published by the Free 00009 // Software Foundation; either version 2.1 of the License, or (at your option) 00010 // any later version. 00011 // 00012 // Seldon is distributed in the hope that it will be useful, but WITHOUT ANY 00013 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS 00014 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for 00015 // more details. 00016 // 00017 // You should have received a copy of the GNU Lesser General Public License 00018 // along with Seldon. If not, see http://www.gnu.org/licenses/. 00019 00020 00021 #ifndef SELDON_FILE_LAPACK_EIGENVALUES_CXX 00022 00023 00024 #include "Lapack_Eigenvalues.hxx" 00025 00026 00027 namespace Seldon 00028 { 00029 00030 00032 // STANDARD EIGENVALUE PROBLEM // 00033 00034 00035 /* RowMajor */ 00036 00037 00038 template<class Prop, class Allocator1, class Allocator2, class Allocator3> 00039 void GetEigenvalues(Matrix<float, Prop, RowMajor, Allocator1>& A, 00040 Vector<float, VectFull, Allocator2>& wr, 00041 Vector<float, VectFull, Allocator3>& wi, 00042 LapackInfo& info) 00043 { 00044 int n = A.GetM(), lwork = 6*n; 00045 char jobvl('N'); 00046 char jobvr('N'); 00047 Vector<float> work(lwork); 00048 wr.Reallocate(n); 00049 wi.Reallocate(n); 00050 sgeev_(&jobvl, &jobvr, &n, A.GetData(), &n, wr.GetData(), wi.GetData(), 00051 A.GetData(), &n, A.GetData(), &n, work.GetData(), 00052 &lwork, &info.GetInfoRef()); 00053 00054 #ifdef SELDON_LAPACK_CHECK_INFO 00055 if (info.GetInfo() != 0) 00056 throw LapackError(info.GetInfo(), "GetEigenvalues", 00057 "Failed to find eigenvalues "); 00058 #endif 00059 00060 } 00061 00062 00063 template<class Prop, class Allocator1, class Allocator2, 00064 class Allocator3, class Allocator4> 00065 void GetEigenvaluesEigenvectors(Matrix<float, Prop, RowMajor, Allocator1>& A, 00066 Vector<float, VectFull, Allocator2>& wr, 00067 Vector<float, VectFull, Allocator3>& wi, 00068 Matrix<float, General, RowMajor, Allocator4>& zr, 00069 LapackInfo& info) 00070 { 00071 int n = A.GetM(), lwork = 6*n; 00072 char jobvl('V'); 00073 char jobvr('N'); 00074 Vector<float> work(lwork); 00075 wr.Reallocate(n); 00076 wi.Reallocate(n); 00077 zr.Reallocate(n, n); 00078 sgeev_(&jobvl, &jobvr, &n, A.GetData(), &n, wr.GetData(), wi.GetData(), 00079 zr.GetData(), &n, zr.GetData(), &n, work.GetData(), &lwork, 00080 &info.GetInfoRef()); 00081 00082 #ifdef SELDON_LAPACK_CHECK_INFO 00083 if (info.GetInfo() != 0) 00084 throw LapackError(info.GetInfo(), "GetEigenvalues", 00085 "Failed to find eigenvalues "); 00086 #endif 00087 00088 Transpose(zr); 00089 // conjugate if necessary 00090 int i = 0; 00091 while (i < n) 00092 { 00093 if (i < (n-1)) 00094 if (wr(i) == wr(i+1)) 00095 { 00096 for (int j = 0; j < n; j++) 00097 zr(j,i+1) = -zr(j,i+1); 00098 00099 i++; 00100 } 00101 00102 i++; 00103 } 00104 } 00105 00106 00107 template<class Prop, class Allocator1, class Allocator2> 00108 void GetEigenvalues(Matrix<complex<float>, Prop, RowMajor, Allocator1>& A, 00109 Vector<complex<float>, VectFull, Allocator2>& w, 00110 LapackInfo& info) 00111 { 00112 int n = A.GetM(); 00113 char jobl('N'), jobr('N'); int lwork = 3*n; 00114 Vector<complex<float> > work(lwork); 00115 Vector<float> rwork(2*n); 00116 w.Reallocate(n); 00117 cgeev_(&jobl, &jobr, &n, A.GetDataVoid(), &n, w.GetDataVoid(), 00118 A.GetDataVoid(), &n, A.GetData(), &n, work.GetDataVoid(), &lwork, 00119 rwork.GetData(), &info.GetInfoRef()); 00120 00121 #ifdef SELDON_LAPACK_CHECK_INFO 00122 if (info.GetInfo() != 0) 00123 throw LapackError(info.GetInfo(), "GetEigenvalues", 00124 "Failed to find eigenvalues "); 00125 #endif 00126 00127 } 00128 00129 00130 template<class Prop, class Allocator1, class Allocator2, class Allocator3> 00131 void GetEigenvaluesEigenvectors(Matrix<complex<float>, 00132 Prop, RowMajor, Allocator1>& A, 00133 Vector<complex<float>, VectFull, Allocator2>& w, 00134 Matrix<complex<float>, 00135 General, RowMajor, Allocator3>& z, 00136 LapackInfo& info) 00137 { 00138 int n = A.GetM(); 00139 char jobl('V'), jobr('N'); int lwork = 3*n; 00140 Vector<complex<float> > work(lwork); 00141 Vector<float> rwork(2*n); 00142 w.Reallocate(n); 00143 z.Reallocate(n, n); 00144 cgeev_(&jobl, &jobr, &n, A.GetDataVoid(), &n, w.GetDataVoid(), 00145 z.GetDataVoid(), &n, z.GetData(), &n, work.GetDataVoid(), &lwork, 00146 rwork.GetData(), &info.GetInfoRef()); 00147 00148 #ifdef SELDON_LAPACK_CHECK_INFO 00149 if (info.GetInfo() != 0) 00150 throw LapackError(info.GetInfo(), "GetEigenvalues", 00151 "Failed to find eigenvalues "); 00152 #endif 00153 00154 TransposeConj(z); 00155 } 00156 00157 00158 template<class Prop, class Allocator1, class Allocator2, class Allocator3> 00159 void GetEigenvalues(Matrix<double, Prop, RowMajor, Allocator1>& A, 00160 Vector<double, VectFull, Allocator2>& wr, 00161 Vector<double, VectFull, Allocator3>& wi, 00162 LapackInfo& info) 00163 { 00164 int n = A.GetM(), lwork = 6*n; 00165 char jobvl('N'), jobvr('N'); Vector<double> work(lwork); 00166 wr.Reallocate(n); 00167 wi.Reallocate(n); 00168 dgeev_(&jobvl, &jobvr, &n, A.GetData(), &n, wr.GetData(), wi.GetData(), 00169 A.GetData(), &n, A.GetData(), &n, work.GetData(), 00170 &lwork, &info.GetInfoRef()); 00171 00172 #ifdef SELDON_LAPACK_CHECK_INFO 00173 if (info.GetInfo() != 0) 00174 throw LapackError(info.GetInfo(), "GetEigenvalues", 00175 "Failed to find eigenvalues "); 00176 #endif 00177 00178 } 00179 00180 00181 template<class Prop, class Allocator1, class Allocator2, 00182 class Allocator3, class Allocator4> 00183 void GetEigenvaluesEigenvectors(Matrix<double, Prop, RowMajor, Allocator1>& A, 00184 Vector<double, VectFull, Allocator2>& wr, 00185 Vector<double, VectFull, Allocator3>& wi, 00186 Matrix<double, General, RowMajor, Allocator4>& zr, 00187 LapackInfo& info) 00188 { 00189 int n = A.GetM(), lwork = 6*n; 00190 char jobvl('V'), jobvr('N'); Vector<double> work(lwork); 00191 wr.Reallocate(n); 00192 wi.Reallocate(n); 00193 zr.Reallocate(n, n); 00194 dgeev_(&jobvl, &jobvr, &n, A.GetData(), &n, wr.GetData(), wi.GetData(), 00195 zr.GetData(), &n, zr.GetData(), &n, work.GetData(), &lwork, 00196 &info.GetInfoRef()); 00197 00198 #ifdef SELDON_LAPACK_CHECK_INFO 00199 if (info.GetInfo() != 0) 00200 throw LapackError(info.GetInfo(), "GetEigenvalues", 00201 "Failed to find eigenvalues "); 00202 #endif 00203 00204 Transpose(zr); 00205 // conjugate if necessary 00206 int i = 0; 00207 while (i < n) 00208 { 00209 if (i < (n-1)) 00210 if (wr(i) == wr(i+1)) 00211 { 00212 for (int j = 0; j < n; j++) 00213 zr(j,i+1) = -zr(j,i+1); 00214 00215 i++; 00216 } 00217 00218 i++; 00219 } 00220 } 00221 00222 00223 template<class Prop, class Allocator1, class Allocator2> 00224 void GetEigenvalues(Matrix<complex<double>, Prop, RowMajor, Allocator1>& A, 00225 Vector<complex<double>, VectFull, Allocator2>& w, 00226 LapackInfo& info) 00227 { 00228 int n = A.GetM(); 00229 char jobl('N'), jobr('N'); int lwork = 3*n; 00230 Vector<complex<double> > work(lwork); 00231 Vector<double> rwork(2*n); 00232 w.Reallocate(n); 00233 zgeev_(&jobl, &jobr, &n, A.GetDataVoid(), &n, w.GetDataVoid(), 00234 A.GetDataVoid(), &n, A.GetData(), &n, work.GetDataVoid(), &lwork, 00235 rwork.GetData(), &info.GetInfoRef()); 00236 00237 #ifdef SELDON_LAPACK_CHECK_INFO 00238 if (info.GetInfo() != 0) 00239 throw LapackError(info.GetInfo(), "GetEigenvalues", 00240 "Failed to find eigenvalues "); 00241 #endif 00242 00243 } 00244 00245 00246 template<class Prop, class Allocator1, class Allocator2, class Allocator3> 00247 void GetEigenvaluesEigenvectors(Matrix<complex<double>, 00248 Prop, RowMajor, Allocator1>& A, 00249 Vector<complex<double>, 00250 VectFull, Allocator2>& w, 00251 Matrix<complex<double>, 00252 General, RowMajor, Allocator3>& z, 00253 LapackInfo& info) 00254 { 00255 int n = A.GetM(); 00256 char jobl('V'), jobr('N'); int lwork = 3*n; 00257 Vector<complex<double> > work(lwork); 00258 Vector<double> rwork(2*n); 00259 w.Reallocate(n); 00260 z.Reallocate(n, n); 00261 zgeev_(&jobl, &jobr, &n, A.GetDataVoid(), &n, w.GetDataVoid(), 00262 z.GetDataVoid(), &n, z.GetData(), &n, work.GetDataVoid(), 00263 &lwork, rwork.GetData(), &info.GetInfoRef()); 00264 00265 #ifdef SELDON_LAPACK_CHECK_INFO 00266 if (info.GetInfo() != 0) 00267 throw LapackError(info.GetInfo(), "GetEigenvalues", 00268 "Failed to find eigenvalues "); 00269 #endif 00270 00271 TransposeConj(z); 00272 } 00273 00274 00275 /* ColMajor */ 00276 00277 00278 template<class Prop, class Allocator1, class Allocator2, class Allocator3> 00279 void GetEigenvalues(Matrix<float, Prop, ColMajor, Allocator1>& A, 00280 Vector<float, VectFull, Allocator2>& wr, 00281 Vector<float, VectFull, Allocator3>& wi, 00282 LapackInfo& info) 00283 { 00284 int n = A.GetM(), lwork = 6*n; 00285 char jobvl('N'), jobvr('N'); Vector<float> work(lwork); 00286 wr.Reallocate(n); 00287 wi.Reallocate(n); 00288 sgeev_(&jobvl, &jobvr, &n, A.GetData(), &n, wr.GetData(), wi.GetData(), 00289 A.GetData(), &n, A.GetData(), &n, work.GetData(), &lwork, 00290 &info.GetInfoRef()); 00291 00292 #ifdef SELDON_LAPACK_CHECK_INFO 00293 if (info.GetInfo() != 0) 00294 throw LapackError(info.GetInfo(), "GetEigenvalues", 00295 "Failed to find eigenvalues "); 00296 #endif 00297 00298 } 00299 00300 template<class Prop, class Allocator1, class Allocator2, 00301 class Allocator3, class Allocator4> 00302 void GetEigenvaluesEigenvectors(Matrix<float, Prop, ColMajor, Allocator1>& A, 00303 Vector<float, VectFull, Allocator2>& wr, 00304 Vector<float, VectFull, Allocator3>& wi, 00305 Matrix<float, General, ColMajor, Allocator4>&zr, 00306 LapackInfo& info) 00307 { 00308 int n = A.GetM(), lwork = 6*n; 00309 char jobvl('N'), jobvr('V'); 00310 Vector<float> work(lwork); 00311 wr.Reallocate(n); 00312 wi.Reallocate(n); 00313 zr.Reallocate(n, n); 00314 sgeev_(&jobvl, &jobvr, &n, A.GetData(), &n, wr.GetData(), wi.GetData(), 00315 zr.GetData(), &n, zr.GetData(), &n, work.GetData(), 00316 &lwork, &info.GetInfoRef()); 00317 00318 #ifdef SELDON_LAPACK_CHECK_INFO 00319 if (info.GetInfo() != 0) 00320 throw LapackError(info.GetInfo(), "GetEigenvalues", 00321 "Failed to find eigenvalues "); 00322 #endif 00323 00324 } 00325 00326 00327 template<class Prop, class Allocator1, class Allocator2> 00328 void GetEigenvalues(Matrix<complex<float>, Prop, ColMajor, Allocator1>& A, 00329 Vector<complex<float>, VectFull, Allocator2>& w, 00330 LapackInfo& info) 00331 { 00332 int n = A.GetM(); 00333 char jobl('N'), jobr('N'); int lwork = 3*n; 00334 Vector<complex<float> > work(lwork); 00335 Vector<float> rwork(2*n); 00336 w.Reallocate(n); 00337 cgeev_(&jobl, &jobr, &n, A.GetDataVoid(), &n, w.GetDataVoid(), 00338 A.GetDataVoid(), &n, A.GetData(), &n, work.GetDataVoid(), 00339 &lwork, rwork.GetData(), &info.GetInfoRef()); 00340 00341 #ifdef SELDON_LAPACK_CHECK_INFO 00342 if (info.GetInfo() != 0) 00343 throw LapackError(info.GetInfo(), "GetEigenvalues", 00344 "Failed to find eigenvalues "); 00345 #endif 00346 00347 } 00348 00349 00350 template<class Prop, class Allocator1, class Allocator2, class Allocator3> 00351 void GetEigenvaluesEigenvectors(Matrix<complex<float>, 00352 Prop, ColMajor, Allocator1>& A, 00353 Vector<complex<float>, 00354 VectFull, Allocator2>& w, 00355 Matrix<complex<float>, 00356 General, ColMajor, Allocator3>& z, 00357 LapackInfo& info) 00358 { 00359 int n = A.GetM(); 00360 char jobl('N'), jobr('V'); int lwork = 3*n; 00361 Vector<complex<float> > work(lwork); 00362 Vector<float> rwork(2*n); 00363 w.Reallocate(n); 00364 z.Reallocate(n, n); 00365 cgeev_(&jobl, &jobr, &n, A.GetDataVoid(), &n, w.GetDataVoid(), 00366 z.GetDataVoid(), &n, z.GetData(), &n, work.GetDataVoid(), 00367 &lwork, rwork.GetData(), &info.GetInfoRef()); 00368 00369 #ifdef SELDON_LAPACK_CHECK_INFO 00370 if (info.GetInfo() != 0) 00371 throw LapackError(info.GetInfo(), "GetEigenvalues", 00372 "Failed to find eigenvalues "); 00373 #endif 00374 00375 } 00376 00377 template<class Prop, class Allocator1, class Allocator2, class Allocator3> 00378 void GetEigenvalues(Matrix<double, Prop, ColMajor, Allocator1>& A, 00379 Vector<double, VectFull, Allocator2>& wr, 00380 Vector<double, VectFull, Allocator3>& wi, 00381 LapackInfo& info) 00382 { 00383 int n = A.GetM(), lwork = 6*n; 00384 char jobvl('N'), jobvr('N'); Vector<double> work(lwork); 00385 wr.Reallocate(n); 00386 wi.Reallocate(n); 00387 dgeev_(&jobvl, &jobvr, &n, A.GetData(), &n, wr.GetData(), wi.GetData(), 00388 A.GetData(), &n, A.GetData(), &n, work.GetData(), &lwork, 00389 &info.GetInfoRef()); 00390 00391 #ifdef SELDON_LAPACK_CHECK_INFO 00392 if (info.GetInfo() != 0) 00393 throw LapackError(info.GetInfo(), "GetEigenvalues", 00394 "Failed to find eigenvalues "); 00395 #endif 00396 00397 } 00398 00399 00400 template<class Prop, class Allocator1, class Allocator2, 00401 class Allocator3, class Allocator4> 00402 void GetEigenvaluesEigenvectors(Matrix<double, Prop, ColMajor, Allocator1>& A, 00403 Vector<double, VectFull, Allocator2>& wr, 00404 Vector<double, VectFull, Allocator3>& wi, 00405 Matrix<double, General, ColMajor, Allocator4>&zr, 00406 LapackInfo& info) 00407 { 00408 int n = A.GetM(), lwork = 6*n; 00409 char jobvl('N'), jobvr('V'); 00410 Vector<double> work(lwork); 00411 wr.Reallocate(n); 00412 wi.Reallocate(n); 00413 zr.Reallocate(n, n); 00414 dgeev_(&jobvl, &jobvr, &n, A.GetData(), &n, wr.GetData(), wi.GetData(), 00415 zr.GetData(), &n, zr.GetData(), &n, work.GetData(), 00416 &lwork, &info.GetInfoRef()); 00417 00418 #ifdef SELDON_LAPACK_CHECK_INFO 00419 if (info.GetInfo() != 0) 00420 throw LapackError(info.GetInfo(), "GetEigenvalues", 00421 "Failed to find eigenvalues "); 00422 #endif 00423 00424 } 00425 00426 00427 template<class Prop, class Allocator1, class Allocator2> 00428 void GetEigenvalues(Matrix<complex<double>, Prop, ColMajor, Allocator1>& A, 00429 Vector<complex<double>, VectFull, Allocator2>& w, 00430 LapackInfo& info) 00431 { 00432 int n = A.GetM(); 00433 char jobl('N'), jobr('N'); int lwork = 3*n; 00434 Vector<complex<double> > work(lwork); 00435 Vector<double> rwork(2*n); 00436 w.Reallocate(n); 00437 zgeev_(&jobl, &jobr, &n, A.GetDataVoid(), &n, w.GetDataVoid(), 00438 A.GetDataVoid(), &n, A.GetData(), &n, work.GetDataVoid(), 00439 &lwork, rwork.GetData(), &info.GetInfoRef()); 00440 00441 #ifdef SELDON_LAPACK_CHECK_INFO 00442 if (info.GetInfo() != 0) 00443 throw LapackError(info.GetInfo(), "GetEigenvalues", 00444 "Failed to find eigenvalues "); 00445 #endif 00446 00447 } 00448 00449 00450 template<class Prop, class Allocator1, class Allocator2, class Allocator3> 00451 void GetEigenvaluesEigenvectors(Matrix<complex<double>, 00452 Prop, ColMajor, Allocator1>& A, 00453 Vector<complex<double>, 00454 VectFull, Allocator2>& w, 00455 Matrix<complex<double>, 00456 General, ColMajor, Allocator3>& z, 00457 LapackInfo& info) 00458 { 00459 int n = A.GetM(); 00460 char jobl('N'), jobr('V'); int lwork = 3*n; 00461 Vector<complex<double> > work(lwork); 00462 Vector<double> rwork(2*n); 00463 w.Reallocate(n); 00464 z.Reallocate(n, n); 00465 zgeev_(&jobl, &jobr, &n, A.GetDataVoid(), &n, w.GetDataVoid(), 00466 z.GetDataVoid(), &n, z.GetData(), &n, work.GetDataVoid(), 00467 &lwork, rwork.GetData(), &info.GetInfoRef()); 00468 00469 #ifdef SELDON_LAPACK_CHECK_INFO 00470 if (info.GetInfo() != 0) 00471 throw LapackError(info.GetInfo(), "GetEigenvalues", 00472 "Failed to find eigenvalues "); 00473 #endif 00474 00475 } 00476 00477 00478 /* RowSym */ 00479 00480 00481 template<class Prop, class Allocator1, class Allocator2> 00482 void GetEigenvalues(Matrix<float, Prop, RowSym, Allocator1>& A, 00483 Vector<float, VectFull, Allocator2>& w, 00484 LapackInfo& info) 00485 { 00486 int n = A.GetM(); 00487 char uplo('L'); char job('N'); 00488 int lwork = 3*n; Vector<float> work(lwork); 00489 w.Reallocate(n); 00490 ssyev_(&job, &uplo, &n, A.GetData(), &n, w.GetData(), work.GetData(), 00491 &lwork, &info.GetInfoRef()); 00492 00493 #ifdef SELDON_LAPACK_CHECK_INFO 00494 if (info.GetInfo() != 0) 00495 throw LapackError(info.GetInfo(), "GetEigenvalues", 00496 "Failed to find eigenvalues "); 00497 #endif 00498 00499 } 00500 00501 00502 template<class Prop, class Allocator1, class Allocator2, class Allocator3> 00503 void GetEigenvaluesEigenvectors(Matrix<float, Prop, RowSym, Allocator1>& A, 00504 Vector<float, VectFull, Allocator2>& w, 00505 Matrix<float, General, RowMajor, Allocator3>& z, 00506 LapackInfo& info) 00507 { 00508 int n = A.GetM(); 00509 char uplo('L'); char job('V'); 00510 int lwork = 3*n; Vector<float> work(lwork); 00511 w.Reallocate(n); 00512 z.Reallocate(n, n); 00513 for (int i = 0; i < n; i++) 00514 for (int j = 0; j < n; j++) 00515 z(i,j) = A(i,j); 00516 00517 ssyev_(&job, &uplo, &n, z.GetData(), &n, w.GetData(), work.GetData(), 00518 &lwork, &info.GetInfoRef()); 00519 00520 #ifdef SELDON_LAPACK_CHECK_INFO 00521 if (info.GetInfo() != 0) 00522 throw LapackError(info.GetInfo(), "GetEigenvalues", 00523 "Failed to find eigenvalues "); 00524 #endif 00525 00526 Transpose(z); 00527 } 00528 00529 template<class Prop, class Allocator1, class Allocator2> 00530 void GetEigenvalues(Matrix<complex<float>, Prop, RowSym, Allocator1>& A, 00531 Vector<complex<float>, VectFull, Allocator2>& w, 00532 LapackInfo& info) 00533 { 00534 int n = A.GetM(); 00535 w.Reallocate(n); 00536 Matrix<complex<float>, General, ColMajor> B(n,n); 00537 for (int i = 0; i < n; i++) 00538 for (int j = 0; j < n; j++) 00539 B(i,j) = A(i,j); 00540 00541 GetEigenvalues(B, w); 00542 } 00543 00544 template<class Prop, class Allocator1, class Allocator2, class Allocator3> 00545 void GetEigenvaluesEigenvectors(Matrix<complex<float>, 00546 Prop, RowSym, Allocator1>& A, 00547 Vector<complex<float>, 00548 VectFull, Allocator2>& w, 00549 Matrix<complex<float>, 00550 General, RowMajor, Allocator3>& z, 00551 LapackInfo& info) 00552 { 00553 int n = A.GetM(); 00554 w.Reallocate(n); 00555 z.Reallocate(n, n); 00556 Matrix<complex<float>, General, RowMajor> B(n,n); 00557 for (int i = 0; i < n; i++) 00558 for (int j = 0; j < n; j++) 00559 B(i,j) = A(i,j); 00560 00561 GetEigenvaluesEigenvectors(B, w, z); 00562 } 00563 00564 00565 template<class Prop, class Allocator1, class Allocator2> 00566 void GetEigenvalues(Matrix<double, Prop, RowSym, Allocator1>& A, 00567 Vector<double, VectFull, Allocator2>& w, 00568 LapackInfo& info) 00569 { 00570 int n = A.GetM(); 00571 char uplo('L'); char job('N'); 00572 int lwork = 3*n; Vector<double> work(lwork); 00573 w.Reallocate(n); 00574 dsyev_(&job, &uplo, &n, A.GetData(), &n, w.GetData(), work.GetData(), 00575 &lwork, &info.GetInfoRef()); 00576 00577 #ifdef SELDON_LAPACK_CHECK_INFO 00578 if (info.GetInfo() != 0) 00579 throw LapackError(info.GetInfo(), "GetEigenvalues", 00580 "Failed to find eigenvalues "); 00581 #endif 00582 00583 } 00584 00585 00586 template<class Prop, class Allocator1, class Allocator2, class Allocator3> 00587 void GetEigenvaluesEigenvectors(Matrix<double, Prop, RowSym, Allocator1>& A, 00588 Vector<double, VectFull, Allocator2>& w, 00589 Matrix<double, General, RowMajor, Allocator3>& z, 00590 LapackInfo& info) 00591 { 00592 int n = A.GetM(); 00593 char uplo('L'); char job('V'); 00594 int lwork = 3*n; Vector<double> work(lwork); 00595 w.Reallocate(n); 00596 z.Reallocate(n, n); 00597 for (int i = 0; i < n; i++) 00598 for (int j = 0; j < n; j++) 00599 z(i,j) = A(i,j); 00600 00601 dsyev_(&job, &uplo, &n, z.GetData(), &n, w.GetData(), work.GetData(), 00602 &lwork, &info.GetInfoRef()); 00603 00604 #ifdef SELDON_LAPACK_CHECK_INFO 00605 if (info.GetInfo() != 0) 00606 throw LapackError(info.GetInfo(), "GetEigenvalues", 00607 "Failed to find eigenvalues "); 00608 #endif 00609 00610 Transpose(z); 00611 } 00612 00613 template<class Prop, class Allocator1, class Allocator2> 00614 void GetEigenvalues(Matrix<complex<double>, Prop, RowSym, Allocator1>& A, 00615 Vector<complex<double>, VectFull, Allocator2>& w, 00616 LapackInfo& info) 00617 { 00618 int n = A.GetM(); 00619 w.Reallocate(n); 00620 Matrix<complex<double>, General, ColMajor> B(n,n); 00621 for (int i = 0; i < n; i++) 00622 for (int j = 0; j < n; j++) 00623 B(i,j) = A(i,j); 00624 00625 GetEigenvalues(B, w); 00626 } 00627 00628 template<class Prop, class Allocator1, class Allocator2, class Allocator3> 00629 void GetEigenvaluesEigenvectors(Matrix<complex<double>, 00630 Prop, RowSym, Allocator1>& A, 00631 Vector<complex<double>, 00632 VectFull, Allocator2>& w, 00633 Matrix<complex<double>, 00634 General, RowMajor, Allocator3>& z, 00635 LapackInfo& info) 00636 { 00637 int n = A.GetM(); 00638 Matrix<complex<double>, General, RowMajor> B(n,n); 00639 for (int i = 0; i < n; i++) 00640 for (int j = 0; j < n; j++) 00641 B(i,j) = A(i,j); 00642 00643 w.Reallocate(n); 00644 z.Reallocate(n, n); 00645 GetEigenvaluesEigenvectors(B, w, z); 00646 } 00647 00648 00649 /* ColSym */ 00650 00651 00652 template<class Prop, class Allocator1, class Allocator2> 00653 void GetEigenvalues(Matrix<float, Prop, ColSym, Allocator1>& A, 00654 Vector<float, VectFull, Allocator2>& w, 00655 LapackInfo& info) 00656 { 00657 int n = A.GetM(); 00658 char uplo('U'); char job('N'); int lwork = 3*n; Vector<float> work(lwork); 00659 w.Reallocate(n); 00660 ssyev_(&job, &uplo, &n, A.GetData(), &n, w.GetData(), work.GetData(), 00661 &lwork, &info.GetInfoRef()); 00662 00663 #ifdef SELDON_LAPACK_CHECK_INFO 00664 if (info.GetInfo() != 0) 00665 throw LapackError(info.GetInfo(), "GetEigenvalues", 00666 "Failed to find eigenvalues "); 00667 #endif 00668 00669 } 00670 00671 00672 template<class Prop, class Allocator1, class Allocator2, class Allocator3> 00673 void GetEigenvaluesEigenvectors(Matrix<float, Prop, ColSym, Allocator1>& A, 00674 Vector<float, VectFull, Allocator2>& w, 00675 Matrix<float, General, ColMajor, Allocator3>& z, 00676 LapackInfo& info) 00677 { 00678 int n = A.GetM(); 00679 char uplo('U'); char job('V'); 00680 int lwork = 3*n; Vector<float> work(lwork); 00681 w.Reallocate(n); 00682 z.Reallocate(n, n); 00683 for (int i = 0; i < n; i++) 00684 for (int j = 0; j < n; j++) 00685 z(i,j) = A(i,j); 00686 00687 ssyev_(&job, &uplo, &n, z.GetData(), &n, w.GetData(), 00688 work.GetData(), &lwork, &info.GetInfoRef()); 00689 00690 #ifdef SELDON_LAPACK_CHECK_INFO 00691 if (info.GetInfo() != 0) 00692 throw LapackError(info.GetInfo(), "GetEigenvalues", 00693 "Failed to find eigenvalues "); 00694 #endif 00695 00696 } 00697 00698 00699 template<class Prop, class Allocator1, class Allocator2> 00700 void GetEigenvalues(Matrix<complex<float>, Prop, ColSym, Allocator1>& A, 00701 Vector<complex<float>, VectFull, Allocator2>& w, 00702 LapackInfo& info) 00703 { 00704 int n = A.GetM(); 00705 w.Reallocate(n); 00706 Matrix<complex<float>, General, ColMajor> B(n,n); 00707 for (int i = 0; i < n; i++) 00708 for (int j = 0; j < n; j++) 00709 B(i,j) = A(i,j); 00710 00711 GetEigenvalues(B, w); 00712 } 00713 00714 00715 template<class Prop, class Allocator1, class Allocator2, class Allocator3> 00716 void GetEigenvaluesEigenvectors(Matrix<complex<float>, 00717 Prop, ColSym, Allocator1>& A, 00718 Vector<complex<float>, 00719 VectFull, Allocator2>& w, 00720 Matrix<complex<float>, 00721 General, ColMajor, Allocator3>& z, 00722 LapackInfo& info) 00723 { 00724 int n = A.GetM(); 00725 Matrix<complex<float>, General, ColMajor> B(n,n); 00726 for (int i = 0; i < n; i++) 00727 for (int j = 0; j < n; j++) 00728 B(i,j) = A(i,j); 00729 00730 w.Reallocate(n); 00731 z.Reallocate(n, n); 00732 GetEigenvaluesEigenvectors(B, w, z); 00733 } 00734 00735 00736 template<class Prop, class Allocator1, class Allocator2> 00737 void GetEigenvalues(Matrix<double, Prop, ColSym, Allocator1>& A, 00738 Vector<double, VectFull, Allocator2>& w, 00739 LapackInfo& info) 00740 { 00741 int n = A.GetM(); 00742 char uplo('U'); char job('N'); int lwork = 3*n; Vector<double> work(lwork); 00743 w.Reallocate(n); 00744 dsyev_(&job, &uplo, &n, A.GetData(), &n, w.GetData(), work.GetData(), 00745 &lwork, &info.GetInfoRef()); 00746 00747 #ifdef SELDON_LAPACK_CHECK_INFO 00748 if (info.GetInfo() != 0) 00749 throw LapackError(info.GetInfo(), "GetEigenvalues", 00750 "Failed to find eigenvalues "); 00751 #endif 00752 00753 } 00754 00755 00756 template<class Prop, class Allocator1, class Allocator2, class Allocator3> 00757 void GetEigenvaluesEigenvectors(Matrix<double, Prop, ColSym, Allocator1>& A, 00758 Vector<double, VectFull, Allocator2>& w, 00759 Matrix<double, General, ColMajor, Allocator3>& z, 00760 LapackInfo& info) 00761 { 00762 int n = A.GetM(); 00763 char uplo('U'); char job('V'); 00764 int lwork = 3*n; Vector<double> work(lwork); 00765 w.Reallocate(n); 00766 z.Reallocate(n, n); 00767 for (int i = 0; i < n; i++) 00768 for (int j = 0; j < n; j++) 00769 z(i,j) = A(i,j); 00770 00771 dsyev_(&job, &uplo, &n, z.GetData(), &n, w.GetData(), 00772 work.GetData(), &lwork, &info.GetInfoRef()); 00773 00774 #ifdef SELDON_LAPACK_CHECK_INFO 00775 if (info.GetInfo() != 0) 00776 throw LapackError(info.GetInfo(), "GetEigenvalues", 00777 "Failed to find eigenvalues "); 00778 #endif 00779 00780 } 00781 00782 00783 template<class Prop, class Allocator1, class Allocator2> 00784 void GetEigenvalues(Matrix<complex<double>, Prop, ColSym, Allocator1>& A, 00785 Vector<complex<double>, VectFull, Allocator2>& w, 00786 LapackInfo& info) 00787 { 00788 int n = A.GetM(); 00789 Matrix<complex<double>, General, ColMajor> B(n,n); 00790 w.Reallocate(n); 00791 for (int i = 0; i < n; i++) 00792 for (int j = 0; j < n; j++) 00793 B(i,j) = A(i,j); 00794 00795 GetEigenvalues(B, w); 00796 } 00797 00798 00799 template<class Prop, class Allocator1, class Allocator2, class Allocator3> 00800 void GetEigenvaluesEigenvectors(Matrix<complex<double>, 00801 Prop, ColSym, Allocator1>& A, 00802 Vector<complex<double>, 00803 VectFull, Allocator2>& w, 00804 Matrix<complex<double>, 00805 General, ColMajor, Allocator3>& z, 00806 LapackInfo& info) 00807 { 00808 int n = A.GetM(); 00809 Matrix<complex<double>, General, ColMajor> B(n,n); 00810 for (int i = 0; i < n; i++) 00811 for (int j = 0; j < n; j++) 00812 B(i,j) = A(i,j); 00813 00814 w.Reallocate(n); 00815 z.Reallocate(n, n); 00816 GetEigenvaluesEigenvectors(B, w, z); 00817 } 00818 00819 00820 /* RowHerm */ 00821 00822 00823 template<class Prop, class Allocator1, class Allocator2> 00824 void GetEigenvalues(Matrix<complex<float>, Prop, RowHerm, Allocator1>& A, 00825 Vector<float, VectFull, Allocator2>& w, 00826 LapackInfo& info) 00827 { 00828 int n = A.GetM(); 00829 char uplo('L'); char job('N'); 00830 int lwork = 2*n; Vector<complex<float> > work(lwork); 00831 Vector<float> rwork(3*n); 00832 w.Reallocate(n); 00833 cheev_(&job, &uplo, &n, A.GetDataVoid(), &n, w.GetData(), 00834 work.GetDataVoid(), &lwork, rwork.GetData(), 00835 &info.GetInfoRef()); 00836 00837 #ifdef SELDON_LAPACK_CHECK_INFO 00838 if (info.GetInfo() != 0) 00839 throw LapackError(info.GetInfo(), "GetEigenvalues", 00840 "Failed to find eigenvalues "); 00841 #endif 00842 00843 } 00844 00845 00846 template<class Prop, class Allocator1, class Allocator2, class Allocator3> 00847 void GetEigenvaluesEigenvectors(Matrix<complex<float>, 00848 Prop, RowHerm, Allocator1>& A, 00849 Vector<float, VectFull, Allocator2>& w, 00850 Matrix<complex<float>, 00851 General, RowMajor, Allocator3>& z, 00852 LapackInfo& info) 00853 { 00854 int n = A.GetM(); 00855 char uplo('L'); char job('V'); 00856 int lwork = 2*n; Vector<complex<float> > work(lwork); 00857 Vector<float> rwork(3*n); 00858 w.Reallocate(n); 00859 z.Reallocate(n,n); 00860 for (int i = 0; i < n; i++) 00861 for (int j = 0; j < n; j++) 00862 z(i,j) = A(i,j); 00863 00864 cheev_(&job, &uplo,&n, z.GetDataVoid(),&n, w.GetData(), work.GetDataVoid(), 00865 &lwork, rwork.GetData(), &info.GetInfoRef()); 00866 00867 #ifdef SELDON_LAPACK_CHECK_INFO 00868 if (info.GetInfo() != 0) 00869 throw LapackError(info.GetInfo(), "GetEigenvalues", 00870 "Failed to find eigenvalues "); 00871 #endif 00872 00873 Transpose(z); 00874 } 00875 00876 template<class Prop, class Allocator1, class Allocator2> 00877 void GetEigenvalues(Matrix<complex<double>, Prop, RowHerm, Allocator1>& A, 00878 Vector<double, VectFull, Allocator2>& w, 00879 LapackInfo& info) 00880 { 00881 int n = A.GetM(); 00882 char uplo('L'); char job('N'); 00883 int lwork = 2*n; Vector<complex<double> > work(lwork); 00884 Vector<double> rwork(3*n); 00885 w.Reallocate(n); 00886 zheev_(&job, &uplo, &n, A.GetDataVoid(), &n, w.GetData(), 00887 work.GetDataVoid(), &lwork, rwork.GetData(), 00888 &info.GetInfoRef()); 00889 00890 #ifdef SELDON_LAPACK_CHECK_INFO 00891 if (info.GetInfo() != 0) 00892 throw LapackError(info.GetInfo(), "GetEigenvalues", 00893 "Failed to find eigenvalues "); 00894 #endif 00895 00896 } 00897 00898 00899 template<class Prop, class Allocator1, class Allocator2, class Allocator3> 00900 void GetEigenvaluesEigenvectors(Matrix<complex<double>, 00901 Prop, RowHerm, Allocator1>& A, 00902 Vector<double, VectFull, Allocator2>& w, 00903 Matrix<complex<double>, 00904 General, RowMajor, Allocator3>& z, 00905 LapackInfo& info) 00906 { 00907 int n = A.GetM(); 00908 char uplo('L'); char job('V'); 00909 int lwork = 2*n; Vector<complex<double> > work(lwork); 00910 Vector<double> rwork(3*n); 00911 w.Reallocate(n); 00912 z.Reallocate(n,n); 00913 for (int i = 0; i < n; i++) 00914 for (int j = 0; j < n; j++) 00915 z(i,j) = A(i,j); 00916 00917 zheev_(&job, &uplo,&n, z.GetDataVoid(),&n, w.GetData(), work.GetDataVoid(), 00918 &lwork, rwork.GetData() , &info.GetInfoRef()); 00919 00920 #ifdef SELDON_LAPACK_CHECK_INFO 00921 if (info.GetInfo() != 0) 00922 throw LapackError(info.GetInfo(), "GetEigenvalues", 00923 "Failed to find eigenvalues "); 00924 #endif 00925 00926 Transpose(z); 00927 } 00928 00929 00930 /* ColHerm */ 00931 00932 00933 template<class Prop, class Allocator1, class Allocator2> 00934 void GetEigenvalues(Matrix<complex<float>, Prop, ColHerm, Allocator1>& A, 00935 Vector<float, VectFull, Allocator2>& w, 00936 LapackInfo& info) 00937 { 00938 int n = A.GetM(); 00939 char uplo('U'); char job('N'); int lwork = 2*n; 00940 Vector<complex<float> > work(lwork); 00941 Vector<float> rwork(3*n); 00942 w.Reallocate(n); 00943 cheev_(&job, &uplo, &n, A.GetDataVoid(), &n, w.GetData(), 00944 work.GetDataVoid(), &lwork, rwork.GetData(), 00945 &info.GetInfoRef()); 00946 00947 #ifdef SELDON_LAPACK_CHECK_INFO 00948 if (info.GetInfo() != 0) 00949 throw LapackError(info.GetInfo(), "GetEigenvalues", 00950 "Failed to find eigenvalues "); 00951 #endif 00952 00953 } 00954 00955 00956 template<class Prop, class Allocator1, class Allocator2, class Allocator3> 00957 void GetEigenvaluesEigenvectors(Matrix<complex<float>, 00958 Prop, ColHerm, Allocator1>& A, 00959 Vector<float, VectFull, Allocator2>& w, 00960 Matrix<complex<float>, 00961 General, ColMajor, Allocator3>& z, 00962 LapackInfo& info) 00963 { 00964 int n = A.GetM(); 00965 char uplo('U'); char job('V'); int lwork = 2*n; 00966 Vector<complex<float> > work(lwork); 00967 Vector<float> rwork(3*n); 00968 w.Reallocate(n); 00969 z.Reallocate(n,n); 00970 for (int i = 0; i < n; i++) 00971 for (int j = 0; j < n; j++) 00972 z(i,j) = A(i,j); 00973 00974 cheev_(&job, &uplo, &n, z.GetDataVoid(), &n, 00975 w.GetData(), work.GetDataVoid(), 00976 &lwork, rwork.GetData() , &info.GetInfoRef()); 00977 00978 #ifdef SELDON_LAPACK_CHECK_INFO 00979 if (info.GetInfo() != 0) 00980 throw LapackError(info.GetInfo(), "GetEigenvalues", 00981 "Failed to find eigenvalues "); 00982 #endif 00983 00984 } 00985 00986 00987 template<class Prop, class Allocator1, class Allocator2> 00988 void GetEigenvalues(Matrix<complex<double>, Prop, ColHerm, Allocator1>& A, 00989 Vector<double, VectFull, Allocator2>& w, 00990 LapackInfo& info) 00991 { 00992 int n = A.GetM(); 00993 char uplo('U'); char job('N'); int lwork = 2*n; 00994 Vector<complex<double> > work(lwork); 00995 Vector<double> rwork(3*n); 00996 w.Reallocate(n); 00997 zheev_(&job, &uplo, &n, A.GetDataVoid(), &n, w.GetData(), 00998 work.GetDataVoid(), &lwork, rwork.GetData(), 00999 &info.GetInfoRef()); 01000 01001 #ifdef SELDON_LAPACK_CHECK_INFO 01002 if (info.GetInfo() != 0) 01003 throw LapackError(info.GetInfo(), "GetEigenvalues", 01004 "Failed to find eigenvalues "); 01005 #endif 01006 01007 } 01008 01009 01010 template<class Prop, class Allocator1, class Allocator2, class Allocator3> 01011 void GetEigenvaluesEigenvectors(Matrix<complex<double>, 01012 Prop, ColHerm, Allocator1>& A, 01013 Vector<double, VectFull, Allocator2>& w, 01014 Matrix<complex<double>, 01015 General, ColMajor, Allocator3>& z, 01016 LapackInfo& info) 01017 { 01018 int n = A.GetM(); 01019 char uplo('U'); char job('V'); int lwork = 2*n; 01020 Vector<complex<double> > work(lwork); 01021 Vector<double> rwork(3*n); 01022 w.Reallocate(n); 01023 z.Reallocate(n,n); 01024 for (int i = 0; i < n; i++) 01025 for (int j = 0; j < n; j++) 01026 z(i,j) = A(i,j); 01027 01028 zheev_(&job, &uplo, &n, z.GetDataVoid(), &n, 01029 w.GetData(), work.GetDataVoid(), 01030 &lwork, rwork.GetData() , &info.GetInfoRef()); 01031 01032 #ifdef SELDON_LAPACK_CHECK_INFO 01033 if (info.GetInfo() != 0) 01034 throw LapackError(info.GetInfo(), "GetEigenvalues", 01035 "Failed to find eigenvalues "); 01036 #endif 01037 01038 } 01039 01040 01041 /* RowSymPacked */ 01042 01043 01044 template<class Prop, class Allocator1, class Allocator2> 01045 void GetEigenvalues(Matrix<float, Prop, RowSymPacked, Allocator1>& A, 01046 Vector<float, VectFull, Allocator2>& w, 01047 LapackInfo& info) 01048 { 01049 int n = A.GetM(); 01050 char uplo('L'); char job('N'); 01051 Vector<float> work(3*n); 01052 w.Reallocate(n); 01053 sspev_(&job, &uplo, &n, A.GetData(), w.GetData(), A.GetData(), &n, 01054 work.GetData() , &info.GetInfoRef()); 01055 01056 #ifdef SELDON_LAPACK_CHECK_INFO 01057 if (info.GetInfo() != 0) 01058 throw LapackError(info.GetInfo(), "GetEigenvalues", 01059 "Failed to find eigenvalues "); 01060 #endif 01061 01062 } 01063 01064 01065 template<class Prop, class Allocator1, class Allocator2, class Allocator3> 01066 void GetEigenvaluesEigenvectors(Matrix<float, Prop,RowSymPacked, Allocator1>& A, 01067 Vector<float, VectFull, Allocator2>& w, 01068 Matrix<float, General, RowMajor, Allocator3>&z, 01069 LapackInfo& info) 01070 { 01071 int n = A.GetM(); 01072 char uplo('L'); char job('V'); 01073 Vector<float> work(3*n); 01074 w.Reallocate(n); 01075 z.Reallocate(n,n); 01076 sspev_(&job, &uplo, &n, A.GetData(), w.GetData(), z.GetData(), &n, 01077 work.GetData() , &info.GetInfoRef()); 01078 01079 #ifdef SELDON_LAPACK_CHECK_INFO 01080 if (info.GetInfo() != 0) 01081 throw LapackError(info.GetInfo(), "GetEigenvalues", 01082 "Failed to find eigenvalues "); 01083 #endif 01084 01085 Transpose(z); 01086 } 01087 01088 01089 template<class Prop, class Allocator1, class Allocator2> 01090 void GetEigenvalues(Matrix<complex<float>, Prop, RowSymPacked, Allocator1>& A, 01091 Vector<complex<float>, VectFull, Allocator2>& w, 01092 LapackInfo& info) 01093 { 01094 int n = A.GetM(); 01095 Matrix<complex<float>, General, ColMajor> B(n,n); 01096 for (int i = 0; i < n; i++) 01097 for (int j = 0; j < n; j++) 01098 B(i,j) = A(i,j); 01099 01100 w.Reallocate(n); 01101 GetEigenvalues(B, w); 01102 } 01103 01104 template<class Prop, class Allocator1, class Allocator2, class Allocator3> 01105 void GetEigenvaluesEigenvectors(Matrix<complex<float>, Prop, 01106 RowSymPacked, Allocator1>& A, 01107 Vector<complex<float>,VectFull, Allocator2>& w, 01108 Matrix<complex<float>, 01109 General, RowMajor, Allocator3>& z, 01110 LapackInfo& info) 01111 { 01112 int n = A.GetM(); 01113 Matrix<complex<float>, General, RowMajor> B(n,n); 01114 for (int i = 0; i < n; i++) 01115 for (int j = 0; j < n; j++) 01116 B(i,j) = A(i,j); 01117 01118 w.Reallocate(n); 01119 z.Reallocate(n,n); 01120 GetEigenvaluesEigenvectors(B, w, z); 01121 } 01122 01123 01124 template<class Prop, class Allocator1, class Allocator2> 01125 void GetEigenvalues(Matrix<double, Prop, RowSymPacked, Allocator1>& A, 01126 Vector<double, VectFull, Allocator2>& w, 01127 LapackInfo& info) 01128 { 01129 int n = A.GetM(); 01130 char uplo('L'); char job('N'); Vector<double> work(3*n); 01131 w.Reallocate(n); 01132 dspev_(&job, &uplo, &n, A.GetData(), w.GetData(), A.GetData(), &n, 01133 work.GetData() , &info.GetInfoRef()); 01134 01135 #ifdef SELDON_LAPACK_CHECK_INFO 01136 if (info.GetInfo() != 0) 01137 throw LapackError(info.GetInfo(), "GetEigenvalues", 01138 "Failed to find eigenvalues "); 01139 #endif 01140 01141 } 01142 01143 01144 template<class Prop, class Allocator1, class Allocator2, class Allocator3> 01145 void GetEigenvaluesEigenvectors(Matrix<double,Prop,RowSymPacked, Allocator1>& A, 01146 Vector<double, VectFull, Allocator2>& w, 01147 Matrix<double, General, RowMajor, Allocator3>&z, 01148 LapackInfo& info) 01149 { 01150 int n = A.GetM(); 01151 char uplo('L'); char job('V'); Vector<double> work(3*n); 01152 w.Reallocate(n); 01153 z.Reallocate(n,n); 01154 dspev_(&job, &uplo, &n, A.GetData(), w.GetData(), z.GetData(), &n, 01155 work.GetData() , &info.GetInfoRef()); 01156 01157 #ifdef SELDON_LAPACK_CHECK_INFO 01158 if (info.GetInfo() != 0) 01159 throw LapackError(info.GetInfo(), "GetEigenvalues", 01160 "Failed to find eigenvalues "); 01161 #endif 01162 01163 Transpose(z); 01164 } 01165 01166 01167 template<class Prop, class Allocator1, class Allocator2> 01168 void GetEigenvalues(Matrix<complex<double>, Prop, RowSymPacked, Allocator1>& A, 01169 Vector<complex<double>, VectFull, Allocator2>& w, 01170 LapackInfo& info) 01171 { 01172 int n = A.GetM(); 01173 Matrix<complex<double>, General, ColMajor> B(n,n); 01174 for (int i = 0; i < n; i++) 01175 for (int j = 0; j < n; j++) 01176 B(i,j) = A(i,j); 01177 01178 w.Reallocate(n); 01179 GetEigenvalues(B, w); 01180 } 01181 01182 01183 template<class Prop, class Allocator1, class Allocator2, class Allocator3> 01184 void GetEigenvaluesEigenvectors(Matrix<complex<double>, Prop, 01185 RowSymPacked, Allocator1>& A, 01186 Vector<complex<double>,VectFull, Allocator2>& w, 01187 Matrix<complex<double>, 01188 General, RowMajor, Allocator3>& z, 01189 LapackInfo& info) 01190 { 01191 int n = A.GetM(); 01192 Matrix<complex<double>, General, RowMajor> B(n,n); 01193 for (int i = 0; i < n; i++) 01194 for (int j = 0; j < n; j++) 01195 B(i,j) = A(i,j); 01196 01197 w.Reallocate(n); 01198 z.Reallocate(n,n); 01199 GetEigenvaluesEigenvectors(B, w, z); 01200 } 01201 01202 01203 /* ColSymPacked */ 01204 01205 01206 template<class Prop, class Allocator1, class Allocator2> 01207 void GetEigenvalues(Matrix<float, Prop, ColSymPacked, Allocator1>& A, 01208 Vector<float, VectFull, Allocator2>& w, 01209 LapackInfo& info) 01210 { 01211 int n = A.GetM(); 01212 char uplo('U'); char job('N'); 01213 Vector<float> work(3*n); 01214 w.Reallocate(n); 01215 sspev_(&job, &uplo, &n, A.GetData(), w.GetData(), A.GetData(), 01216 &n, work.GetData() , &info.GetInfoRef()); 01217 01218 #ifdef SELDON_LAPACK_CHECK_INFO 01219 if (info.GetInfo() != 0) 01220 throw LapackError(info.GetInfo(), "GetEigenvalues", 01221 "Failed to find eigenvalues "); 01222 #endif 01223 01224 } 01225 01226 01227 template<class Prop, class Allocator1, class Allocator2, class Allocator3> 01228 void GetEigenvaluesEigenvectors(Matrix<float, Prop,ColSymPacked, Allocator1>& A, 01229 Vector<float, VectFull, Allocator2>& w, 01230 Matrix<float, General, ColMajor, Allocator3>&z, 01231 LapackInfo& info) 01232 { 01233 int n = A.GetM(); 01234 char uplo('U'); char job('V'); 01235 Vector<float> work(3*n); 01236 w.Reallocate(n); 01237 z.Reallocate(n,n); 01238 sspev_(&job, &uplo, &n, A.GetData(), w.GetData(), z.GetData(), 01239 &n, work.GetData() , &info.GetInfoRef()); 01240 01241 #ifdef SELDON_LAPACK_CHECK_INFO 01242 if (info.GetInfo() != 0) 01243 throw LapackError(info.GetInfo(), "GetEigenvalues", 01244 "Failed to find eigenvalues "); 01245 #endif 01246 01247 } 01248 01249 01250 template<class Prop, class Allocator1, class Allocator2> 01251 void GetEigenvalues(Matrix<complex<float>, 01252 Prop, ColSymPacked, Allocator1>& A, 01253 Vector<complex<float>, VectFull, Allocator2>& w, 01254 LapackInfo& info) 01255 { 01256 int n = A.GetM(); 01257 Matrix<complex<float>, General, ColMajor> B(n,n); 01258 for (int i = 0; i < n; i++) 01259 for (int j = 0; j < n; j++) 01260 B(i,j) = A(i,j); 01261 01262 w.Reallocate(n); 01263 GetEigenvalues(B, w); 01264 } 01265 01266 01267 template<class Prop, class Allocator1, class Allocator2, class Allocator3> 01268 void GetEigenvaluesEigenvectors(Matrix<complex<float>, 01269 Prop, ColSymPacked, Allocator1>& A, 01270 Vector<complex<float>, VectFull, Allocator2>& w, 01271 Matrix<complex<float>, 01272 General, ColMajor, Allocator3>& z, 01273 LapackInfo& info) 01274 { 01275 int n = A.GetM(); 01276 Matrix<complex<float>, General, ColMajor> B(n,n); 01277 for (int i = 0; i < n; i++) 01278 for (int j = 0; j < n; j++) 01279 B(i,j) = A(i,j); 01280 01281 w.Reallocate(n); 01282 z.Reallocate(n,n); 01283 GetEigenvaluesEigenvectors(B, w, z); 01284 } 01285 01286 01287 template<class Prop, class Allocator1, class Allocator2> 01288 void GetEigenvalues(Matrix<double, Prop, ColSymPacked, Allocator1>& A, 01289 Vector<double, VectFull, Allocator2>& w, 01290 LapackInfo& info) 01291 { 01292 int n = A.GetM(); 01293 char uplo('U'); char job('N'); Vector<double> work(3*n); 01294 w.Reallocate(n); 01295 dspev_(&job, &uplo, &n, A.GetData(), w.GetData(), A.GetData(), 01296 &n, work.GetData() , &info.GetInfoRef()); 01297 01298 #ifdef SELDON_LAPACK_CHECK_INFO 01299 if (info.GetInfo() != 0) 01300 throw LapackError(info.GetInfo(), "GetEigenvalues", 01301 "Failed to find eigenvalues "); 01302 #endif 01303 01304 } 01305 01306 01307 template<class Prop, class Allocator1, class Allocator2, class Allocator3> 01308 void GetEigenvaluesEigenvectors(Matrix<double,Prop,ColSymPacked, Allocator1>& A, 01309 Vector<double, VectFull, Allocator2>& w, 01310 Matrix<double, General, ColMajor, Allocator3>&z, 01311 LapackInfo& info) 01312 { 01313 int n = A.GetM(); 01314 char uplo('U'); char job('V'); 01315 Vector<double> work(3*n); 01316 w.Reallocate(n); 01317 z.Reallocate(n,n); 01318 dspev_(&job, &uplo, &n, A.GetData(), w.GetData(), z.GetData(), 01319 &n, work.GetData() , &info.GetInfoRef()); 01320 01321 #ifdef SELDON_LAPACK_CHECK_INFO 01322 if (info.GetInfo() != 0) 01323 throw LapackError(info.GetInfo(), "GetEigenvalues", 01324 "Failed to find eigenvalues "); 01325 #endif 01326 01327 } 01328 01329 01330 template<class Prop, class Allocator1, class Allocator2> 01331 void GetEigenvalues(Matrix<complex<double>, 01332 Prop, ColSymPacked, Allocator1>& A, 01333 Vector<complex<double>, VectFull, Allocator2>& w, 01334 LapackInfo& info) 01335 { 01336 int n = A.GetM(); 01337 Matrix<complex<double>, General, ColMajor> B(n,n); 01338 for (int i = 0; i < n; i++) 01339 for (int j = 0; j < n; j++) 01340 B(i,j) = A(i,j); 01341 01342 w.Reallocate(n); 01343 GetEigenvalues(B, w); 01344 } 01345 01346 01347 template<class Prop, class Allocator1, class Allocator2, class Allocator3> 01348 void GetEigenvaluesEigenvectors(Matrix<complex<double>, 01349 Prop, ColSymPacked, Allocator1>& A, 01350 Vector<complex<double>, VectFull, Allocator2>& w, 01351 Matrix<complex<double>, 01352 General, ColMajor, Allocator3>& z, 01353 LapackInfo& info) 01354 { 01355 int n = A.GetM(); 01356 Matrix<complex<double>, General, ColMajor> B(n,n); 01357 for (int i = 0; i < n; i++) 01358 for (int j = 0; j < n; j++) 01359 B(i,j) = A(i,j); 01360 01361 w.Reallocate(n); 01362 z.Reallocate(n,n); 01363 GetEigenvaluesEigenvectors(B, w, z); 01364 } 01365 01366 01367 /* RowHermPacked */ 01368 01369 01370 template<class Prop, class Allocator1, class Allocator2> 01371 void GetEigenvalues(Matrix<complex<float>, 01372 Prop, RowHermPacked, Allocator1>& A, 01373 Vector<float, VectFull, Allocator2>& w, 01374 LapackInfo& info) 01375 { 01376 int n = A.GetM(); 01377 char uplo('L'); char job('N'); 01378 Vector<complex<float> > work(2*n); 01379 Vector<float> rwork(3*n); 01380 w.Reallocate(n); 01381 chpev_(&job, &uplo, &n, A.GetDataVoid(), w.GetData(), A.GetDataVoid(), &n, 01382 work.GetDataVoid(), rwork.GetData(), &info.GetInfoRef()); 01383 01384 #ifdef SELDON_LAPACK_CHECK_INFO 01385 if (info.GetInfo() != 0) 01386 throw LapackError(info.GetInfo(), "GetEigenvalues", 01387 "Failed to find eigenvalues "); 01388 #endif 01389 01390 } 01391 01392 01393 template<class Prop, class Allocator1, class Allocator2, class Allocator3> 01394 void GetEigenvaluesEigenvectors(Matrix<complex<float>, 01395 Prop, RowHermPacked, Allocator1>& A, 01396 Vector<float, VectFull, Allocator2>& w, 01397 Matrix<complex<float>, 01398 General, RowMajor, Allocator3>&z, 01399 LapackInfo& info) 01400 { 01401 int n = A.GetM(); 01402 char uplo('L'); char job('V'); 01403 Vector<complex<float> > work(2*n); 01404 Vector<float> rwork(3*n); 01405 w.Reallocate(n); 01406 z.Reallocate(n,n); 01407 chpev_(&job, &uplo, &n, A.GetDataVoid(), w.GetData(), z.GetDataVoid(), 01408 &n, work.GetDataVoid(), rwork.GetData(), &info.GetInfoRef()); 01409 01410 #ifdef SELDON_LAPACK_CHECK_INFO 01411 if (info.GetInfo() != 0) 01412 throw LapackError(info.GetInfo(), "GetEigenvalues", 01413 "Failed to find eigenvalues "); 01414 #endif 01415 01416 Transpose(z); 01417 } 01418 01419 01420 template<class Prop, class Allocator1, class Allocator2> 01421 void GetEigenvalues(Matrix<complex<double>, 01422 Prop, RowHermPacked, Allocator1>& A, 01423 Vector<double, VectFull, Allocator2>& w, 01424 LapackInfo& info) 01425 { 01426 int n = A.GetM(); 01427 char uplo('L'); char job('N'); 01428 Vector<complex<double> > work(2*n); 01429 Vector<double> rwork(3*n); 01430 w.Reallocate(n); 01431 zhpev_(&job, &uplo, &n, A.GetDataVoid(), w.GetData(), A.GetDataVoid(), &n, 01432 work.GetDataVoid(), rwork.GetData(), &info.GetInfoRef()); 01433 01434 #ifdef SELDON_LAPACK_CHECK_INFO 01435 if (info.GetInfo() != 0) 01436 throw LapackError(info.GetInfo(), "GetEigenvalues", 01437 "Failed to find eigenvalues "); 01438 #endif 01439 01440 } 01441 01442 01443 template<class Prop, class Allocator1, class Allocator2, class Allocator3> 01444 void GetEigenvaluesEigenvectors(Matrix<complex<double>, 01445 Prop, RowHermPacked, Allocator1>& A, 01446 Vector<double, VectFull, Allocator2>& w, 01447 Matrix<complex<double>, 01448 General, RowMajor, Allocator3>&z, 01449 LapackInfo& info) 01450 { 01451 int n = A.GetM(); 01452 char uplo('L'); char job('V'); 01453 Vector<complex<double> > work(2*n); 01454 Vector<double> rwork(3*n); 01455 w.Reallocate(n); 01456 z.Reallocate(n,n); 01457 zhpev_(&job, &uplo, &n, A.GetDataVoid(), w.GetData(), z.GetDataVoid(), 01458 &n, work.GetDataVoid(), rwork.GetData(), &info.GetInfoRef()); 01459 01460 #ifdef SELDON_LAPACK_CHECK_INFO 01461 if (info.GetInfo() != 0) 01462 throw LapackError(info.GetInfo(), "GetEigenvalues", 01463 "Failed to find eigenvalues "); 01464 #endif 01465 01466 Transpose(z); 01467 } 01468 01469 01470 /* ColHermPacked */ 01471 01472 01473 template<class Prop, class Allocator1, class Allocator2> 01474 void GetEigenvalues(Matrix<complex<float>, 01475 Prop, ColHermPacked, Allocator1>& A, 01476 Vector<float, VectFull, Allocator2>& w, 01477 LapackInfo& info) 01478 { 01479 int n = A.GetM(); 01480 char uplo('U'); 01481 char job('N'); 01482 Vector<complex<float> > work(2*n); 01483 Vector<float> rwork(3*n); 01484 w.Reallocate(n); 01485 chpev_(&job, &uplo, &n, A.GetDataVoid(), w.GetData(), A.GetDataVoid(), &n, 01486 work.GetDataVoid(), rwork.GetData(), &info.GetInfoRef()); 01487 01488 #ifdef SELDON_LAPACK_CHECK_INFO 01489 if (info.GetInfo() != 0) 01490 throw LapackError(info.GetInfo(), "GetEigenvalues", 01491 "Failed to find eigenvalues "); 01492 #endif 01493 01494 } 01495 01496 01497 template<class Prop, class Allocator1, class Allocator2, class Allocator3> 01498 void GetEigenvaluesEigenvectors(Matrix<complex<float>, 01499 Prop, ColHermPacked, Allocator1>& A, 01500 Vector<float, VectFull, Allocator2>& w, 01501 Matrix<complex<float>, 01502 General, ColMajor, Allocator3>&z, 01503 LapackInfo& info) 01504 { 01505 int n = A.GetM(); 01506 char uplo('U'); char job('V'); 01507 Vector<complex<float> > work(2*n); 01508 Vector<float> rwork(3*n); 01509 w.Reallocate(n); 01510 z.Reallocate(n,n); 01511 chpev_(&job, &uplo, &n, A.GetDataVoid(), w.GetData(), z.GetDataVoid(), &n, 01512 work.GetDataVoid(), rwork.GetData(), &info.GetInfoRef()); 01513 01514 #ifdef SELDON_LAPACK_CHECK_INFO 01515 if (info.GetInfo() != 0) 01516 throw LapackError(info.GetInfo(), "GetEigenvalues", 01517 "Failed to find eigenvalues "); 01518 #endif 01519 01520 } 01521 01522 01523 template<class Prop, class Allocator1, class Allocator2> 01524 void GetEigenvalues(Matrix<complex<double>, 01525 Prop, ColHermPacked, Allocator1>& A, 01526 Vector<double, VectFull, Allocator2>& w, 01527 LapackInfo& info) 01528 { 01529 int n = A.GetM(); 01530 char uplo('U'); 01531 char job('N'); 01532 Vector<complex<double> > work(2*n); 01533 Vector<double> rwork(3*n); 01534 w.Reallocate(n); 01535 zhpev_(&job, &uplo, &n, A.GetDataVoid(), w.GetData(), A.GetDataVoid(), &n, 01536 work.GetDataVoid(), rwork.GetData(), &info.GetInfoRef()); 01537 01538 #ifdef SELDON_LAPACK_CHECK_INFO 01539 if (info.GetInfo() != 0) 01540 throw LapackError(info.GetInfo(), "GetEigenvalues", 01541 "Failed to find eigenvalues "); 01542 #endif 01543 01544 } 01545 01546 01547 template<class Prop, class Allocator1, class Allocator2, class Allocator3> 01548 void GetEigenvaluesEigenvectors(Matrix<complex<double>, 01549 Prop, ColHermPacked, Allocator1>& A, 01550 Vector<double, VectFull, Allocator2>& w, 01551 Matrix<complex<double>, 01552 General, ColMajor, Allocator3>&z, 01553 LapackInfo& info) 01554 { 01555 int n = A.GetM(); 01556 char uplo('U'); char job('V'); 01557 Vector<complex<double> > work(2*n); 01558 Vector<double> rwork(3*n); 01559 w.Reallocate(n); 01560 z.Reallocate(n,n); 01561 zhpev_(&job, &uplo, &n, A.GetDataVoid(), w.GetData(), z.GetDataVoid(), &n, 01562 work.GetDataVoid(), rwork.GetData(), &info.GetInfoRef()); 01563 01564 #ifdef SELDON_LAPACK_CHECK_INFO 01565 if (info.GetInfo() != 0) 01566 throw LapackError(info.GetInfo(), "GetEigenvalues", 01567 "Failed to find eigenvalues "); 01568 #endif 01569 01570 } 01571 01572 01573 // STANDARD EIGENVALUE PROBLEM // 01575 01576 01578 // GENERALIZED EIGENVALUE PROBLEM // 01579 01580 01581 /* RowSym */ 01582 01583 01584 template<class Prop1, class Prop2, class Allocator1, 01585 class Allocator2, class Allocator3> 01586 void GetEigenvalues(Matrix<float, Prop1, RowSym, Allocator1>& A, 01587 Matrix<float, Prop2, RowSym, Allocator2>& B, 01588 Vector<float, VectFull, Allocator3>& w, 01589 LapackInfo& info) 01590 { 01591 int itype = 1; 01592 int n = A.GetM(); 01593 char uplo('L'); char job('N'); 01594 int lwork = 3*n; Vector<float> work(lwork); 01595 w.Reallocate(n); 01596 ssygv_(&itype, &job, &uplo, &n, A.GetData(), &n, B.GetData(), &n, 01597 w.GetData(), work.GetData(), &lwork, &info.GetInfoRef()); 01598 01599 #ifdef SELDON_LAPACK_CHECK_INFO 01600 if (info.GetInfo() != 0) 01601 throw LapackError(info.GetInfo(), "GetEigenvalues", 01602 "Failed to find eigenvalues "); 01603 #endif 01604 01605 } 01606 01607 01608 template<class Prop1, class Prop2, class Allocator1, class Allocator2, 01609 class Allocator3, class Allocator4> 01610 void GetEigenvaluesEigenvectors(Matrix<float, Prop1, RowSym, Allocator1>& A, 01611 Matrix<float, Prop2, RowSym, Allocator2>& B, 01612 Vector<float, VectFull, Allocator3>& w, 01613 Matrix<float, General, RowMajor, Allocator4>& z, 01614 LapackInfo& info) 01615 { 01616 int itype = 1; 01617 int n = A.GetM(); 01618 char uplo('L'); char job('V'); 01619 int lwork = 3*n; Vector<float> work(lwork); 01620 w.Reallocate(n); 01621 z.Reallocate(n, n); 01622 for (int i = 0; i < n; i++) 01623 for (int j = 0; j < n; j++) 01624 z(i,j) = A(i,j); 01625 01626 ssygv_(&itype, &job, &uplo, &n, z.GetData(), &n, B.GetData(), &n, 01627 w.GetData(), work.GetData(), &lwork, &info.GetInfoRef()); 01628 01629 #ifdef SELDON_LAPACK_CHECK_INFO 01630 if (info.GetInfo() != 0) 01631 throw LapackError(info.GetInfo(), "GetEigenvalues", 01632 "Failed to find eigenvalues "); 01633 #endif 01634 01635 } 01636 01637 01638 template<class Prop1, class Prop2, class Allocator1, 01639 class Allocator2, class Allocator4, class Allocator5> 01640 void GetEigenvalues(Matrix<complex<float>, Prop1, RowSym, Allocator1>& A, 01641 Matrix<complex<float>, Prop2, RowSym, Allocator2>& B, 01642 Vector<complex<float>, VectFull, Allocator4>& alpha, 01643 Vector<complex<float>, VectFull, Allocator5>& beta, 01644 LapackInfo& info) 01645 { 01646 int n = A.GetM(); 01647 char jobvl('N'), jobvr('N'); 01648 int lwork = 2*n; Vector<complex<float> > work(lwork); 01649 Vector<float> rwork(8*n); 01650 alpha.Reallocate(n); 01651 beta.Reallocate(n); 01652 cggev_(&jobvl, &jobvr, &n, A.GetData(), &n, B.GetData(), &n, 01653 alpha.GetData(), beta.GetData(), A.GetData(), &n, A.GetData(), &n, 01654 work.GetData(), &lwork, rwork.GetData(), &info.GetInfoRef()); 01655 01656 #ifdef SELDON_LAPACK_CHECK_INFO 01657 if (info.GetInfo() != 0) 01658 throw LapackError(info.GetInfo(), "GetEigenvalues", 01659 "Failed to find eigenvalues "); 01660 #endif 01661 01662 } 01663 01664 01665 template<class Prop1, class Prop2, class Prop3, class Allocator1, 01666 class Allocator2, class Allocator4, 01667 class Allocator5, class Allocator6> 01668 void GetEigenvaluesEigenvectors(Matrix<complex<float>, 01669 Prop1, RowSym, Allocator1>& A, 01670 Matrix<complex<float>, 01671 Prop2, RowSym, Allocator2>& B, 01672 Vector<complex<float>, 01673 VectFull, Allocator4>& alpha, 01674 Vector<complex<float>, 01675 VectFull, Allocator5>& beta, 01676 Matrix<complex<float>, 01677 Prop3, RowMajor, Allocator6>& V, 01678 LapackInfo& info) 01679 { 01680 int n = A.GetM(); 01681 char jobvl('V'), jobvr('N'); 01682 int lwork = 2*n; Vector<complex<float> > work(lwork); 01683 Vector<float> rwork(8*n); 01684 alpha.Reallocate(n); 01685 beta.Reallocate(n); 01686 V.Reallocate(n); 01687 cggev_(&jobvl, &jobvr, &n, A.GetData(), &n, B.GetData(), &n, 01688 alpha.GetData(), beta.GetData(), V.GetData(), &n, V.GetData(), &n, 01689 work.GetData(), &lwork, rwork.GetData(), &info.GetInfoRef()); 01690 01691 #ifdef SELDON_LAPACK_CHECK_INFO 01692 if (info.GetInfo() != 0) 01693 throw LapackError(info.GetInfo(), "GetEigenvalues", 01694 "Failed to find eigenvalues "); 01695 #endif 01696 01697 TransposeConj(V); 01698 } 01699 01700 01701 template<class Prop1, class Prop2, class Allocator1, 01702 class Allocator2, class Allocator3> 01703 void GetEigenvalues(Matrix<double, Prop1, RowSym, Allocator1>& A, 01704 Matrix<double, Prop2, RowSym, Allocator2>& B, 01705 Vector<double, VectFull, Allocator3>& w, 01706 LapackInfo& info) 01707 { 01708 int itype = 1; 01709 int n = A.GetM(); 01710 char uplo('L'); char job('N'); 01711 int lwork = 3*n; Vector<double> work(lwork); 01712 w.Reallocate(n); 01713 dsygv_(&itype, &job, &uplo, &n, A.GetData(), &n, B.GetData(), &n, 01714 w.GetData(), work.GetData(), &lwork, &info.GetInfoRef()); 01715 01716 #ifdef SELDON_LAPACK_CHECK_INFO 01717 if (info.GetInfo() != 0) 01718 throw LapackError(info.GetInfo(), "GetEigenvalues", 01719 "Failed to find eigenvalues "); 01720 #endif 01721 01722 } 01723 01724 01725 template<class Prop1, class Prop2, class Allocator1, class Allocator2, 01726 class Allocator3, class Allocator4> 01727 void GetEigenvaluesEigenvectors(Matrix<double, Prop1, RowSym, Allocator1>& A, 01728 Matrix<double, Prop2, RowSym, Allocator2>& B, 01729 Vector<double, VectFull, Allocator3>& w, 01730 Matrix<double, General, RowMajor, Allocator4>& z, 01731 LapackInfo& info) 01732 { 01733 int itype = 1; 01734 int n = A.GetM(); 01735 char uplo('L'); char job('V'); 01736 int lwork = 3*n; Vector<double> work(lwork); 01737 w.Reallocate(n); 01738 z.Reallocate(n,n); 01739 for (int i = 0; i < n; i++) 01740 for (int j = 0; j < n; j++) 01741 z(i,j) = A(i,j); 01742 01743 dsygv_(&itype, &job, &uplo, &n, z.GetData(), &n, B.GetData(), &n, 01744 w.GetData(), work.GetData(), &lwork, &info.GetInfoRef()); 01745 01746 #ifdef SELDON_LAPACK_CHECK_INFO 01747 if (info.GetInfo() != 0) 01748 throw LapackError(info.GetInfo(), "GetEigenvalues", 01749 "Failed to find eigenvalues "); 01750 #endif 01751 01752 } 01753 01754 01755 template<class Prop1, class Prop2, class Allocator1, 01756 class Allocator2, class Allocator4, class Allocator5> 01757 void GetEigenvalues(Matrix<complex<double>, Prop1, RowSym, Allocator1>& A, 01758 Matrix<complex<double>, Prop2, RowSym, Allocator2>& B, 01759 Vector<complex<double>, VectFull, Allocator4>& alpha, 01760 Vector<complex<double>, VectFull, Allocator5>& beta, 01761 LapackInfo& info) 01762 { 01763 int n = A.GetM(); 01764 char jobvl('N'), jobvr('N'); 01765 int lwork = 2*n; Vector<complex<double> > work(lwork); 01766 Vector<double> rwork(8*n); 01767 alpha.Reallocate(n); 01768 beta.Reallocate(n); 01769 zggev_(&jobvl, &jobvr, &n, A.GetData(), &n, B.GetData(), &n, 01770 alpha.GetData(), beta.GetData(), A.GetData(), &n, A.GetData(), &n, 01771 work.GetData(), &lwork, rwork.GetData(), &info.GetInfoRef()); 01772 01773 #ifdef SELDON_LAPACK_CHECK_INFO 01774 if (info.GetInfo() != 0) 01775 throw LapackError(info.GetInfo(), "GetEigenvalues", 01776 "Failed to find eigenvalues "); 01777 #endif 01778 01779 } 01780 01781 01782 template<class Prop1, class Prop2, class Prop3, class Allocator1, 01783 class Allocator2, class Allocator4, 01784 class Allocator5, class Allocator6> 01785 void GetEigenvaluesEigenvectors(Matrix<complex<double>, 01786 Prop1, RowSym, Allocator1>& A, 01787 Matrix<complex<double>, 01788 Prop2, RowSym, Allocator2>& B, 01789 Vector<complex<double>, 01790 VectFull, Allocator4>& alpha, 01791 Vector<complex<double>, 01792 VectFull, Allocator5>& beta, 01793 Matrix<complex<double>, 01794 Prop3, RowMajor, Allocator6>& V, 01795 LapackInfo& info) 01796 { 01797 int n = A.GetM(); 01798 char jobvl('V'), jobvr('N'); 01799 int lwork = 2*n; Vector<complex<double> > work(lwork); 01800 Vector<double> rwork(8*n); 01801 alpha.Reallocate(n); 01802 beta.Reallocate(n); 01803 V.Reallocate(n,n); 01804 zggev_(&jobvl, &jobvr, &n, A.GetData(), &n, B.GetData(), &n, 01805 alpha.GetData(), beta.GetData(), V.GetData(), &n, V.GetData(), &n, 01806 work.GetData(), &lwork, rwork.GetData(), &info.GetInfoRef()); 01807 01808 #ifdef SELDON_LAPACK_CHECK_INFO 01809 if (info.GetInfo() != 0) 01810 throw LapackError(info.GetInfo(), "GetEigenvalues", 01811 "Failed to find eigenvalues "); 01812 #endif 01813 01814 TransposeConj(V); 01815 } 01816 01817 01818 /* ColSym */ 01819 01820 01821 template<class Prop1, class Prop2, class Allocator1, 01822 class Allocator2, class Allocator3> 01823 void GetEigenvalues(Matrix<float, Prop1, ColSym, Allocator1>& A, 01824 Matrix<float, Prop2, ColSym, Allocator2>& B, 01825 Vector<float, VectFull, Allocator3>& w, 01826 LapackInfo& info) 01827 { 01828 int itype = 1; 01829 int n = A.GetM(); 01830 char uplo('U'); char job('N'); 01831 int lwork = 3*n; Vector<float> work(lwork); 01832 w.Reallocate(n); 01833 ssygv_(&itype, &job, &uplo, &n, A.GetData(), &n, B.GetData(), &n, 01834 w.GetData(), work.GetData(), &lwork, &info.GetInfoRef()); 01835 01836 #ifdef SELDON_LAPACK_CHECK_INFO 01837 if (info.GetInfo() != 0) 01838 throw LapackError(info.GetInfo(), "GetEigenvalues", 01839 "Failed to find eigenvalues "); 01840 #endif 01841 01842 } 01843 01844 01845 template<class Prop1, class Prop2, class Allocator1, class Allocator2, 01846 class Allocator3, class Allocator4> 01847 void GetEigenvaluesEigenvectors(Matrix<float, Prop1, ColSym, Allocator1>& A, 01848 Matrix<float, Prop2, ColSym, Allocator2>& B, 01849 Vector<float, VectFull, Allocator3>& w, 01850 Matrix<float, General, ColMajor, Allocator4>& z, 01851 LapackInfo& info) 01852 { 01853 int itype = 1; 01854 int n = A.GetM(); 01855 char uplo('U'); char job('V'); 01856 int lwork = 3*n; Vector<float> work(lwork); 01857 w.Reallocate(n); 01858 z.Reallocate(n, n); 01859 for (int i = 0; i < n; i++) 01860 for (int j = 0; j < n; j++) 01861 z(i,j) = A(i,j); 01862 01863 ssygv_(&itype, &job, &uplo, &n, z.GetData(), &n, B.GetData(), &n, 01864 w.GetData(), work.GetData(), &lwork, &info.GetInfoRef()); 01865 01866 #ifdef SELDON_LAPACK_CHECK_INFO 01867 if (info.GetInfo() != 0) 01868 throw LapackError(info.GetInfo(), "GetEigenvalues", 01869 "Failed to find eigenvalues "); 01870 #endif 01871 01872 } 01873 01874 01875 template<class Prop1, class Prop2, class Allocator1, 01876 class Allocator2, class Allocator4, class Allocator5> 01877 void GetEigenvalues(Matrix<complex<float>, Prop1, ColSym, Allocator1>& A, 01878 Matrix<complex<float>, Prop2, ColSym, Allocator2>& B, 01879 Vector<complex<float>, VectFull, Allocator4>& alpha, 01880 Vector<complex<float>, VectFull, Allocator5>& beta, 01881 LapackInfo& info) 01882 { 01883 int n = A.GetM(); 01884 char jobvl('N'), jobvr('N'); int lwork = 2*n; 01885 Vector<complex<float> > work(lwork); 01886 Vector<float> rwork(8*n); 01887 alpha.Reallocate(n); 01888 beta.Reallocate(n); 01889 cggev_(&jobvl, &jobvr, &n, A.GetData(), &n, B.GetData(), &n, 01890 alpha.GetData(), beta.GetData(), A.GetData(), &n, A.GetData(), &n, 01891 work.GetData(), &lwork, rwork.GetData(), &info.GetInfoRef()); 01892 01893 #ifdef SELDON_LAPACK_CHECK_INFO 01894 if (info.GetInfo() != 0) 01895 throw LapackError(info.GetInfo(), "GetEigenvalues", 01896 "Failed to find eigenvalues "); 01897 #endif 01898 01899 } 01900 01901 01902 template<class Prop1, class Prop2, class Prop3, class Alloc1, 01903 class Alloc2, class Alloc4, class Alloc5, class Alloc6> 01904 void GetEigenvaluesEigenvectors(Matrix<complex<float>, 01905 Prop1, ColSym, Alloc1>& A, 01906 Matrix<complex<float>, 01907 Prop2, ColSym, Alloc2>& B, 01908 Vector<complex<float>, VectFull, Alloc4>& alpha, 01909 Vector<complex<float>, VectFull, Alloc5>& beta, 01910 Matrix<complex<float>, 01911 Prop3, ColMajor, Alloc6>& V, 01912 LapackInfo& info) 01913 { 01914 int n = A.GetM(); 01915 char jobvl('N'), jobvr('V'); 01916 int lwork = 2*n; Vector<complex<float> > work(lwork); 01917 Vector<float> rwork(8*n); 01918 alpha.Reallocate(n); 01919 beta.Reallocate(n); 01920 V.Reallocate(n,n); 01921 cggev_(&jobvl, &jobvr, &n, A.GetData(), &n, B.GetData(), &n, 01922 alpha.GetData(), beta.GetData(), V.GetData(), &n, V.GetData(), &n, 01923 work.GetData(), &lwork, rwork.GetData(), &info.GetInfoRef()); 01924 01925 #ifdef SELDON_LAPACK_CHECK_INFO 01926 if (info.GetInfo() != 0) 01927 throw LapackError(info.GetInfo(), "GetEigenvalues", 01928 "Failed to find eigenvalues "); 01929 #endif 01930 01931 } 01932 01933 template<class Prop1, class Prop2, class Allocator1, 01934 class Allocator2, class Allocator3> 01935 void GetEigenvalues(Matrix<double, Prop1, ColSym, Allocator1>& A, 01936 Matrix<double, Prop2, ColSym, Allocator2>& B, 01937 Vector<double, VectFull, Allocator3>& w, 01938 LapackInfo& info) 01939 { 01940 int itype = 1; 01941 int n = A.GetM(); 01942 char uplo('U'); char job('N'); 01943 int lwork = 3*n; Vector<double> work(lwork); 01944 w.Reallocate(n); 01945 dsygv_(&itype, &job, &uplo, &n, A.GetData(), &n, B.GetData(), &n, 01946 w.GetData(), work.GetData(), &lwork, &info.GetInfoRef()); 01947 01948 #ifdef SELDON_LAPACK_CHECK_INFO 01949 if (info.GetInfo() != 0) 01950 throw LapackError(info.GetInfo(), "GetEigenvalues", 01951 "Failed to find eigenvalues "); 01952 #endif 01953 01954 } 01955 01956 01957 template<class Prop1, class Prop2, class Allocator1, class Allocator2, 01958 class Allocator3, class Allocator4> 01959 void GetEigenvaluesEigenvectors(Matrix<double, Prop1, ColSym, Allocator1>& A, 01960 Matrix<double, Prop2, ColSym, Allocator2>& B, 01961 Vector<double, VectFull, Allocator3>& w, 01962 Matrix<double, General, ColMajor, Allocator4>& z, 01963 LapackInfo& info) 01964 { 01965 int itype = 1; 01966 int n = A.GetM(); 01967 char uplo('U'); char job('V'); 01968 int lwork = 3*n; Vector<double> work(lwork); 01969 w.Reallocate(n); 01970 z.Reallocate(n,n); 01971 for (int i = 0; i < n; i++) 01972 for (int j = 0; j < n; j++) 01973 z(i,j) = A(i,j); 01974 01975 dsygv_(&itype, &job, &uplo, &n, z.GetData(), &n, B.GetData(), &n, 01976 w.GetData(), work.GetData(), &lwork, &info.GetInfoRef()); 01977 01978 #ifdef SELDON_LAPACK_CHECK_INFO 01979 if (info.GetInfo() != 0) 01980 throw LapackError(info.GetInfo(), "GetEigenvalues", 01981 "Failed to find eigenvalues "); 01982 #endif 01983 01984 } 01985 01986 01987 template<class Prop1, class Prop2, class Allocator1, 01988 class Allocator2, class Allocator4, class Allocator5> 01989 void GetEigenvalues(Matrix<complex<double>, Prop1, ColSym, Allocator1>& A, 01990 Matrix<complex<double>, Prop2, ColSym, Allocator2>& B, 01991 Vector<complex<double>, VectFull, Allocator4>& alpha, 01992 Vector<complex<double>, VectFull, Allocator5>& beta, 01993 LapackInfo& info) 01994 { 01995 int n = A.GetM(); 01996 char jobvl('N'), jobvr('N'); int lwork = 2*n; 01997 Vector<complex<double> > work(lwork); 01998 Vector<double> rwork(8*n); 01999 alpha.Reallocate(n); 02000 beta.Reallocate(n); 02001 zggev_(&jobvl, &jobvr, &n, A.GetData(), &n, B.GetData(), &n, 02002 alpha.GetData(), beta.GetData(), A.GetData(), &n, A.GetData(), &n, 02003 work.GetData(), &lwork, rwork.GetData(), &info.GetInfoRef()); 02004 02005 #ifdef SELDON_LAPACK_CHECK_INFO 02006 if (info.GetInfo() != 0) 02007 throw LapackError(info.GetInfo(), "GetEigenvalues", 02008 "Failed to find eigenvalues "); 02009 #endif 02010 02011 } 02012 02013 02014 template<class Prop1, class Prop2, class Prop3, class Alloc1, 02015 class Alloc2, class Alloc4, class Alloc5, class Alloc6> 02016 void GetEigenvaluesEigenvectors(Matrix<complex<double>, 02017 Prop1, ColSym, Alloc1>& A, 02018 Matrix<complex<double>, 02019 Prop2, ColSym, Alloc2>& B, 02020 Vector<complex<double>, VectFull, Alloc4>& alpha, 02021 Vector<complex<double>, VectFull, Alloc5>& beta, 02022 Matrix<complex<double>, 02023 Prop3, ColMajor, Alloc6>& V, 02024 LapackInfo& info) 02025 { 02026 int n = A.GetM(); 02027 char jobvl('N'), jobvr('V'); 02028 int lwork = 2*n; Vector<complex<double> > work(lwork); 02029 Vector<double> rwork(8*n); 02030 alpha.Reallocate(n); 02031 beta.Reallocate(n); 02032 V.Reallocate(n,n); 02033 zggev_(&jobvl, &jobvr, &n, A.GetData(), &n, B.GetData(), &n, 02034 alpha.GetData(), beta.GetData(), V.GetData(), &n, V.GetData(), &n, 02035 work.GetData(), &lwork, rwork.GetData(), &info.GetInfoRef()); 02036 02037 #ifdef SELDON_LAPACK_CHECK_INFO 02038 if (info.GetInfo() != 0) 02039 throw LapackError(info.GetInfo(), "GetEigenvalues", 02040 "Failed to find eigenvalues "); 02041 #endif 02042 02043 } 02044 02045 02046 /* RowHerm */ 02047 02048 02049 template<class Prop1, class Prop2, class Allocator1, 02050 class Allocator2, class Allocator3> 02051 void GetEigenvalues(Matrix<complex<float>, Prop1, RowHerm, Allocator1>& A, 02052 Matrix<complex<float>, Prop2, RowHerm, Allocator2>& B, 02053 Vector<float, VectFull, Allocator3>& w, 02054 LapackInfo& info) 02055 { 02056 int itype = 1; 02057 int n = A.GetM(); 02058 char uplo('L'); char job('N'); 02059 int lwork = 2*n; Vector<float> work(lwork); 02060 Vector<float> rwork(3*n); 02061 w.Reallocate(n); 02062 chegv_(&itype, &job, &uplo, &n, A.GetData(), &n, B.GetData(), &n, 02063 w.GetData(), work.GetData(), &lwork, rwork.GetData(), 02064 &info.GetInfoRef()); 02065 02066 #ifdef SELDON_LAPACK_CHECK_INFO 02067 if (info.GetInfo() != 0) 02068 throw LapackError(info.GetInfo(), "GetEigenvalues", 02069 "Failed to find eigenvalues "); 02070 #endif 02071 02072 } 02073 02074 02075 template<class Prop1, class Prop2, class Allocator1, class Allocator2, 02076 class Allocator3, class Allocator4> 02077 void GetEigenvaluesEigenvectors(Matrix<complex<float>, 02078 Prop1, RowHerm, Allocator1>& A, 02079 Matrix<complex<float>, 02080 Prop2, RowHerm, Allocator2>& B, 02081 Vector<float, VectFull, Allocator3>& w, 02082 Matrix<complex<float>, 02083 General, RowMajor, Allocator4>& z, 02084 LapackInfo& info) 02085 { 02086 int itype = 1; 02087 int n = A.GetM(); 02088 char uplo('L'); char job('V'); 02089 int lwork = 2*n; Vector<float> work(lwork); 02090 w.Reallocate(n); 02091 z.Reallocate(n,n); 02092 for (int i = 0; i < n; i++) 02093 for (int j = 0; j < n; j++) 02094 z(i,j) = A(i,j); 02095 02096 Vector<float> rwork(3*n); 02097 chegv_(&itype, &job, &uplo, &n, z.GetData(), &n, B.GetData(), &n, 02098 w.GetData(), work.GetData(), &lwork, rwork.GetData(), 02099 &info.GetInfoRef()); 02100 02101 #ifdef SELDON_LAPACK_CHECK_INFO 02102 if (info.GetInfo() != 0) 02103 throw LapackError(info.GetInfo(), "GetEigenvalues", 02104 "Failed to find eigenvalues "); 02105 #endif 02106 02107 } 02108 02109 02110 template<class Prop1, class Prop2, class Allocator1, 02111 class Allocator2, class Allocator3> 02112 void GetEigenvalues(Matrix<complex<double>, Prop1, RowHerm, Allocator1>& A, 02113 Matrix<complex<double>, Prop2, RowHerm, Allocator2>& B, 02114 Vector<double, VectFull, Allocator3>& w, 02115 LapackInfo& info) 02116 { 02117 int itype = 1; 02118 int n = A.GetM(); 02119 char uplo('L'); char job('N'); 02120 int lwork = 3*n; Vector<double> work(lwork); 02121 Vector<double> rwork(3*n); 02122 w.Reallocate(n); 02123 zhegv_(&itype, &job, &uplo, &n, A.GetData(), &n, B.GetData(), &n, 02124 w.GetData(), work.GetData(), &lwork, rwork.GetData(), 02125 &info.GetInfoRef()); 02126 02127 #ifdef SELDON_LAPACK_CHECK_INFO 02128 if (info.GetInfo() != 0) 02129 throw LapackError(info.GetInfo(), "GetEigenvalues", 02130 "Failed to find eigenvalues "); 02131 #endif 02132 02133 } 02134 02135 02136 template<class Prop1, class Prop2, class Allocator1, class Allocator2, 02137 class Allocator3, class Allocator4> 02138 void GetEigenvaluesEigenvectors(Matrix<complex<double>, 02139 Prop1, RowHerm, Allocator1>& A, 02140 Matrix<complex<double>, 02141 Prop2, RowHerm, Allocator2>& B, 02142 Vector<double, VectFull, Allocator3>& w, 02143 Matrix<complex<double>, 02144 General, RowMajor, Allocator4>& z, 02145 LapackInfo& info) 02146 { 02147 int itype = 1; 02148 int n = A.GetM(); 02149 char uplo('L'); char job('V'); 02150 int lwork = 3*n; Vector<double> work(lwork); 02151 w.Reallocate(n); 02152 z.Reallocate(n,n); 02153 for (int i = 0; i < n; i++) 02154 for (int j = 0; j < n; j++) 02155 z(i,j) = A(i,j); 02156 02157 Vector<double> rwork(3*n); 02158 zhegv_(&itype, &job, &uplo, &n, z.GetData(), &n, B.GetData(), &n, 02159 w.GetData(), work.GetData(), &lwork, rwork.GetData(), 02160 &info.GetInfoRef()); 02161 02162 #ifdef SELDON_LAPACK_CHECK_INFO 02163 if (info.GetInfo() != 0) 02164 throw LapackError(info.GetInfo(), "GetEigenvalues", 02165 "Failed to find eigenvalues "); 02166 #endif 02167 02168 } 02169 02170 02171 /* ColHerm */ 02172 02173 02174 template<class Prop1, class Prop2, class Allocator1, 02175 class Allocator2, class Allocator3> 02176 void GetEigenvalues(Matrix<complex<float>, Prop1, ColHerm, Allocator1>& A, 02177 Matrix<complex<float>, Prop2, ColHerm, Allocator2>& B, 02178 Vector<float, VectFull, Allocator3>& w, 02179 LapackInfo& info) 02180 { 02181 int itype = 1; 02182 int n = A.GetM(); 02183 char uplo('U'); char job('N'); 02184 int lwork = 2*n; Vector<float> work(lwork); 02185 Vector<float> rwork(3*n); 02186 w.Reallocate(n); 02187 chegv_(&itype, &job, &uplo, &n, A.GetData(), &n, B.GetData(), &n, 02188 w.GetData(), work.GetData(), &lwork, rwork.GetData(), 02189 &info.GetInfoRef()); 02190 02191 #ifdef SELDON_LAPACK_CHECK_INFO 02192 if (info.GetInfo() != 0) 02193 throw LapackError(info.GetInfo(), "GetEigenvalues", 02194 "Failed to find eigenvalues "); 02195 #endif 02196 02197 } 02198 02199 02200 template<class Prop1, class Prop2, class Allocator1, class Allocator2, 02201 class Allocator3, class Allocator4> 02202 void GetEigenvaluesEigenvectors(Matrix<complex<float>, 02203 Prop1, ColHerm, Allocator1>& A, 02204 Matrix<complex<float>, 02205 Prop2, ColHerm, Allocator2>& B, 02206 Vector<float, VectFull, Allocator3>& w, 02207 Matrix<complex<float>, 02208 General, ColMajor, Allocator4>& z, 02209 LapackInfo& info) 02210 { 02211 int itype = 1; 02212 int n = A.GetM(); 02213 char uplo('U'); char job('V'); 02214 int lwork = 3*n; Vector<float> work(lwork); 02215 w.Reallocate(n); 02216 z.Reallocate(n,n); 02217 for (int i = 0; i < n; i++) 02218 for (int j = 0; j < n; j++) 02219 z(i,j) = A(i,j); 02220 02221 Vector<float> rwork(3*n); 02222 chegv_(&itype, &job, &uplo, &n, z.GetData(), &n, B.GetData(), &n, 02223 w.GetData(), work.GetData(), &lwork, 02224 rwork.GetData(), &info.GetInfoRef()); 02225 02226 #ifdef SELDON_LAPACK_CHECK_INFO 02227 if (info.GetInfo() != 0) 02228 throw LapackError(info.GetInfo(), "GetEigenvalues", 02229 "Failed to find eigenvalues "); 02230 #endif 02231 02232 } 02233 02234 02235 template<class Prop1, class Prop2, class Allocator1, 02236 class Allocator2, class Allocator3> 02237 void GetEigenvalues(Matrix<complex<double>, Prop1, ColHerm, Allocator1>& A, 02238 Matrix<complex<double>, Prop2, ColHerm, Allocator2>& B, 02239 Vector<double, VectFull, Allocator3>& w, 02240 LapackInfo& info) 02241 { 02242 int itype = 1; 02243 int n = A.GetM(); 02244 char uplo('U'); char job('N'); 02245 int lwork = 3*n; Vector<double> work(lwork); 02246 Vector<double> rwork(3*n); 02247 w.Reallocate(n); 02248 zhegv_(&itype, &job, &uplo, &n, A.GetData(), &n, B.GetData(), &n, 02249 w.GetData(), work.GetData(), &lwork, rwork.GetData(), 02250 &info.GetInfoRef()); 02251 02252 #ifdef SELDON_LAPACK_CHECK_INFO 02253 if (info.GetInfo() != 0) 02254 throw LapackError(info.GetInfo(), "GetEigenvalues", 02255 "Failed to find eigenvalues "); 02256 #endif 02257 02258 } 02259 02260 02261 template<class Prop1, class Prop2, class Allocator1, class Allocator2, 02262 class Allocator3, class Allocator4> 02263 void GetEigenvaluesEigenvectors(Matrix<complex<double>, 02264 Prop1, ColHerm, Allocator1>& A, 02265 Matrix<complex<double>, 02266 Prop2, ColHerm, Allocator2>& B, 02267 Vector<double, VectFull, Allocator3>& w, 02268 Matrix<complex<double>, 02269 General, ColMajor, Allocator4>& z, 02270 LapackInfo& info) 02271 { 02272 int itype = 1; 02273 int n = A.GetM(); 02274 char uplo('U'); char job('V'); 02275 int lwork = 3*n; Vector<double> work(lwork); 02276 w.Reallocate(n); 02277 z.Reallocate(n,n); 02278 for (int i = 0; i < n; i++) 02279 for (int j = 0; j < n; j++) 02280 z(i,j) = A(i,j); 02281 02282 Vector<double> rwork(3*n); 02283 zhegv_(&itype, &job, &uplo, &n, z.GetData(), &n, B.GetData(), &n, 02284 w.GetData(), work.GetData(), &lwork, 02285 rwork.GetData(), &info.GetInfoRef()); 02286 02287 #ifdef SELDON_LAPACK_CHECK_INFO 02288 if (info.GetInfo() != 0) 02289 throw LapackError(info.GetInfo(), "GetEigenvalues", 02290 "Failed to find eigenvalues "); 02291 #endif 02292 02293 } 02294 02295 02296 /* RowSymPacked */ 02297 02298 02299 template<class Prop1, class Prop2, class Allocator1, 02300 class Allocator2, class Allocator3> 02301 void GetEigenvalues(Matrix<float, Prop1, RowSymPacked, Allocator1>& A, 02302 Matrix<float, Prop2, RowSymPacked, Allocator2>& B, 02303 Vector<float, VectFull, Allocator3>& w, 02304 LapackInfo& info) 02305 { 02306 int itype = 1; 02307 int n = A.GetM(); 02308 char uplo('L'); char job('N'); 02309 int lwork = 3*n; Vector<float> work(lwork); 02310 w.Reallocate(n); 02311 sspgv_(&itype, &job, &uplo, &n, A.GetData(), B.GetData(), w.GetData(), 02312 A.GetData(), &n, work.GetData(), &info.GetInfoRef()); 02313 02314 #ifdef SELDON_LAPACK_CHECK_INFO 02315 if (info.GetInfo() != 0) 02316 throw LapackError(info.GetInfo(), "GetEigenvalues", 02317 "Failed to find eigenvalues "); 02318 #endif 02319 02320 } 02321 02322 02323 template<class Prop1, class Prop2, class Allocator1, class Allocator2, 02324 class Allocator3, class Allocator4> 02325 void GetEigenvaluesEigenvectors(Matrix<float, 02326 Prop1, RowSymPacked, Allocator1>& A, 02327 Matrix<float, 02328 Prop2, RowSymPacked, Allocator2>& B, 02329 Vector<float, VectFull, Allocator3>& w, 02330 Matrix<float, 02331 General, RowMajor, Allocator4>& z, 02332 LapackInfo& info) 02333 { 02334 int itype = 1; 02335 int n = A.GetM(); 02336 char uplo('L'); char job('V'); 02337 int lwork = 3*n; Vector<float> work(lwork); 02338 w.Reallocate(n); 02339 z.Reallocate(n,n); 02340 sspgv_(&itype, &job, &uplo, &n, A.GetData(), B.GetData(), w.GetData(), 02341 z.GetData(), &n, work.GetData(), &info.GetInfoRef()); 02342 02343 #ifdef SELDON_LAPACK_CHECK_INFO 02344 if (info.GetInfo() != 0) 02345 throw LapackError(info.GetInfo(), "GetEigenvalues", 02346 "Failed to find eigenvalues "); 02347 #endif 02348 02349 } 02350 02351 02352 template<class Prop1, class Prop2, class Allocator1, class Allocator2, 02353 class Allocator3, class Allocator4> 02354 void GetEigenvalues(Matrix<complex<float>, 02355 Prop1, RowSymPacked, Allocator1>& A, 02356 Matrix<complex<float>, 02357 Prop2, RowSymPacked, Allocator2>& B, 02358 Vector<complex<float>, VectFull, Allocator3>& alpha, 02359 Vector<complex<float>, VectFull, Allocator4>& beta, 02360 LapackInfo& info) 02361 { 02362 int n = A.GetM(); 02363 Matrix<complex<float>, General, ColMajor> C(n,n), D(n,n); 02364 for (int i = 0; i < n; i++) 02365 for (int j = 0; j < n; j++) 02366 { 02367 C(i,j) = A(i,j); 02368 D(i,j) = B(i,j); 02369 } 02370 02371 alpha.Reallocate(n); 02372 beta.Reallocate(n); 02373 GetEigenvalues(C, D, alpha, beta); 02374 } 02375 02376 02377 template<class Prop1, class Prop2, class Allocator1, class Allocator2, 02378 class Allocator3, class Allocator4, class Allocator5> 02379 void GetEigenvaluesEigenvectors(Matrix<complex<float>, 02380 Prop1, RowSymPacked, Allocator1>& A, 02381 Matrix<complex<float>, 02382 Prop2, RowSymPacked, Allocator2>& B, 02383 Vector<complex<float>, 02384 VectFull, Allocator3>& alpha, 02385 Vector<complex<float>, 02386 VectFull, Allocator4>& beta, 02387 Matrix<complex<float>, 02388 General, RowMajor, Allocator5>& z, 02389 LapackInfo& info) 02390 { 02391 int n = A.GetM(); 02392 Matrix<complex<float>, General, RowMajor> C(n,n), D(n,n); 02393 for (int i = 0; i < n; i++) 02394 for (int j = 0; j < n; j++) 02395 { 02396 C(i,j) = A(i,j); 02397 D(i,j) = B(i,j); 02398 } 02399 02400 alpha.Reallocate(n); 02401 beta.Reallocate(n); 02402 z.Reallocate(n,n); 02403 GetEigenvaluesEigenvectors(C, D, alpha, beta, z); 02404 } 02405 02406 02407 template<class Prop1, class Prop2, class Allocator1, 02408 class Allocator2, class Allocator3> 02409 void GetEigenvalues(Matrix<double, Prop1, RowSymPacked, Allocator1>& A, 02410 Matrix<double, Prop2, RowSymPacked, Allocator2>& B, 02411 Vector<double, VectFull, Allocator3>& w, 02412 LapackInfo& info) 02413 { 02414 int itype = 1; 02415 int n = A.GetM(); 02416 char uplo('L'); char job('N'); 02417 int lwork = 3*n; Vector<double> work(lwork); 02418 w.Reallocate(n); 02419 dspgv_(&itype, &job, &uplo, &n, A.GetData(), B.GetData(), w.GetData(), 02420 A.GetData(), &n, work.GetData(), &info.GetInfoRef()); 02421 02422 #ifdef SELDON_LAPACK_CHECK_INFO 02423 if (info.GetInfo() != 0) 02424 throw LapackError(info.GetInfo(), "GetEigenvalues", 02425 "Failed to find eigenvalues "); 02426 #endif 02427 02428 } 02429 02430 02431 template<class Prop1, class Prop2, class Allocator1, class Allocator2, 02432 class Allocator3, class Allocator4> 02433 void GetEigenvaluesEigenvectors(Matrix<double, 02434 Prop1, RowSymPacked, Allocator1>& A, 02435 Matrix<double, 02436 Prop2, RowSymPacked, Allocator2>& B, 02437 Vector<double, VectFull, Allocator3>& w, 02438 Matrix<double, 02439 General, RowMajor, Allocator4>& z, 02440 LapackInfo& info) 02441 { 02442 int itype = 1; 02443 int n = A.GetM(); 02444 char uplo('L'); char job('V'); 02445 int lwork = 3*n; Vector<double> work(lwork); 02446 w.Reallocate(n); 02447 z.Reallocate(n,n); 02448 dspgv_(&itype, &job, &uplo, &n, A.GetData(), B.GetData(), w.GetData(), 02449 z.GetData(), &n, work.GetData(), &info.GetInfoRef()); 02450 02451 #ifdef SELDON_LAPACK_CHECK_INFO 02452 if (info.GetInfo() != 0) 02453 throw LapackError(info.GetInfo(), "GetEigenvalues", 02454 "Failed to find eigenvalues "); 02455 #endif 02456 02457 } 02458 02459 02460 template<class Prop1, class Prop2, class Allocator1, class Allocator2, 02461 class Allocator3, class Allocator4> 02462 void GetEigenvalues(Matrix<complex<double>, 02463 Prop1, RowSymPacked, Allocator1>& A, 02464 Matrix<complex<double>, 02465 Prop2, RowSymPacked, Allocator2>& B, 02466 Vector<complex<double>, VectFull, Allocator3>& alpha, 02467 Vector<complex<double>, VectFull, Allocator4>& beta, 02468 LapackInfo& info) 02469 { 02470 int n = A.GetM(); 02471 Matrix<complex<double>, General, ColMajor> C(n,n), D(n,n); 02472 for (int i = 0; i < n; i++) 02473 for (int j = 0; j < n; j++) 02474 { 02475 C(i,j) = A(i,j); 02476 D(i,j) = B(i,j); 02477 } 02478 02479 alpha.Reallocate(n); 02480 beta.Reallocate(n); 02481 GetEigenvalues(C, D, alpha, beta); 02482 } 02483 02484 02485 template<class Prop1, class Prop2, class Allocator1, class Allocator2, 02486 class Allocator3, class Allocator4, class Allocator5> 02487 void GetEigenvaluesEigenvectors(Matrix<complex<double>, 02488 Prop1, RowSymPacked, Allocator1>& A, 02489 Matrix<complex<double>, 02490 Prop2, RowSymPacked, Allocator2>& B, 02491 Vector<complex<double>, 02492 VectFull, Allocator3>& alpha, 02493 Vector<complex<double>, 02494 VectFull, Allocator4>& beta, 02495 Matrix<complex<double>, 02496 General, RowMajor, Allocator5>& z, 02497 LapackInfo& info) 02498 { 02499 int n = A.GetM(); 02500 Matrix<complex<double>, General, RowMajor> C(n,n), D(n,n); 02501 for (int i = 0; i < n; i++) 02502 for (int j = 0; j < n; j++) 02503 { 02504 C(i,j) = A(i,j); 02505 D(i,j) = B(i,j); 02506 } 02507 02508 alpha.Reallocate(n); 02509 beta.Reallocate(n); 02510 z.Reallocate(n,n); 02511 GetEigenvaluesEigenvectors(C, D, alpha, beta, z); 02512 } 02513 02514 02515 /* ColSymPacked */ 02516 02517 02518 template<class Prop1, class Prop2, class Allocator1, 02519 class Allocator2, class Allocator3> 02520 void GetEigenvalues(Matrix<float, Prop1, ColSymPacked, Allocator1>& A, 02521 Matrix<float, Prop2, ColSymPacked, Allocator2>& B, 02522 Vector<float, VectFull, Allocator3>& w, 02523 LapackInfo& info) 02524 { 02525 int itype = 1; 02526 int n = A.GetM(); 02527 char uplo('U'); char job('N'); 02528 int lwork = 3*n; Vector<float> work(lwork); 02529 w.Reallocate(n); 02530 sspgv_(&itype, &job, &uplo, &n, A.GetData(), B.GetData(), w.GetData(), 02531 A.GetData(), &n, work.GetData(), &info.GetInfoRef()); 02532 02533 #ifdef SELDON_LAPACK_CHECK_INFO 02534 if (info.GetInfo() != 0) 02535 throw LapackError(info.GetInfo(), "GetEigenvalues", 02536 "Failed to find eigenvalues "); 02537 #endif 02538 02539 } 02540 02541 02542 template<class Prop1, class Prop2, class Allocator1, class Allocator2, 02543 class Allocator3, class Allocator4> 02544 void GetEigenvaluesEigenvectors(Matrix<float, 02545 Prop1, ColSymPacked, Allocator1>& A, 02546 Matrix<float, 02547 Prop2, ColSymPacked, Allocator2>& B, 02548 Vector<float, VectFull, Allocator3>& w, 02549 Matrix<float, 02550 General, ColMajor, Allocator4>& z, 02551 LapackInfo& info) 02552 { 02553 int itype = 1; 02554 int n = A.GetM(); 02555 char uplo('U'); char job('V'); int lwork = 3*n; 02556 Vector<float> work(lwork); 02557 w.Reallocate(n); 02558 z.Reallocate(n,n); 02559 sspgv_(&itype, &job, &uplo, &n, A.GetData(), B.GetData(), w.GetData(), 02560 z.GetData(), &n, work.GetData(), &info.GetInfoRef()); 02561 02562 #ifdef SELDON_LAPACK_CHECK_INFO 02563 if (info.GetInfo() != 0) 02564 throw LapackError(info.GetInfo(), "GetEigenvalues", 02565 "Failed to find eigenvalues "); 02566 #endif 02567 02568 } 02569 02570 02571 template<class Prop1, class Prop2, class Allocator1, class Allocator2, 02572 class Allocator3, class Allocator4> 02573 void GetEigenvalues(Matrix<complex<float>, 02574 Prop1, ColSymPacked, Allocator1>& A, 02575 Matrix<complex<float>, 02576 Prop2, ColSymPacked, Allocator2>& B, 02577 Vector<complex<float>, VectFull, Allocator3>& alpha, 02578 Vector<complex<float>, VectFull, Allocator4>& beta, 02579 LapackInfo& info) 02580 { 02581 int n = A.GetM(); 02582 Matrix<complex<float>, General, ColMajor> C(n,n), D(n,n); 02583 for (int i = 0; i < n; i++) 02584 for (int j = 0; j < n; j++) 02585 { 02586 C(i,j) = A(i,j); 02587 D(i,j) = B(i,j); 02588 } 02589 02590 alpha.Reallocate(n); 02591 beta.Reallocate(n); 02592 GetEigenvalues(C, D, alpha, beta); 02593 } 02594 02595 02596 template<class Prop1, class Prop2, class Allocator1, class Allocator2, 02597 class Allocator3, class Allocator4, class Allocator5> 02598 void GetEigenvaluesEigenvectors(Matrix<complex<float>, 02599 Prop1, ColSymPacked, Allocator1>& A, 02600 Matrix<complex<float>, 02601 Prop2, ColSymPacked, Allocator2>& B, 02602 Vector<complex<float>, 02603 VectFull, Allocator3>& alpha, 02604 Vector<complex<float>, 02605 VectFull, Allocator4>& beta, 02606 Matrix<complex<float>, 02607 General, ColMajor, Allocator5>& z, 02608 LapackInfo& info) 02609 { 02610 int n = A.GetM(); 02611 Matrix<complex<float>, General, ColMajor> C(n,n), D(n,n); 02612 for (int i = 0; i < n; i++) 02613 for (int j = 0; j < n; j++) 02614 { 02615 C(i,j) = A(i,j); 02616 D(i,j) = B(i,j); 02617 } 02618 02619 alpha.Reallocate(n); 02620 beta.Reallocate(n); 02621 z.Reallocate(n,n); 02622 GetEigenvaluesEigenvectors(C, D, alpha, beta, z); 02623 } 02624 02625 02626 template<class Prop1, class Prop2, class Allocator1, 02627 class Allocator2, class Allocator3> 02628 void GetEigenvalues(Matrix<double, Prop1, ColSymPacked, Allocator1>& A, 02629 Matrix<double, Prop2, ColSymPacked, Allocator2>& B, 02630 Vector<double, VectFull, Allocator3>& w, 02631 LapackInfo& info) 02632 { 02633 int itype = 1; 02634 int n = A.GetM(); 02635 char uplo('U'); char job('N'); 02636 int lwork = 3*n; Vector<double> work(lwork); 02637 w.Reallocate(n); 02638 dspgv_(&itype, &job, &uplo, &n, A.GetData(), B.GetData(), w.GetData(), 02639 A.GetData(), &n, work.GetData(), &info.GetInfoRef()); 02640 02641 #ifdef SELDON_LAPACK_CHECK_INFO 02642 if (info.GetInfo() != 0) 02643 throw LapackError(info.GetInfo(), "GetEigenvalues", 02644 "Failed to find eigenvalues "); 02645 #endif 02646 02647 } 02648 02649 02650 template<class Prop1, class Prop2, class Allocator1, class Allocator2, 02651 class Allocator3, class Allocator4> 02652 void GetEigenvaluesEigenvectors(Matrix<double, 02653 Prop1, ColSymPacked, Allocator1>& A, 02654 Matrix<double, 02655 Prop2, ColSymPacked, Allocator2>& B, 02656 Vector<double, VectFull, Allocator3>& w, 02657 Matrix<double, 02658 General, ColMajor, Allocator4>& z, 02659 LapackInfo& info) 02660 { 02661 int itype = 1; 02662 int n = A.GetM(); 02663 char uplo('U'); char job('V'); int lwork = 3*n; 02664 Vector<double> work(lwork); 02665 w.Reallocate(n); 02666 z.Reallocate(n,n); 02667 dspgv_(&itype, &job, &uplo, &n, A.GetData(), B.GetData(), w.GetData(), 02668 z.GetData(), &n, work.GetData(), &info.GetInfoRef()); 02669 02670 #ifdef SELDON_LAPACK_CHECK_INFO 02671 if (info.GetInfo() != 0) 02672 throw LapackError(info.GetInfo(), "GetEigenvalues", 02673 "Failed to find eigenvalues "); 02674 #endif 02675 02676 } 02677 02678 02679 template<class Prop1, class Prop2, class Allocator1, class Allocator2, 02680 class Allocator3, class Allocator4> 02681 void GetEigenvalues(Matrix<complex<double>, 02682 Prop1, ColSymPacked, Allocator1>& A, 02683 Matrix<complex<double>, 02684 Prop2, ColSymPacked, Allocator2>& B, 02685 Vector<complex<double>, VectFull, Allocator3>& alpha, 02686 Vector<complex<double>, VectFull, Allocator4>& beta, 02687 LapackInfo& info) 02688 { 02689 int n = A.GetM(); 02690 Matrix<complex<double>, General, ColMajor> C(n,n), D(n,n); 02691 for (int i = 0; i < n; i++) 02692 for (int j = 0; j < n; j++) 02693 { 02694 C(i,j) = A(i,j); 02695 D(i,j) = B(i,j); 02696 } 02697 02698 alpha.Reallocate(n); 02699 beta.Reallocate(n); 02700 GetEigenvalues(C, D, alpha, beta); 02701 } 02702 02703 02704 template<class Prop1, class Prop2, class Allocator1, class Allocator2, 02705 class Allocator3, class Allocator4, class Allocator5> 02706 void GetEigenvaluesEigenvectors(Matrix<complex<double>, 02707 Prop1, ColSymPacked, Allocator1>& A, 02708 Matrix<complex<double>, 02709 Prop2, ColSymPacked, Allocator2>& B, 02710 Vector<complex<double>, 02711 VectFull, Allocator3>& alpha, 02712 Vector<complex<double>, 02713 VectFull, Allocator4>& beta, 02714 Matrix<complex<double>, 02715 General, ColMajor, Allocator5>& z, 02716 LapackInfo& info) 02717 { 02718 int n = A.GetM(); 02719 Matrix<complex<double>, General, ColMajor> C(n,n), D(n,n); 02720 for (int i = 0; i < n; i++) 02721 for (int j = 0; j < n; j++) 02722 { 02723 C(i,j) = A(i,j); 02724 D(i,j) = B(i,j); 02725 } 02726 02727 alpha.Reallocate(n); 02728 beta.Reallocate(n); 02729 z.Reallocate(n,n); 02730 GetEigenvaluesEigenvectors(C, D, alpha, beta, z); 02731 } 02732 02733 02734 /* RowHermPacked */ 02735 02736 02737 template<class Prop1, class Prop2, class Allocator1, 02738 class Allocator2, class Allocator3> 02739 void GetEigenvalues(Matrix<complex<float>, 02740 Prop1, RowHermPacked, Allocator1>& A, 02741 Matrix<complex<float>, 02742 Prop2, RowHermPacked, Allocator2>& B, 02743 Vector<float, VectFull, Allocator3>& w, 02744 LapackInfo& info) 02745 { 02746 int itype = 1; 02747 int n = A.GetM(); 02748 char uplo('L'); char job('N'); 02749 int lwork = 2*n; 02750 Vector<float> work(lwork); 02751 Vector<float> rwork(3*n); 02752 w.Reallocate(n); 02753 chpgv_(&itype, &job, &uplo, &n, A.GetData(), B.GetData(), w.GetData(), 02754 A.GetData(), &n, work.GetData(), rwork.GetData(), 02755 &info.GetInfoRef()); 02756 02757 #ifdef SELDON_LAPACK_CHECK_INFO 02758 if (info.GetInfo() != 0) 02759 throw LapackError(info.GetInfo(), "GetEigenvalues", 02760 "Failed to find eigenvalues "); 02761 #endif 02762 02763 } 02764 02765 02766 template<class Prop1, class Prop2, class Allocator1, class Allocator2, 02767 class Allocator3, class Allocator4> 02768 void GetEigenvaluesEigenvectors(Matrix<complex<float>, 02769 Prop1, RowHermPacked, Allocator1>& A, 02770 Matrix<complex<float>, 02771 Prop2, RowHermPacked, Allocator2>& B, 02772 Vector<float, VectFull, Allocator3>& w, 02773 Matrix<complex<float>, 02774 General, RowMajor, Allocator4>& z, 02775 LapackInfo& info) 02776 { 02777 int itype = 1; 02778 int n = A.GetM(); 02779 char uplo('L'); char job('V'); int lwork = 2*n; 02780 Vector<float> work(lwork); 02781 Vector<float> rwork(3*n); 02782 w.Reallocate(n); 02783 z.Reallocate(n,n); 02784 chpgv_(&itype, &job, &uplo, &n, A.GetData(), B.GetData(), w.GetData(), 02785 z.GetData(), &n, work.GetData(), rwork.GetData(), 02786 &info.GetInfoRef()); 02787 02788 #ifdef SELDON_LAPACK_CHECK_INFO 02789 if (info.GetInfo() != 0) 02790 throw LapackError(info.GetInfo(), "GetEigenvalues", 02791 "Failed to find eigenvalues "); 02792 #endif 02793 02794 } 02795 02796 02797 template<class Prop1, class Prop2, class Allocator1, 02798 class Allocator2, class Allocator3> 02799 void GetEigenvalues(Matrix<complex<double>, 02800 Prop1, RowHermPacked, Allocator1>& A, 02801 Matrix<complex<double>, 02802 Prop2, RowHermPacked, Allocator2>& B, 02803 Vector<double, VectFull, Allocator3>& w, 02804 LapackInfo& info) 02805 { 02806 int itype = 1; 02807 int n = A.GetM(); 02808 char uplo('L'); char job('N'); 02809 int lwork = 2*n; 02810 Vector<double> work(lwork); 02811 Vector<double> rwork(3*n); 02812 w.Reallocate(n); 02813 zhpgv_(&itype, &job, &uplo, &n, A.GetData(), B.GetData(), w.GetData(), 02814 A.GetData(), &n, work.GetData(), rwork.GetData(), 02815 &info.GetInfoRef()); 02816 02817 #ifdef SELDON_LAPACK_CHECK_INFO 02818 if (info.GetInfo() != 0) 02819 throw LapackError(info.GetInfo(), "GetEigenvalues", 02820 "Failed to find eigenvalues "); 02821 #endif 02822 02823 } 02824 02825 02826 template<class Prop1, class Prop2, class Allocator1, class Allocator2, 02827 class Allocator3, class Allocator4> 02828 void GetEigenvaluesEigenvectors(Matrix<complex<double>, 02829 Prop1, RowHermPacked, Allocator1>& A, 02830 Matrix<complex<double>, 02831 Prop2, RowHermPacked, Allocator2>& B, 02832 Vector<double, VectFull, Allocator3>& w, 02833 Matrix<complex<double>, 02834 General, RowMajor, Allocator4>& z, 02835 LapackInfo& info) 02836 { 02837 int itype = 1; 02838 int n = A.GetM(); 02839 char uplo('L'); char job('V'); int lwork = 2*n; 02840 Vector<double> work(lwork); 02841 Vector<double> rwork(3*n); 02842 w.Reallocate(n); 02843 z.Reallocate(n,n); 02844 zhpgv_(&itype, &job, &uplo, &n, A.GetData(), B.GetData(), w.GetData(), 02845 z.GetData(), &n, work.GetData(), rwork.GetData(), 02846 &info.GetInfoRef()); 02847 02848 #ifdef SELDON_LAPACK_CHECK_INFO 02849 if (info.GetInfo() != 0) 02850 throw LapackError(info.GetInfo(), "GetEigenvalues", 02851 "Failed to find eigenvalues "); 02852 #endif 02853 02854 } 02855 02856 02857 /* ColHermPacked */ 02858 02859 02860 template<class Prop1, class Prop2, class Allocator1, 02861 class Allocator2, class Allocator3> 02862 void GetEigenvalues(Matrix<complex<float>, 02863 Prop1, ColHermPacked, Allocator1>& A, 02864 Matrix<complex<float>, 02865 Prop2, ColHermPacked, Allocator2>& B, 02866 Vector<float, VectFull, Allocator3>& w, 02867 LapackInfo& info) 02868 { 02869 int itype = 1; 02870 int n = A.GetM(); 02871 char uplo('U'); char job('N'); 02872 int lwork = 2*n; 02873 Vector<float> work(lwork); 02874 Vector<float> rwork(3*n); 02875 w.Reallocate(n); 02876 chpgv_(&itype, &job, &uplo, &n, A.GetData(), B.GetData(), 02877 w.GetData(), A.GetData(), &n, work.GetData(), 02878 rwork.GetData(), &info.GetInfoRef()); 02879 02880 #ifdef SELDON_LAPACK_CHECK_INFO 02881 if (info.GetInfo() != 0) 02882 throw LapackError(info.GetInfo(), "GetEigenvalues", 02883 "Failed to find eigenvalues "); 02884 #endif 02885 02886 } 02887 02888 02889 template<class Prop1, class Prop2, class Allocator1, class Allocator2, 02890 class Allocator3, class Allocator4> 02891 void GetEigenvaluesEigenvectors(Matrix<complex<float>, 02892 Prop1, ColHermPacked, Allocator1>& A, 02893 Matrix<complex<float>, 02894 Prop2, ColHermPacked, Allocator2>& B, 02895 Vector<float, VectFull, Allocator3>& w, 02896 Matrix<complex<float>, 02897 General, ColMajor, Allocator4>& z, 02898 LapackInfo& info) 02899 { 02900 int itype = 1; 02901 int n = A.GetM(); 02902 char uplo('U'); char job('V'); 02903 int lwork = 3*n; 02904 Vector<float> work(lwork); 02905 Vector<float> rwork(3*n); 02906 w.Reallocate(n); 02907 z.Reallocate(n,n); 02908 chpgv_(&itype, &job, &uplo, &n, A.GetData(), B.GetData(), w.GetData(), 02909 z.GetData(), &n, work.GetData(), rwork.GetData(), 02910 &info.GetInfoRef()); 02911 02912 #ifdef SELDON_LAPACK_CHECK_INFO 02913 if (info.GetInfo() != 0) 02914 throw LapackError(info.GetInfo(), "GetEigenvalues", 02915 "Failed to find eigenvalues "); 02916 #endif 02917 02918 } 02919 02920 02921 template<class Prop1, class Prop2, class Allocator1, 02922 class Allocator2, class Allocator3> 02923 void GetEigenvalues(Matrix<complex<double>, 02924 Prop1, ColHermPacked, Allocator1>& A, 02925 Matrix<complex<double>, 02926 Prop2, ColHermPacked, Allocator2>& B, 02927 Vector<double, VectFull, Allocator3>& w, 02928 LapackInfo& info) 02929 { 02930 int itype = 1; 02931 int n = A.GetM(); 02932 char uplo('U'); char job('N'); 02933 int lwork = 2*n; 02934 Vector<double> work(lwork); 02935 Vector<double> rwork(3*n); 02936 w.Reallocate(n); 02937 zhpgv_(&itype, &job, &uplo, &n, A.GetData(), B.GetData(), 02938 w.GetData(), A.GetData(), &n, work.GetData(), 02939 rwork.GetData(), &info.GetInfoRef()); 02940 02941 #ifdef SELDON_LAPACK_CHECK_INFO 02942 if (info.GetInfo() != 0) 02943 throw LapackError(info.GetInfo(), "GetEigenvalues", 02944 "Failed to find eigenvalues "); 02945 #endif 02946 02947 } 02948 02949 02950 template<class Prop1, class Prop2, class Allocator1, class Allocator2, 02951 class Allocator3, class Allocator4> 02952 void GetEigenvaluesEigenvectors(Matrix<complex<double>, 02953 Prop1, ColHermPacked, Allocator1>& A, 02954 Matrix<complex<double>, 02955 Prop2, ColHermPacked, Allocator2>& B, 02956 Vector<double, VectFull, Allocator3>& w, 02957 Matrix<complex<double>, 02958 General, ColMajor, Allocator4>& z, 02959 LapackInfo& info) 02960 { 02961 int itype = 1; 02962 int n = A.GetM(); 02963 char uplo('U'); char job('V'); 02964 int lwork = 3*n; 02965 Vector<double> work(lwork); 02966 Vector<double> rwork(3*n); 02967 w.Reallocate(n); 02968 z.Reallocate(n,n); 02969 zhpgv_(&itype, &job, &uplo, &n, A.GetData(), B.GetData(), w.GetData(), 02970 z.GetData(), &n, work.GetData(), rwork.GetData(), 02971 &info.GetInfoRef()); 02972 02973 #ifdef SELDON_LAPACK_CHECK_INFO 02974 if (info.GetInfo() != 0) 02975 throw LapackError(info.GetInfo(), "GetEigenvalues", 02976 "Failed to find eigenvalues "); 02977 #endif 02978 02979 } 02980 02981 02982 /* RowMajor */ 02983 02984 02985 template<class Prop1, class Prop2, class Allocator1, 02986 class Allocator2, class Allocator3, 02987 class Allocator4, class Allocator5> 02988 void GetEigenvalues(Matrix<float, Prop1, RowMajor, Allocator1>& A, 02989 Matrix<float, Prop2, RowMajor, Allocator2>& B, 02990 Vector<float, VectFull, Allocator3>& alpha_real, 02991 Vector<float, VectFull, Allocator4>& alpha_imag, 02992 Vector<float, VectFull, Allocator5>& beta, 02993 LapackInfo& info) 02994 { 02995 int n = A.GetM(); 02996 char jobvl('N'), jobvr('N'); 02997 int lwork = 8*n+16; Vector<float> work(lwork); 02998 alpha_real.Reallocate(n); 02999 alpha_imag.Reallocate(n); 03000 beta.Reallocate(n); 03001 sggev_(&jobvl, &jobvr, &n, A.GetData(), &n, B.GetData(), &n, 03002 alpha_real.GetData(), alpha_imag.GetData(), beta.GetData(), 03003 A.GetData(), &n, A.GetData(), &n, 03004 work.GetData(), &lwork, &info.GetInfoRef()); 03005 03006 #ifdef SELDON_LAPACK_CHECK_INFO 03007 if (info.GetInfo() != 0) 03008 throw LapackError(info.GetInfo(), "GetEigenvalues", 03009 "Failed to find eigenvalues "); 03010 #endif 03011 03012 } 03013 03014 03015 template<class Prop1, class Prop2, class Prop3, class Allocator1, 03016 class Allocator2, class Allocator3, class Allocator4, 03017 class Allocator5, class Allocator6> 03018 void GetEigenvaluesEigenvectors(Matrix<float, Prop1, RowMajor, Allocator1>& A, 03019 Matrix<float, Prop2, RowMajor, Allocator2>& B, 03020 Vector<float, VectFull, Allocator3>& alpha_real, 03021 Vector<float, VectFull, Allocator4>& alpha_imag, 03022 Vector<float, VectFull, Allocator5>& beta, 03023 Matrix<float, Prop3, RowMajor, Allocator6>& V, 03024 LapackInfo& info) 03025 { 03026 int n = A.GetM(); 03027 char jobvl('V'), jobvr('N'); 03028 int lwork = 8*n+16; Vector<float> work(lwork); 03029 alpha_real.Reallocate(n); 03030 alpha_imag.Reallocate(n); 03031 beta.Reallocate(n); 03032 V.Reallocate(n,n); 03033 sggev_(&jobvl, &jobvr, &n, A.GetData(), &n, B.GetData(), &n, 03034 alpha_real.GetData(), alpha_imag.GetData(), beta.GetData(), 03035 V.GetData(), &n, V.GetData(), &n, 03036 work.GetData(), &lwork, &info.GetInfoRef()); 03037 03038 03039 #ifdef SELDON_LAPACK_CHECK_INFO 03040 if (info.GetInfo() != 0) 03041 throw LapackError(info.GetInfo(), "GetEigenvalues", 03042 "Failed to find eigenvalues "); 03043 #endif 03044 03045 Transpose(V); 03046 // conjugate if necessary 03047 int i = 0; 03048 while (i < n) 03049 { 03050 if (i < (n-1)) 03051 if (alpha_real(i) == alpha_real(i+1)) 03052 { 03053 for (int j = 0; j < n; j++) 03054 V(j,i+1) = -V(j,i+1); 03055 03056 i++; 03057 } 03058 03059 i++; 03060 } 03061 } 03062 03063 03064 template<class Prop1, class Prop2, class Allocator1, 03065 class Allocator2, class Allocator4, class Allocator5> 03066 void GetEigenvalues(Matrix<complex<float>, Prop1, RowMajor, Allocator1>& A, 03067 Matrix<complex<float>, Prop2, RowMajor, Allocator2>& B, 03068 Vector<complex<float>, VectFull, Allocator4>& alpha, 03069 Vector<complex<float>, VectFull, Allocator5>& beta, 03070 LapackInfo& info) 03071 { 03072 int n = A.GetM(); 03073 char jobvl('N'), jobvr('N'); int lwork = 2*n; 03074 Vector<complex<float> > work(lwork); 03075 Vector<float> rwork(8*n); 03076 alpha.Reallocate(n); 03077 beta.Reallocate(n); 03078 cggev_(&jobvl, &jobvr, &n, A.GetData(), &n, B.GetData(), &n, 03079 alpha.GetData(), beta.GetData(), A.GetData(), &n, A.GetData(), &n, 03080 work.GetData(), &lwork, rwork.GetData(), &info.GetInfoRef()); 03081 03082 #ifdef SELDON_LAPACK_CHECK_INFO 03083 if (info.GetInfo() != 0) 03084 throw LapackError(info.GetInfo(), "GetEigenvalues", 03085 "Failed to find eigenvalues "); 03086 #endif 03087 03088 } 03089 03090 03091 template<class Prop1, class Prop2, class Prop3, class Allocator1, 03092 class Allocator2, class Allocator4, 03093 class Allocator5, class Allocator6> 03094 void GetEigenvaluesEigenvectors(Matrix<complex<float>, 03095 Prop1, RowMajor, Allocator1>& A, 03096 Matrix<complex<float>, 03097 Prop2, RowMajor, Allocator2>& B, 03098 Vector<complex<float>, 03099 VectFull, Allocator4>& alpha, 03100 Vector<complex<float>, 03101 VectFull, Allocator5>& beta, 03102 Matrix<complex<float>, 03103 Prop3, RowMajor, Allocator6>& V, 03104 LapackInfo& info) 03105 { 03106 int n = A.GetM(); 03107 char jobvl('V'), jobvr('N'); 03108 int lwork = 2*n; Vector<complex<float> > work(lwork); 03109 Vector<float> rwork(8*n); 03110 alpha.Reallocate(n); 03111 beta.Reallocate(n); 03112 V.Reallocate(n,n); 03113 cggev_(&jobvl, &jobvr, &n, A.GetData(), &n, B.GetData(), &n, alpha.GetData(), 03114 beta.GetData(), V.GetData(), &n, V.GetData(), &n, work.GetData(), 03115 &lwork, rwork.GetData(), &info.GetInfoRef()); 03116 03117 03118 #ifdef SELDON_LAPACK_CHECK_INFO 03119 if (info.GetInfo() != 0) 03120 throw LapackError(info.GetInfo(), "GetEigenvalues", 03121 "Failed to find eigenvalues "); 03122 #endif 03123 03124 TransposeConj(V); 03125 } 03126 03127 03128 template<class Prop1, class Prop2, class Allocator1, 03129 class Allocator2, class Allocator3, 03130 class Allocator4, class Allocator5> 03131 void GetEigenvalues(Matrix<double, Prop1, RowMajor, Allocator1>& A, 03132 Matrix<double, Prop2, RowMajor, Allocator2>& B, 03133 Vector<double, VectFull, Allocator3>& alpha_real, 03134 Vector<double, VectFull, Allocator4>& alpha_imag, 03135 Vector<double, VectFull, Allocator5>& beta, 03136 LapackInfo& info) 03137 { 03138 int n = A.GetM(); 03139 char jobvl('N'), jobvr('N'); 03140 int lwork = 8*n+16; Vector<double> work(lwork); 03141 alpha_real.Reallocate(n); 03142 alpha_imag.Reallocate(n); 03143 beta.Reallocate(n); 03144 dggev_(&jobvl, &jobvr, &n, A.GetData(), &n, B.GetData(), &n, 03145 alpha_real.GetData(), alpha_imag.GetData(), beta.GetData(), 03146 A.GetData(), &n, A.GetData(), &n, 03147 work.GetData(), &lwork, &info.GetInfoRef()); 03148 03149 #ifdef SELDON_LAPACK_CHECK_INFO 03150 if (info.GetInfo() != 0) 03151 throw LapackError(info.GetInfo(), "GetEigenvalues", 03152 "Failed to find eigenvalues "); 03153 #endif 03154 03155 } 03156 03157 03158 template<class Prop1, class Prop2, class Prop3, class Allocator1, 03159 class Allocator2, class Allocator3, class Allocator4, 03160 class Allocator5, class Allocator6> 03161 void GetEigenvaluesEigenvectors(Matrix<double, Prop1, RowMajor, Allocator1>& A, 03162 Matrix<double, Prop2, RowMajor, Allocator2>& B, 03163 Vector<double, VectFull, Allocator3>& alpha_real, 03164 Vector<double, VectFull, Allocator4>& alpha_imag, 03165 Vector<double, VectFull, Allocator5>& beta, 03166 Matrix<double, Prop3, RowMajor, Allocator6>& V, 03167 LapackInfo& info) 03168 { 03169 int n = A.GetM(); 03170 char jobvl('V'), jobvr('N'); 03171 int lwork = 8*n+16; Vector<double> work(lwork); 03172 alpha_real.Reallocate(n); 03173 alpha_imag.Reallocate(n); 03174 beta.Reallocate(n); 03175 V.Reallocate(n, n); 03176 dggev_(&jobvl, &jobvr, &n, A.GetData(), &n, B.GetData(), &n, 03177 alpha_real.GetData(), alpha_imag.GetData(), beta.GetData(), 03178 V.GetData(), &n, V.GetData(), &n, 03179 work.GetData(), &lwork, &info.GetInfoRef()); 03180 03181 #ifdef SELDON_LAPACK_CHECK_INFO 03182 if (info.GetInfo() != 0) 03183 throw LapackError(info.GetInfo(), "GetEigenvalues", 03184 "Failed to find eigenvalues "); 03185 #endif 03186 03187 Transpose(V); 03188 // conjugate if necessary 03189 int i = 0; 03190 while (i < n) 03191 { 03192 if (i < (n-1)) 03193 if (alpha_real(i) == alpha_real(i+1)) 03194 { 03195 for (int j = 0; j < n; j++) 03196 V(j,i+1) = -V(j,i+1); 03197 03198 i++; 03199 } 03200 03201 i++; 03202 } 03203 } 03204 03205 03206 template<class Prop1, class Prop2, class Allocator1, 03207 class Allocator2, class Allocator4, class Allocator5> 03208 void GetEigenvalues(Matrix<complex<double>, Prop1, RowMajor, Allocator1>& A, 03209 Matrix<complex<double>, Prop2, RowMajor, Allocator2>& B, 03210 Vector<complex<double>, VectFull, Allocator4>& alpha, 03211 Vector<complex<double>, VectFull, Allocator5>& beta, 03212 LapackInfo& info) 03213 { 03214 int n = A.GetM(); 03215 char jobvl('N'), jobvr('N'); int lwork = 2*n; 03216 Vector<complex<double> > work(lwork); 03217 Vector<double> rwork(8*n); 03218 alpha.Reallocate(n); 03219 beta.Reallocate(n); 03220 zggev_(&jobvl, &jobvr, &n, A.GetData(), &n, B.GetData(), &n, 03221 alpha.GetData(), beta.GetData(), A.GetData(), &n, A.GetData(), &n, 03222 work.GetData(), &lwork, rwork.GetData(), &info.GetInfoRef()); 03223 03224 #ifdef SELDON_LAPACK_CHECK_INFO 03225 if (info.GetInfo() != 0) 03226 throw LapackError(info.GetInfo(), "GetEigenvalues", 03227 "Failed to find eigenvalues "); 03228 #endif 03229 03230 } 03231 03232 03233 template<class Prop1, class Prop2, class Prop3, class Allocator1, 03234 class Allocator2, class Allocator4, 03235 class Allocator5, class Allocator6> 03236 void GetEigenvaluesEigenvectors(Matrix<complex<double>, 03237 Prop1, RowMajor, Allocator1>& A, 03238 Matrix<complex<double>, 03239 Prop2, RowMajor, Allocator2>& B, 03240 Vector<complex<double>, 03241 VectFull, Allocator4>& alpha, 03242 Vector<complex<double>, 03243 VectFull, Allocator5>& beta, 03244 Matrix<complex<double>, 03245 Prop3, RowMajor, Allocator6>& V, 03246 LapackInfo& info) 03247 { 03248 int n = A.GetM(); 03249 char jobvl('V'), jobvr('N'); 03250 int lwork = 2*n; Vector<complex<double> > work(lwork); 03251 Vector<double> rwork(8*n); 03252 alpha.Reallocate(n); 03253 beta.Reallocate(n); 03254 V.Reallocate(n,n); 03255 zggev_(&jobvl, &jobvr, &n, A.GetData(), &n, B.GetData(), &n, alpha.GetData(), 03256 beta.GetData(), V.GetData(), &n, V.GetData(), &n, work.GetData(), 03257 &lwork, rwork.GetData(), &info.GetInfoRef()); 03258 03259 #ifdef SELDON_LAPACK_CHECK_INFO 03260 if (info.GetInfo() != 0) 03261 throw LapackError(info.GetInfo(), "GetEigenvalues", 03262 "Failed to find eigenvalues "); 03263 #endif 03264 03265 TransposeConj(V); 03266 } 03267 03268 03269 /* ColMajor */ 03270 03271 03272 template<class Prop1, class Prop2, class Allocator1, 03273 class Allocator2, class Allocator3, 03274 class Allocator4, class Allocator5> 03275 void GetEigenvalues(Matrix<float, Prop1, ColMajor, Allocator1>& A, 03276 Matrix<float, Prop2, ColMajor, Allocator2>& B, 03277 Vector<float, VectFull, Allocator3>& alpha_real, 03278 Vector<float, VectFull, Allocator4>& alpha_imag, 03279 Vector<float, VectFull, Allocator5>& beta, 03280 LapackInfo& info) 03281 { 03282 int n = A.GetM(); 03283 char jobvl('N'), jobvr('N'); 03284 int lwork = 8*n+16; Vector<float> work(lwork); 03285 alpha_real.Reallocate(n); 03286 alpha_imag.Reallocate(n); 03287 beta.Reallocate(n); 03288 03289 sggev_(&jobvl, &jobvr, &n, A.GetData(), &n, B.GetData(), &n, 03290 alpha_real.GetData(), alpha_imag.GetData(), beta.GetData(), 03291 A.GetData(), &n, A.GetData(), &n, 03292 work.GetData(), &lwork, &info.GetInfoRef()); 03293 03294 #ifdef SELDON_LAPACK_CHECK_INFO 03295 if (info.GetInfo() != 0) 03296 throw LapackError(info.GetInfo(), "GetEigenvalues", 03297 "Failed to find eigenvalues "); 03298 #endif 03299 03300 } 03301 03302 03303 template<class Prop1, class Prop2, class Prop3, class Allocator1, 03304 class Allocator2, class Allocator3, class Allocator4, 03305 class Allocator5, class Allocator6> 03306 void GetEigenvaluesEigenvectors(Matrix<float, Prop1, ColMajor, Allocator1>& A, 03307 Matrix<float, Prop2, ColMajor, Allocator2>& B, 03308 Vector<float, 03309 VectFull, Allocator3>& alpha_real, 03310 Vector<float, 03311 VectFull, Allocator4>& alpha_imag, 03312 Vector<float, VectFull, Allocator5>& beta, 03313 Matrix<float, Prop3, ColMajor, Allocator6>& V, 03314 LapackInfo& info) 03315 { 03316 int n = A.GetM(); 03317 char jobvl('V'), jobvr('N'); 03318 int lwork = 8*n+16; Vector<float> work(lwork); 03319 alpha_real.Reallocate(n); 03320 alpha_imag.Reallocate(n); 03321 beta.Reallocate(n); 03322 V.Reallocate(n,n); 03323 sggev_(&jobvl, &jobvr, &n, A.GetData(), &n, B.GetData(), 03324 &n, alpha_real.GetData(), alpha_imag.GetData(), 03325 beta.GetData(), V.GetData(), &n, V.GetData(), &n, 03326 work.GetData(), &lwork, &info.GetInfoRef()); 03327 03328 #ifdef SELDON_LAPACK_CHECK_INFO 03329 if (info.GetInfo() != 0) 03330 throw LapackError(info.GetInfo(), "GetEigenvalues", 03331 "Failed to find eigenvalues "); 03332 #endif 03333 03334 } 03335 03336 03337 template<class Prop1, class Prop2, class Allocator1, 03338 class Allocator2, class Allocator4, class Allocator5> 03339 void GetEigenvalues(Matrix<complex<float>, Prop1, ColMajor, Allocator1>& A, 03340 Matrix<complex<float>, Prop2, ColMajor, Allocator2>& B, 03341 Vector<complex<float>, VectFull, Allocator4>& alpha, 03342 Vector<complex<float>, VectFull, Allocator5>& beta, 03343 LapackInfo& info) 03344 { 03345 int n = A.GetM(); 03346 char jobvl('N'), jobvr('N'); int lwork = 2*n; 03347 Vector<complex<float> > work(lwork); 03348 Vector<float> rwork(8*n); 03349 alpha.Reallocate(n); 03350 beta.Reallocate(n); 03351 03352 cggev_(&jobvl, &jobvr, &n, A.GetData(), &n, B.GetData(), 03353 &n, alpha.GetData(), 03354 beta.GetData(), A.GetData(), &n, A.GetData(), &n, work.GetData(), 03355 &lwork, rwork.GetData(), &info.GetInfoRef()); 03356 03357 #ifdef SELDON_LAPACK_CHECK_INFO 03358 if (info.GetInfo() != 0) 03359 throw LapackError(info.GetInfo(), "GetEigenvalues", 03360 "Failed to find eigenvalues "); 03361 #endif 03362 03363 } 03364 03365 03366 template<class Prop1, class Prop2, class Prop3, class Allocator1, 03367 class Allocator2, class Allocator4, 03368 class Allocator5, class Allocator6> 03369 void GetEigenvaluesEigenvectors(Matrix<complex<float>, 03370 Prop1, ColMajor, Allocator1>& A, 03371 Matrix<complex<float>, 03372 Prop2, ColMajor, Allocator2>& B, 03373 Vector<complex<float>, 03374 VectFull, Allocator4>& alpha, 03375 Vector<complex<float>, 03376 VectFull, Allocator5>& beta, 03377 Matrix<complex<float>, 03378 Prop3, ColMajor, Allocator6>& V, 03379 LapackInfo& info) 03380 { 03381 int n = A.GetM(); 03382 char jobvl('N'), jobvr('V'); int lwork = 2*n; 03383 Vector<complex<float> > work(lwork); 03384 Vector<float> rwork(8*n); 03385 alpha.Reallocate(n); 03386 beta.Reallocate(n); 03387 V.Reallocate(n,n); 03388 cggev_(&jobvl, &jobvr, &n, A.GetData(), &n, B.GetData(), 03389 &n, alpha.GetData(), 03390 beta.GetData(), V.GetData(), &n, V.GetData(), &n, work.GetData(), 03391 &lwork, rwork.GetData(), &info.GetInfoRef()); 03392 03393 #ifdef SELDON_LAPACK_CHECK_INFO 03394 if (info.GetInfo() != 0) 03395 throw LapackError(info.GetInfo(), "GetEigenvalues", 03396 "Failed to find eigenvalues "); 03397 #endif 03398 03399 } 03400 03401 03402 template<class Prop1, class Prop2, class Allocator1, 03403 class Allocator2, class Allocator3, 03404 class Allocator4, class Allocator5> 03405 void GetEigenvalues(Matrix<double, Prop1, ColMajor, Allocator1>& A, 03406 Matrix<double, Prop2, ColMajor, Allocator2>& B, 03407 Vector<double, VectFull, Allocator3>& alpha_real, 03408 Vector<double, VectFull, Allocator4>& alpha_imag, 03409 Vector<double, VectFull, Allocator5>& beta, 03410 LapackInfo& info) 03411 { 03412 int n = A.GetM(); 03413 char jobvl('N'), jobvr('N'); 03414 int lwork = 8*n+16; Vector<double> work(lwork); 03415 alpha_real.Reallocate(n); 03416 alpha_imag.Reallocate(n); 03417 beta.Reallocate(n); 03418 03419 dggev_(&jobvl, &jobvr, &n, A.GetData(), &n, B.GetData(), &n, 03420 alpha_real.GetData(), alpha_imag.GetData(), beta.GetData(), 03421 A.GetData(), &n, A.GetData(), &n, 03422 work.GetData(), &lwork, &info.GetInfoRef()); 03423 03424 #ifdef SELDON_LAPACK_CHECK_INFO 03425 if (info.GetInfo() != 0) 03426 throw LapackError(info.GetInfo(), "GetEigenvalues", 03427 "Failed to find eigenvalues "); 03428 #endif 03429 03430 } 03431 03432 03433 template<class Prop1, class Prop2, class Prop3, class Allocator1, 03434 class Allocator2, class Allocator3, class Allocator4, 03435 class Allocator5, class Allocator6> 03436 void GetEigenvaluesEigenvectors(Matrix<double, Prop1, ColMajor, Allocator1>& A, 03437 Matrix<double, Prop2, ColMajor, Allocator2>& B, 03438 Vector<double, 03439 VectFull, Allocator3>& alpha_real, 03440 Vector<double, 03441 VectFull, Allocator4>& alpha_imag, 03442 Vector<double, VectFull, Allocator5>& beta, 03443 Matrix<double, Prop3, ColMajor, Allocator6>& V, 03444 LapackInfo& info) 03445 { 03446 int n = A.GetM(); 03447 char jobvl('V'), jobvr('N'); 03448 int lwork = 8*n+16; Vector<double> work(lwork); 03449 alpha_real.Reallocate(n); 03450 alpha_imag.Reallocate(n); 03451 beta.Reallocate(n); 03452 V.Reallocate(n,n); 03453 03454 dggev_(&jobvl, &jobvr, &n, A.GetData(), &n, B.GetData(), 03455 &n, alpha_real.GetData(), alpha_imag.GetData(), 03456 beta.GetData(), V.GetData(), &n, V.GetData(), &n, 03457 work.GetData(), &lwork, &info.GetInfoRef()); 03458 03459 #ifdef SELDON_LAPACK_CHECK_INFO 03460 if (info.GetInfo() != 0) 03461 throw LapackError(info.GetInfo(), "GetEigenvalues", 03462 "Failed to find eigenvalues "); 03463 #endif 03464 03465 } 03466 03467 03468 template<class Prop1, class Prop2, class Allocator1, 03469 class Allocator2, class Allocator4, class Allocator5> 03470 void GetEigenvalues(Matrix<complex<double>, Prop1, ColMajor, Allocator1>& A, 03471 Matrix<complex<double>, Prop2, ColMajor, Allocator2>& B, 03472 Vector<complex<double>, VectFull, Allocator4>& alpha, 03473 Vector<complex<double>, VectFull, Allocator5>& beta, 03474 LapackInfo& info) 03475 { 03476 int n = A.GetM(); 03477 char jobvl('N'), jobvr('N'); int lwork = 2*n; 03478 Vector<complex<double> > work(lwork); 03479 Vector<double> rwork(8*n); 03480 alpha.Reallocate(n); 03481 beta.Reallocate(n); 03482 03483 zggev_(&jobvl, &jobvr, &n, A.GetData(), &n, B.GetData(), 03484 &n, alpha.GetData(), 03485 beta.GetData(), A.GetData(), &n, A.GetData(), &n, work.GetData(), 03486 &lwork, rwork.GetData(), &info.GetInfoRef()); 03487 03488 #ifdef SELDON_LAPACK_CHECK_INFO 03489 if (info.GetInfo() != 0) 03490 throw LapackError(info.GetInfo(), "GetEigenvalues", 03491 "Failed to find eigenvalues "); 03492 #endif 03493 03494 } 03495 03496 03497 template<class Prop1, class Prop2, class Prop3, class Allocator1, 03498 class Allocator2, class Allocator4, 03499 class Allocator5, class Allocator6> 03500 void GetEigenvaluesEigenvectors(Matrix<complex<double>, 03501 Prop1, ColMajor, Allocator1>& A, 03502 Matrix<complex<double>, 03503 Prop2, ColMajor, Allocator2>& B, 03504 Vector<complex<double>, 03505 VectFull, Allocator4>& alpha, 03506 Vector<complex<double>, 03507 VectFull, Allocator5>& beta, 03508 Matrix<complex<double>, 03509 Prop3, ColMajor, Allocator6>& V, 03510 LapackInfo& info) 03511 { 03512 int n = A.GetM(); 03513 char jobvl('N'), jobvr('V'); int lwork = 2*n; 03514 Vector<complex<double> > work(lwork); 03515 Vector<double> rwork(8*n); 03516 alpha.Reallocate(n); 03517 beta.Reallocate(n); 03518 V.Reallocate(n,n); 03519 zggev_(&jobvl, &jobvr, &n, A.GetData(), &n, B.GetData(), 03520 &n, alpha.GetData(), 03521 beta.GetData(), V.GetData(), &n, V.GetData(), &n, work.GetData(), 03522 &lwork, rwork.GetData(), &info.GetInfoRef()); 03523 03524 #ifdef SELDON_LAPACK_CHECK_INFO 03525 if (info.GetInfo() != 0) 03526 throw LapackError(info.GetInfo(), "GetEigenvalues", 03527 "Failed to find eigenvalues "); 03528 #endif 03529 03530 } 03531 03532 03533 // GENERALIZED EIGENVALUE PROBLEM // 03535 03536 03538 // SINGULAR VALUE DECOMPOSITION // 03539 03540 03541 /* RowMajor */ 03542 03543 03544 template<class Prop1, class Allocator1, class Allocator2, 03545 class Allocator3, class Allocator4> 03546 void GetSVD(Matrix<float, Prop1, RowMajor, Allocator1>& A, 03547 Vector<float, VectFull, Allocator4>& lambda, 03548 Matrix<float, General, RowMajor, Allocator2>& u, 03549 Matrix<float, General, RowMajor, Allocator3>& v, 03550 LapackInfo& info) 03551 { 03552 int m = A.GetM(); 03553 int n = A.GetN(); 03554 char jobl('A'), jobr('A'); 03555 int lwork = max(3 * min(m, n) + max(m, n), 5 * min(m, n)); 03556 Vector<float> work(lwork); 03557 lambda.Reallocate(min(m, n)); 03558 lambda.Zero(); 03559 u.Reallocate(m, m); 03560 u.Zero(); 03561 v.Reallocate(n, n); 03562 v.Zero(); 03563 sgesvd_(&jobl, &jobr, &n, &m, A.GetData(), &n, lambda.GetData(), 03564 v.GetData(), &n, u.GetData(), &m, work.GetData(), 03565 &lwork, &info.GetInfoRef()); 03566 } 03567 03568 03569 template<class Prop1, class Allocator1, class Allocator2, 03570 class Allocator3, class Allocator4> 03571 void GetSVD(Matrix<complex<float>, Prop1, RowMajor, Allocator1>& A, 03572 Vector<float, VectFull, Allocator4>& lambda, 03573 Matrix<complex<float>, General, RowMajor, Allocator2>& u, 03574 Matrix<complex<float>, General, RowMajor, Allocator3>& v, 03575 LapackInfo& info) 03576 { 03577 int m = A.GetM(); 03578 int n = A.GetN(); 03579 char jobl('A'), jobr('A'); 03580 int lwork = 2 * min(m, n) + max(m, n); 03581 Vector<complex<float> > work(lwork); 03582 Vector<float> rwork(5 * min(m, n)); 03583 lambda.Reallocate(min(m, n)); 03584 lambda.Zero(); 03585 u.Reallocate(m, m); 03586 u.Zero(); 03587 v.Reallocate(n, n); 03588 v.Zero(); 03589 cgesvd_(&jobl, &jobr, &n, &m, A.GetDataVoid(), &n, lambda.GetData(), 03590 v.GetDataVoid(), &n, u.GetDataVoid(), &m, work.GetDataVoid(), 03591 &lwork, rwork.GetData(), &info.GetInfoRef()); 03592 } 03593 03594 03595 template<class Prop1, class Allocator1, class Allocator2, 03596 class Allocator3, class Allocator4> 03597 void GetSVD(Matrix<double, Prop1, RowMajor, Allocator1>& A, 03598 Vector<double, VectFull, Allocator4>& lambda, 03599 Matrix<double, General, RowMajor, Allocator2>& u, 03600 Matrix<double, General, RowMajor, Allocator3>& v, 03601 LapackInfo& info) 03602 { 03603 int m = A.GetM(); 03604 int n = A.GetN(); 03605 char jobl('A'), jobr('A'); 03606 int lwork = max(3 * min(m, n) + max(m, n), 5 * min(m, n)); 03607 Vector<double> work(lwork); 03608 lambda.Reallocate(min(m, n)); 03609 lambda.Zero(); 03610 u.Reallocate(m, m); 03611 u.Zero(); 03612 v.Reallocate(n, n); 03613 v.Zero(); 03614 dgesvd_(&jobl, &jobr, &n, &m, A.GetData(), &n, lambda.GetData(), 03615 v.GetData(), &n, u.GetData(), &m, work.GetData(), 03616 &lwork, &info.GetInfoRef()); 03617 } 03618 03619 03620 template<class Prop1, class Allocator1, class Allocator2, 03621 class Allocator3, class Allocator4> 03622 void GetSVD(Matrix<complex<double>, Prop1, RowMajor, Allocator1>& A, 03623 Vector<double, VectFull, Allocator4>& lambda, 03624 Matrix<complex<double>, General, RowMajor, Allocator2>& u, 03625 Matrix<complex<double>, General, RowMajor, Allocator3>& v, 03626 LapackInfo& info) 03627 { 03628 int m = A.GetM(); 03629 int n = A.GetN(); 03630 char jobl('A'), jobr('A'); 03631 int lwork = 2 * min(m, n) + max(m, n); 03632 Vector<complex<double> > work(lwork); 03633 Vector<double> rwork(5 * min(m, n)); 03634 lambda.Reallocate(min(m, n)); 03635 lambda.Zero(); 03636 u.Reallocate(m, m); 03637 u.Zero(); 03638 v.Reallocate(n, n); 03639 v.Zero(); 03640 zgesvd_(&jobl, &jobr, &n, &m, A.GetDataVoid(), &n, lambda.GetData(), 03641 v.GetDataVoid(), &n, u.GetDataVoid(), &m, work.GetDataVoid(), 03642 &lwork, rwork.GetData(), &info.GetInfoRef()); 03643 } 03644 03645 03646 /* ColMajor */ 03647 03648 03649 template<class Prop1, class Allocator1, class Allocator2, 03650 class Allocator3, class Allocator4> 03651 void GetSVD(Matrix<float, Prop1, ColMajor, Allocator1>& A, 03652 Vector<float, VectFull, Allocator4>& lambda, 03653 Matrix<float, General, ColMajor, Allocator2>& u, 03654 Matrix<float, General, ColMajor, Allocator3>& v, 03655 LapackInfo& info) 03656 { 03657 int m = A.GetM(); 03658 int n = A.GetN(); 03659 char jobl('A'), jobr('A'); 03660 int lwork = max(3 * min(m, n) + max(m, n), 5 * min(m, n)); 03661 Vector<float> work(lwork); 03662 lambda.Reallocate(min(m, n)); 03663 lambda.Zero(); 03664 u.Reallocate(m, m); 03665 u.Zero(); 03666 v.Reallocate(n, n); 03667 v.Zero(); 03668 sgesvd_(&jobl, &jobr, &m, &n, A.GetData(), &m, lambda.GetData(), 03669 u.GetData(), &m, v.GetData(), &n, work.GetData(), 03670 &lwork, &info.GetInfoRef()); 03671 } 03672 03673 03674 template<class Prop1, class Allocator1, class Allocator2, 03675 class Allocator3, class Allocator4> 03676 void GetSVD(Matrix<complex<float>, Prop1, ColMajor, Allocator1>& A, 03677 Vector<float, VectFull, Allocator4>& lambda, 03678 Matrix<complex<float>, General, ColMajor, Allocator2>& u, 03679 Matrix<complex<float>, General, ColMajor, Allocator3>& v, 03680 LapackInfo& info) 03681 { 03682 int m = A.GetM(); 03683 int n = A.GetN(); 03684 char jobl('A'), jobr('A'); 03685 int lwork = 2 * min(m, n) + max(m, n); 03686 Vector<complex<float> > work(lwork); 03687 Vector<float> rwork(5 * min(m, n)); 03688 lambda.Reallocate(min(m, n)); 03689 lambda.Zero(); 03690 u.Reallocate(m, m); 03691 u.Zero(); 03692 v.Reallocate(n, n); 03693 v.Zero(); 03694 cgesvd_(&jobl, &jobr, &m, &n, A.GetDataVoid(), &m, lambda.GetData(), 03695 u.GetDataVoid(), &m, v.GetDataVoid(), &n, work.GetDataVoid(), 03696 &lwork, rwork.GetData(), &info.GetInfoRef()); 03697 } 03698 03699 03700 template<class Prop1, class Allocator1, class Allocator2, 03701 class Allocator3, class Allocator4> 03702 void GetSVD(Matrix<double, Prop1, ColMajor, Allocator1>& A, 03703 Vector<double, VectFull, Allocator4>& lambda, 03704 Matrix<double, General, ColMajor, Allocator2>& u, 03705 Matrix<double, General, ColMajor, Allocator3>& v, 03706 LapackInfo& info) 03707 { 03708 int m = A.GetM(); 03709 int n = A.GetN(); 03710 char jobl('A'), jobr('A'); 03711 int lwork = 10 * max(m, n); 03712 Vector<double> work(lwork); 03713 lambda.Reallocate(min(m, n)); 03714 lambda.Zero(); 03715 u.Reallocate(m, m); 03716 u.Zero(); 03717 v.Reallocate(n, n); 03718 v.Zero(); 03719 dgesvd_(&jobl, &jobr, &m, &n, A.GetData(), &m, lambda.GetData(), 03720 u.GetData(), &m, v.GetData(), &n, work.GetData(), 03721 &lwork, &info.GetInfoRef()); 03722 } 03723 03724 03725 template<class Prop1, class Allocator1, class Allocator2, 03726 class Allocator3, class Allocator4> 03727 void GetSVD(Matrix<complex<double>, Prop1, ColMajor, Allocator1>& A, 03728 Vector<double, VectFull, Allocator4>& sigma, 03729 Matrix<complex<double>, General, ColMajor, Allocator2>& u, 03730 Matrix<complex<double>, General, ColMajor, Allocator3>& v, 03731 LapackInfo& info) 03732 { 03733 int m = A.GetM(); 03734 int n = A.GetN(); 03735 char jobl('A'), jobr('A'); 03736 int lwork = 2 * min(m, n) + max(m, n); 03737 Vector<complex<double> > work(lwork); 03738 Vector<double> rwork(5 * min(m, n)); 03739 sigma.Reallocate(min(m, n)); 03740 sigma.Zero(); 03741 u.Reallocate(m, m); 03742 u.Zero(); 03743 v.Reallocate(n, n); 03744 v.Zero(); 03745 zgesvd_(&jobl, &jobr, &m, &n, A.GetDataVoid(), &m, sigma.GetData(), 03746 u.GetDataVoid(), &m, v.GetDataVoid(), &n, work.GetDataVoid(), 03747 &lwork, rwork.GetData(), &info.GetInfoRef()); 03748 } 03749 03750 03751 // pseudo inverse 03752 template<class Prop1, class Allocator1> 03753 void GetPseudoInverse(Matrix<double, Prop1, ColMajor, Allocator1>& A, 03754 double epsilon, LapackInfo& info) 03755 { 03756 int m = A.GetM(), n = A.GetN(); 03757 Vector<double, VectFull, Allocator1> lambda; 03758 Matrix<double, General, ColMajor, Allocator1> U; 03759 Matrix<double, General, ColMajor, Allocator1> V; 03760 03761 GetSVD(A, lambda, U, V); 03762 03763 A.Reallocate(n, m); A.Fill(0); 03764 // computation of A = V Sigma U^* 03765 for (int k = 0; k < min(m, n); k++) 03766 if (abs(lambda(k)) > epsilon) 03767 { 03768 lambda(k) = 1.0/lambda(k); 03769 for (int i = 0; i < n; i++) 03770 for (int j = 0; j < m; j++) 03771 A(i, j) += V(k, i)*lambda(k)*U(j, k); 03772 } 03773 } 03774 03775 03776 // SINGULAR VALUE DECOMPOSITION // 03778 03779 03781 // RESOLUTION SYLVESTER EQUATION // 03782 03783 03784 void GetHessenberg(Matrix<complex<double>, General, ColMajor>& A, 03785 Matrix<complex<double>, General, ColMajor>& B, 03786 Matrix<complex<double>, General, ColMajor>& Q, 03787 Matrix<complex<double>, General, ColMajor>& Z) 03788 { 03789 char compq('V'), compz('I'); 03790 int n = A.GetM(), ilo = 1, ihi = n, info, lwork = 4 * n; 03791 Vector<complex<double> > tau(n); 03792 Vector<complex<double> > work(lwork); 03793 zgeqrf_(&n, &n, B.GetDataVoid(), &n, tau.GetDataVoid(), 03794 work.GetDataVoid(), &lwork, &info); 03795 03796 Q = B; 03797 zungqr_(&n, &n, &n, Q.GetDataVoid(), &n, tau.GetDataVoid(), 03798 work.GetDataVoid(), &lwork, &info); 03799 03800 char side('L'), trans('C'); 03801 zunmqr_(&side, &trans, &n, &n, &n, B.GetDataVoid(), &n, tau.GetDataVoid(), 03802 A.GetDataVoid(), &n, work.GetData(), &lwork, &info); 03803 03804 for (int i = 0; i < n; i++) 03805 for (int j = 0; j < i; j++) 03806 B(i, j) = 0; 03807 03808 zgghrd_(&compq, &compz, &n, &ilo, &ihi, A.GetDataVoid(), &n, 03809 B.GetDataVoid(), &n, Q.GetDataVoid(), &n, Z.GetDataVoid(), 03810 &n, &info); 03811 03812 } 03813 03814 03815 void GetQZ(Matrix<complex<double>, General, ColMajor>& A, 03816 Matrix<complex<double>, General, ColMajor>& B, 03817 Matrix<complex<double>, General, ColMajor>& Q, 03818 Matrix<complex<double>, General, ColMajor>& Z) 03819 { 03820 char compq('V'), compz('I'); 03821 int n = A.GetM(), ilo = 1, ihi = n, info, lwork = 4*n; 03822 Vector<complex<double> > tau(n); 03823 Vector<complex<double> > work(lwork); 03824 zgeqrf_(&n, &n, B.GetDataVoid(), &n, tau.GetDataVoid(), 03825 work.GetDataVoid(), &lwork, &info); 03826 03827 Q = B; 03828 zungqr_(&n, &n, &n, Q.GetDataVoid(), &n, tau.GetDataVoid(), 03829 work.GetDataVoid(), &lwork, &info); 03830 03831 char side('L'), trans('C'); 03832 zunmqr_(&side, &trans, &n, &n, &n, B.GetDataVoid(), &n, tau.GetDataVoid(), 03833 A.GetDataVoid(), &n, work.GetData(), &lwork, &info); 03834 03835 for (int i = 0; i < n; i++) 03836 for (int j = 0; j < i; j++) 03837 B(i,j) = 0; 03838 03839 zgghrd_(&compq, &compz, &n, &ilo, &ihi, A.GetDataVoid(), &n, 03840 B.GetDataVoid(), &n, Q.GetDataVoid(), &n, Z.GetDataVoid(), 03841 &n, &info); 03842 03843 char job('S'); 03844 compq = 'V'; 03845 compz = 'V'; 03846 Vector<complex<double> > alpha(n), beta(n); 03847 Vector<double> rwork(lwork); 03848 zhgeqz_(&job, &compq, &compz, &n, &ilo, &ihi, A.GetDataVoid(), &n, 03849 B.GetDataVoid(), &n, alpha.GetDataVoid(), beta.GetDataVoid(), 03850 Q.GetDataVoid(), &n, Z.GetDataVoid(), &n, work.GetDataVoid(), 03851 &lwork, rwork.GetData(), &info); 03852 } 03853 03854 03855 void GetHessenberg(Matrix<complex<double>, General, RowMajor>& A, 03856 Matrix<complex<double>, General, RowMajor>& B, 03857 Matrix<complex<double>, General, RowMajor>& Q, 03858 Matrix<complex<double>, General, RowMajor>& Z) 03859 { 03860 char compq('V'), compz('I'); 03861 int n = A.GetM(), ilo = 1, ihi = n, info, lwork = 4 * n; 03862 Transpose(A); Transpose(B); 03863 Vector<complex<double> > tau(n); 03864 Vector<complex<double> > work(lwork); 03865 zgeqrf_(&n, &n, B.GetDataVoid(), &n, tau.GetDataVoid(), 03866 work.GetDataVoid(), &lwork, &info); 03867 03868 Q = B; 03869 zungqr_(&n, &n, &n, Q.GetDataVoid(), &n, tau.GetDataVoid(), 03870 work.GetDataVoid(), &lwork, &info); 03871 03872 char side('L'), trans('C'); 03873 zunmqr_(&side, &trans, &n, &n, &n, B.GetDataVoid(), &n, tau.GetDataVoid(), 03874 A.GetDataVoid(), &n, work.GetData(), &lwork, &info); 03875 03876 for (int i = 0; i < n; i++) 03877 for (int j = 0; j < i; j++) 03878 B(j, i) = 0; 03879 03880 zgghrd_(&compq, &compz, &n, &ilo, &ihi, A.GetDataVoid(), &n, 03881 B.GetDataVoid(), &n, Q.GetDataVoid(), &n, Z.GetDataVoid(), 03882 &n, &info); 03883 03884 Transpose(A); Transpose(B); 03885 Transpose(Q); Transpose(Z); 03886 } 03887 03888 03889 void GetQZ(Matrix<complex<double>, General, RowMajor>& A, 03890 Matrix<complex<double>, General, RowMajor>& B, 03891 Matrix<complex<double>, General, RowMajor>& Q, 03892 Matrix<complex<double>, General, RowMajor>& Z) 03893 { 03894 char compq('V'), compz('I'); 03895 int n = A.GetM(), ilo = 1, ihi = n, info, lwork = 4*n; 03896 Transpose(A); Transpose(B); 03897 Vector<complex<double> > tau(n); 03898 Vector<complex<double> > work(lwork); 03899 zgeqrf_(&n, &n, B.GetDataVoid(), &n, tau.GetDataVoid(), 03900 work.GetDataVoid(), &lwork, &info); 03901 03902 Q = B; 03903 zungqr_(&n, &n, &n, Q.GetDataVoid(), &n, tau.GetDataVoid(), 03904 work.GetDataVoid(), &lwork, &info); 03905 03906 char side('L'), trans('C'); 03907 zunmqr_(&side, &trans, &n, &n, &n, B.GetDataVoid(), &n, tau.GetDataVoid(), 03908 A.GetDataVoid(), &n, work.GetData(), &lwork, &info); 03909 03910 for (int i = 0; i < n; i++) 03911 for (int j = 0; j < i; j++) 03912 B(j,i) = 0; 03913 03914 zgghrd_(&compq, &compz, &n, &ilo, &ihi, A.GetDataVoid(), &n, 03915 B.GetDataVoid(), &n, Q.GetDataVoid(), &n, Z.GetDataVoid(), 03916 &n, &info); 03917 03918 char job('S'); 03919 compq = 'V'; 03920 compz = 'V'; 03921 Vector<complex<double> > alpha(n), beta(n); 03922 Vector<double> rwork(lwork); 03923 zhgeqz_(&job, &compq, &compz, &n, &ilo, &ihi, A.GetDataVoid(), &n, 03924 B.GetDataVoid(), &n, alpha.GetDataVoid(), beta.GetDataVoid(), 03925 Q.GetDataVoid(), &n, Z.GetDataVoid(), &n, work.GetDataVoid(), 03926 &lwork, rwork.GetData(), &info); 03927 03928 Transpose(A); Transpose(B); 03929 Transpose(Q); Transpose(Z); 03930 } 03931 03932 03934 template<class T, class Prop, class Storage, class Allocator, class Vector1> 03935 void SolveHessenberg(Matrix<T, Prop, Storage, Allocator>& A, Vector1& B) 03936 { 03937 int n = A.GetM(); 03938 T tmp, pivot, invDiag; 03939 // loop over rows 03940 for (int i = 0; i < n-1; i++) 03941 { 03942 // pivoting 03943 if (abs(A(i+1, i)) > abs(A(i, i))) 03944 { 03945 // swapping rows 03946 for (int j = i; j < n; j++) 03947 { 03948 tmp = A(i, j); 03949 A(i, j) = A(i+1, j); 03950 A(i+1, j) = tmp; 03951 } 03952 03953 tmp = B(i); 03954 B(i) = B(i+1); 03955 B(i+1) = tmp; 03956 } 03957 03958 // performing elimination of A(i+1, i) 03959 invDiag = 1.0/A(i, i); 03960 pivot = A(i+1, i)*invDiag; 03961 A(i, i) = invDiag; 03962 A(i+1, i) = 0; 03963 for (int j = i+1; j < n; j++) 03964 A(i+1, j) -= pivot*A(i, j); 03965 03966 B(i+1) -= pivot*B(i); 03967 } 03968 03969 // inverting last element 03970 A(n-1, n-1) = 1.0/A(n-1, n-1); 03971 03972 // then solving triangular system 03973 for (int i = n-1; i >= 0; i--) 03974 { 03975 tmp = B(i); 03976 for (int j = i+1; j < n; j++) 03977 tmp -= A(i, j)*B(j); 03978 03979 B(i) = tmp*A(i, i); 03980 } 03981 } 03982 03983 03986 template<class T, class Prop, class Storage, class Allocator, class Vector1> 03987 void SolveHessenbergTwo(Matrix<T, Prop, Storage, Allocator>& A, Vector1& B) 03988 { 03989 int n = A.GetM(); 03990 T tmp, pivot, invDiag; 03991 T a1, a2, a3; 03992 // loop over rows 03993 for (int i = 0; i < n-2; i++) 03994 { 03995 // pivoting 03996 a1 = abs(A(i, i)); 03997 a2 = abs(A(i+1, i)); 03998 a3 = abs(A(i+2, i)); 03999 if (a2 > max(a1, a3)) 04000 { 04001 // swapping rows i and i+1 04002 for (int j = i; j < n; j++) 04003 { 04004 tmp = A(i, j); 04005 A(i, j) = A(i+1, j); 04006 A(i+1, j) = tmp; 04007 } 04008 04009 tmp = B(i); 04010 B(i) = B(i+1); 04011 B(i+1) = tmp; 04012 } 04013 else if (a3 > max(a1, a2)) 04014 { 04015 // swapping rows i+1 and i+2 04016 for (int j = i; j < n; j++) 04017 { 04018 tmp = A(i, j); 04019 A(i, j) = A(i+2, j); 04020 A(i+2, j) = tmp; 04021 } 04022 04023 tmp = B(i); 04024 B(i) = B(i+2); 04025 B(i+2) = tmp; 04026 } 04027 04028 // performing elimination of A(i+1, i) 04029 invDiag = 1.0/A(i, i); 04030 pivot = A(i+1, i)*invDiag; 04031 A(i, i) = invDiag; 04032 A(i+1, i) = 0; 04033 for (int j = i+1; j < n; j++) 04034 A(i+1, j) -= pivot*A(i, j); 04035 04036 B(i+1) -= pivot*B(i); 04037 04038 // then elimination of A(i+2, i) 04039 pivot = A(i+2, i)*invDiag; 04040 A(i+2, i) = 0; 04041 for (int j = i+1; j < n; j++) 04042 A(i+2, j) -= pivot*A(i, j); 04043 04044 B(i+2) -= pivot*B(i); 04045 } 04046 04047 // elimination of A(n, n-1) 04048 if (abs(A(n-1, n-2)) > abs(A(n-2, n-2))) 04049 { 04050 for (int j = n-2; j < n; j++) 04051 { 04052 tmp = A(n-2, j); 04053 A(n-2, j) = A(n-1, j); 04054 A(n-1, j) = tmp; 04055 } 04056 04057 tmp = B(n-2); 04058 B(n-2) = B(n-1); 04059 B(n-1) = tmp; 04060 } 04061 04062 invDiag = 1.0/A(n-2, n-2); 04063 pivot = A(n-1, n-2)*invDiag; 04064 A(n-2, n-2) = invDiag; 04065 A(n-1, n-2) = 0; 04066 A(n-1, n-1) -= pivot*A(n-2, n-1); 04067 B(n-1) -= pivot*B(n-2); 04068 04069 // inverting last element 04070 A(n-1, n-1) = 1.0/A(n-1, n-1); 04071 04072 // then solving triangular system 04073 for (int i = n-1; i >= 0; i--) 04074 { 04075 tmp = B(i); 04076 for (int j = i+1; j < n; j++) 04077 tmp -= A(i, j)*B(j); 04078 04079 B(i) = tmp*A(i, i); 04080 } 04081 } 04082 04083 04086 template<class Prop, class Storage, class Allocator> 04087 void SolveSylvester(Matrix<complex<double>, Prop, Storage, Allocator>& A, 04088 Matrix<complex<double>, Prop, Storage, Allocator>& B, 04089 Matrix<complex<double>, Prop, Storage, Allocator>& C, 04090 Matrix<complex<double>, Prop, Storage, Allocator>& D, 04091 Matrix<complex<double>, Prop, Storage, Allocator>& E) 04092 { 04093 complex<double> one(1), zero(0); 04094 int n = A.GetM(); 04095 04096 if (n <= 0) 04097 return; 04098 04099 if (n == 1) 04100 { 04101 E(0, 0) /= A(0, 0) * conj(B(0, 0)) + C(0, 0) * conj(D(0, 0)); 04102 return; 04103 } 04104 04105 Matrix<complex<double>, Prop, Storage, Allocator> Q1(n, n), Q2(n, n), 04106 Z1(n, n), Z2(n, n); 04107 Matrix<complex<double>, Prop, Storage, Allocator> Y(n, n), F(n, n); 04108 04109 GetHessenberg(A, C, Q1, Z1); 04110 GetQZ(D, B, Q2, Z2); 04111 04112 Y.Zero(); 04113 MltAdd(one, SeldonConjTrans, Q1, SeldonNoTrans, E, zero, Y); 04114 MltAdd(one, SeldonNoTrans, Y, SeldonNoTrans, Q2, zero, F); 04115 04116 // Q1 = A y_j, Q2 = C y_j 04117 Vector<complex<double> > Yvec(n); 04118 Q1.Zero(); Q2.Zero(); E.Zero(); 04119 complex<double> coef_b, coef_d; 04120 for (int k = n-1; k >= 0; k--) 04121 { 04122 // computation of Yvec = F_k 04123 // - \sum_{j=k+1}^n conj(b_{k, j}) A y_j + conj(d_{k, j}) C y_j 04124 for (int j = 0; j < n; j++) 04125 Yvec(j) = F(j, k); 04126 04127 for (int j = k+1; j < n; j++) 04128 { 04129 coef_b = conj(B(k, j)); 04130 coef_d = conj(D(k, j)); 04131 for (int i = 0; i < n; i++) 04132 Yvec(i) -= coef_b * Q1(i, j) + coef_d * Q2(i, j); 04133 } 04134 04135 coef_b = conj(B(k, k)); coef_d = conj(D(k, k)); 04136 for (int i = 0; i < n; i++) 04137 for (int j = max(0, i-1); j < n; j++) 04138 E(i, j) = coef_b * A(i, j) + coef_d * C(i, j); 04139 04140 SolveHessenberg(E, Yvec); 04141 for (int i = 0; i < n; i++) 04142 Y(i, k) = Yvec(i); 04143 04144 // computation of A y_k and C y_k 04145 for (int i = 0; i < n; i++) 04146 for (int m = max(0, i-1); m < n; m++) 04147 Q1(i, k) += A(i, m)*Yvec(m); 04148 04149 for (int i = 0; i < n; i++) 04150 for (int m = i; m < n; m++) 04151 Q2(i, k) += C(i, m)*Yvec(m); 04152 } 04153 04154 MltAdd(one, SeldonNoTrans, Y, SeldonConjTrans, Z2, zero, F); 04155 MltAdd(one, Z1, F, zero, E); 04156 } 04157 04158 04159 void GetHessenberg(Matrix<double, General, ColMajor>& A, 04160 Matrix<double, General, ColMajor>& B, 04161 Matrix<double, General, ColMajor>& Q, 04162 Matrix<double, General, ColMajor>& Z) 04163 { 04164 char compq('V'), compz('I'); 04165 int n = A.GetM(), ilo = 1, ihi = n, info, lwork = 4 * n; 04166 Vector<double> tau(n); 04167 Vector<double> work(lwork); 04168 dgeqrf_(&n, &n, B.GetData(), &n, tau.GetData(), 04169 work.GetData(), &lwork, &info); 04170 04171 Q = B; 04172 dorgqr_(&n, &n, &n, Q.GetData(), &n, tau.GetData(), 04173 work.GetData(), &lwork, &info); 04174 04175 char side('L'), trans('T'); 04176 dormqr_(&side, &trans, &n, &n, &n, B.GetData(), &n, tau.GetData(), 04177 A.GetData(), &n, work.GetData(), &lwork, &info); 04178 04179 for (int i = 0; i < n; i++) 04180 for (int j = 0; j < i; j++) 04181 B(i, j) = 0; 04182 04183 dgghrd_(&compq, &compz, &n, &ilo, &ihi, A.GetData(), &n, 04184 B.GetData(), &n, Q.GetData(), &n, Z.GetData(), 04185 &n, &info); 04186 04187 } 04188 04189 04190 void GetQZ(Matrix<double, General, ColMajor>& A, 04191 Matrix<double, General, ColMajor>& B, 04192 Matrix<double, General, ColMajor>& Q, 04193 Matrix<double, General, ColMajor>& Z) 04194 { 04195 char compq('V'), compz('I'); 04196 int n = A.GetM(), ilo = 1, ihi = n, info, lwork = 4*n; 04197 Vector<double> tau(n); 04198 Vector<double> work(lwork); 04199 dgeqrf_(&n, &n, B.GetData(), &n, tau.GetData(), 04200 work.GetData(), &lwork, &info); 04201 04202 Q = B; 04203 dorgqr_(&n, &n, &n, Q.GetData(), &n, tau.GetData(), 04204 work.GetData(), &lwork, &info); 04205 04206 char side('L'), trans('T'); 04207 dormqr_(&side, &trans, &n, &n, &n, B.GetData(), &n, tau.GetData(), 04208 A.GetData(), &n, work.GetData(), &lwork, &info); 04209 04210 for (int i = 0; i < n; i++) 04211 for (int j = 0; j < i; j++) 04212 B(i,j) = 0; 04213 04214 dgghrd_(&compq, &compz, &n, &ilo, &ihi, A.GetData(), &n, 04215 B.GetData(), &n, Q.GetData(), &n, Z.GetData(), 04216 &n, &info); 04217 04218 char job('S'); 04219 compq = 'V'; 04220 compz = 'V'; 04221 Vector<double> alphar(n), alphai(n), beta(n); 04222 Vector<double> rwork(lwork); 04223 dhgeqz_(&job, &compq, &compz, &n, &ilo, &ihi, A.GetData(), &n, 04224 B.GetData(), &n, alphar.GetData(), alphai.GetData(), beta.GetData(), 04225 Q.GetData(), &n, Z.GetData(), &n, work.GetData(), 04226 &lwork, &info); 04227 } 04228 04229 04230 void GetHessenberg(Matrix<double, General, RowMajor>& A, 04231 Matrix<double, General, RowMajor>& B, 04232 Matrix<double, General, RowMajor>& Q, 04233 Matrix<double, General, RowMajor>& Z) 04234 { 04235 char compq('V'), compz('I'); 04236 int n = A.GetM(), ilo = 1, ihi = n, info, lwork = 4 * n; 04237 Transpose(A); Transpose(B); 04238 Vector<double> tau(n); 04239 Vector<double> work(lwork); 04240 dgeqrf_(&n, &n, B.GetData(), &n, tau.GetData(), 04241 work.GetData(), &lwork, &info); 04242 04243 Q = B; 04244 dorgqr_(&n, &n, &n, Q.GetData(), &n, tau.GetData(), 04245 work.GetData(), &lwork, &info); 04246 04247 char side('L'), trans('T'); 04248 dormqr_(&side, &trans, &n, &n, &n, B.GetData(), &n, tau.GetData(), 04249 A.GetData(), &n, work.GetData(), &lwork, &info); 04250 04251 for (int i = 0; i < n; i++) 04252 for (int j = 0; j < i; j++) 04253 B(j, i) = 0; 04254 04255 dgghrd_(&compq, &compz, &n, &ilo, &ihi, A.GetData(), &n, 04256 B.GetData(), &n, Q.GetData(), &n, Z.GetData(), 04257 &n, &info); 04258 04259 Transpose(A); Transpose(B); 04260 Transpose(Q); Transpose(Z); 04261 04262 } 04263 04264 04265 void GetQZ(Matrix<double, General, RowMajor>& A, 04266 Matrix<double, General, RowMajor>& B, 04267 Matrix<double, General, RowMajor>& Q, 04268 Matrix<double, General, RowMajor>& Z) 04269 { 04270 char compq('V'), compz('I'); 04271 int n = A.GetM(), ilo = 1, ihi = n, info, lwork = 4*n; 04272 Transpose(A); Transpose(B); 04273 Vector<double> tau(n); 04274 Vector<double> work(lwork); 04275 dgeqrf_(&n, &n, B.GetData(), &n, tau.GetData(), 04276 work.GetData(), &lwork, &info); 04277 04278 Q = B; 04279 dorgqr_(&n, &n, &n, Q.GetData(), &n, tau.GetData(), 04280 work.GetData(), &lwork, &info); 04281 04282 char side('L'), trans('T'); 04283 dormqr_(&side, &trans, &n, &n, &n, B.GetData(), &n, tau.GetData(), 04284 A.GetData(), &n, work.GetData(), &lwork, &info); 04285 04286 for (int i = 0; i < n; i++) 04287 for (int j = 0; j < i; j++) 04288 B(j, i) = 0; 04289 04290 dgghrd_(&compq, &compz, &n, &ilo, &ihi, A.GetData(), &n, 04291 B.GetData(), &n, Q.GetData(), &n, Z.GetData(), 04292 &n, &info); 04293 04294 char job('S'); 04295 compq = 'V'; 04296 compz = 'V'; 04297 Vector<double> alphar(n), alphai(n), beta(n); 04298 Vector<double> rwork(lwork); 04299 dhgeqz_(&job, &compq, &compz, &n, &ilo, &ihi, A.GetData(), &n, 04300 B.GetData(), &n, alphar.GetData(), alphai.GetData(), beta.GetData(), 04301 Q.GetData(), &n, Z.GetData(), &n, work.GetData(), 04302 &lwork, &info); 04303 04304 Transpose(A); Transpose(B); 04305 Transpose(Q); Transpose(Z); 04306 } 04307 04308 04311 template<class Prop, class Storage, class Allocator> 04312 void SolveSylvester(Matrix<double, Prop, Storage, Allocator>& A, 04313 Matrix<double, Prop, Storage, Allocator>& B, 04314 Matrix<double, Prop, Storage, Allocator>& C, 04315 Matrix<double, Prop, Storage, Allocator>& D, 04316 Matrix<double, Prop, Storage, Allocator>& E) 04317 { 04318 double one(1), zero(0); 04319 int n = A.GetM(); 04320 if (n <= 0) 04321 return; 04322 04323 if (n == 1) 04324 { 04325 E(0,0) /= A(0, 0) * B(0, 0) + C(0, 0) * D(0, 0); 04326 return; 04327 } 04328 04329 Matrix<double, Prop, Storage, Allocator> Q1(n, n), Z1(n, n); 04330 Matrix<double, Prop, Storage, Allocator> Q2(n, n), Z2(n, n); 04331 Matrix<double, Prop, Storage, Allocator> Y(n, n), F(n, n); 04332 04333 GetHessenberg(A, C, Q1, Z1); 04334 GetQZ(D, B, Q2, Z2); 04335 04336 Y.Zero(); 04337 MltAdd(one, SeldonTrans, Q1, SeldonNoTrans, E, zero, Y); 04338 MltAdd(one, SeldonNoTrans, Y, SeldonNoTrans, Q2, zero, F); 04339 Y.Zero(); 04340 04341 // Q1 = A y_j, Q2 = C y_j 04342 Vector<double> Yvec(n); 04343 Q1.Zero(); Q2.Zero(); E.Zero(); 04344 double coef_b, coef_d, coef_bm1, coef_dm1, coef_b2, coef_d2, coef_d3; 04345 Vector<double> Yr(2*n); 04346 Matrix<double, Prop, Storage, Allocator> Er(2*n, 2*n); 04347 Yr.Zero(); Er.Zero(); 04348 int k = n-1; 04349 while (k >= 0) 04350 { 04351 if ((k==0) || (D(k, k-1) == 0)) 04352 { 04353 // 1x1 block : same case as for complex matrix 04354 04355 // computation of Yvec = F_k 04356 // - \sum_{j=k+1}^n b_{k, j} A y_j + d_{k, j} C y_j 04357 for (int j = 0; j < n; j++) 04358 Yvec(j) = F(j, k); 04359 04360 for (int j = k+1; j < n; j++) 04361 { 04362 coef_b = B(k, j); 04363 coef_d = D(k, j); 04364 for (int i = 0; i < n; i++) 04365 Yvec(i) -= coef_b * Q1(i, j) + coef_d * Q2(i, j); 04366 } 04367 04368 coef_b = B(k, k); coef_d = D(k, k); 04369 for (int i = 0; i < n; i++) 04370 for (int j = max(0, i-1); j < n; j++) 04371 E(i, j) = coef_b * A(i, j) + coef_d * C(i, j); 04372 04373 SolveHessenberg(E, Yvec); 04374 for (int i = 0; i < n; i++) 04375 Y(i, k) = Yvec(i); 04376 04377 // computation of A y_k and C y_k 04378 for (int i = 0; i < n; i++) 04379 for (int m = max(0, i-1); m < n; m++) 04380 Q1(i, k) += A(i, m)*Yvec(m); 04381 04382 for (int i = 0; i < n; i++) 04383 for (int m = i; m < n; m++) 04384 Q2(i, k) += C(i, m)*Yvec(m); 04385 04386 k--; 04387 } 04388 else 04389 { 04390 // 2x2 block, we have to solve 2n x 2n linear system 04391 04392 // computation of Yr = (f_k-1, f_k) 04393 for (int j = 0; j < n; j++) 04394 { 04395 Yr(2*j) = F(j, k-1); 04396 Yr(2*j+1) = F(j, k); 04397 } 04398 04399 for (int j = k+1; j < n; j++) 04400 { 04401 coef_b = B(k-1, j); 04402 coef_d = D(k-1, j); 04403 for (int i = 0; i < n; i++) 04404 Yr(2*i) -= coef_b * Q1(i, j) + coef_d * Q2(i, j); 04405 04406 coef_b = B(k, j); 04407 coef_d = D(k, j); 04408 for (int i = 0; i < n; i++) 04409 Yr(2*i+1) -= coef_b * Q1(i, j) + coef_d * Q2(i, j); 04410 } 04411 04412 // computation of linear system 04413 // b_{k-1, k-1} A + d_{k-1, k-1} C b_{k-1,k} A + d_{k-1,k} C 04414 // d_{k, k-1} C b_{k,k} A + d_{k,k} C 04415 coef_bm1 = B(k-1, k-1); coef_dm1 = D(k-1, k-1); 04416 coef_b2 = B(k-1, k); coef_d2 = D(k-1, k); coef_d3 = D(k, k-1); 04417 coef_b = B(k, k); coef_d = D(k, k); 04418 for (int i = 0; i < n; i++) 04419 for (int j = max(0, i-1); j < n; j++) 04420 { 04421 Er(2*i, 2*j) = coef_bm1 * A(i, j) + coef_dm1 * C(i, j); 04422 Er(2*i, 2*j+1) = coef_b2 * A(i, j) + coef_d2 * C(i, j); 04423 Er(2*i+1, 2*j) = coef_d3 * C(i, j); 04424 Er(2*i+1, 2*j+1) = coef_b * A(i, j) + coef_d * C(i, j); 04425 } 04426 04427 SolveHessenbergTwo(Er, Yr); 04428 04429 for (int i = 0; i < n; i++) 04430 { 04431 Y(i, k-1) = Yr(2*i); 04432 Y(i, k) = Yr(2*i+1); 04433 } 04434 04435 // computation of A y_k and C y_k 04436 for (int i = 0; i < n; i++) 04437 for (int m = max(0, i-1); m < n; m++) 04438 Q1(i, k-1) += A(i, m)*Yr(2*m); 04439 04440 for (int i = 0; i < n; i++) 04441 for (int m = i; m < n; m++) 04442 Q2(i, k-1) += C(i, m)*Yr(2*m); 04443 04444 for (int i = 0; i < n; i++) 04445 for (int m = max(0, i-1); m < n; m++) 04446 Q1(i, k) += A(i, m)*Yr(2*m+1); 04447 04448 for (int i = 0; i < n; i++) 04449 for (int m = i; m < n; m++) 04450 Q2(i, k) += C(i, m)*Yr(2*m+1); 04451 04452 k -= 2; 04453 } 04454 } 04455 04456 MltAdd(one, SeldonNoTrans, Y, SeldonTrans, Z2, zero, F); 04457 MltAdd(one, Z1, F, zero, E); 04458 } 04459 04460 04461 // RESOLUTION SYLVESTER EQUATION // 04463 04464 04465 } // end namespace 04466 04467 #define SELDON_FILE_LAPACK_EIGENVALUES_CXX 04468 #endif
