00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
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
00033
00034
00035
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
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
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
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
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
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
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
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
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
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
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
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
01575
01576
01578
01579
01580
01581
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
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
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
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
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
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
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
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
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
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
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
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
03535
03536
03538
03539
03540
03541
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
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
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
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
03778
03779
03781
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
03940 for (int i = 0; i < n-1; i++)
03941 {
03942
03943 if (abs(A(i+1, i)) > abs(A(i, i)))
03944 {
03945
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
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
03970 A(n-1, n-1) = 1.0/A(n-1, n-1);
03971
03972
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
03993 for (int i = 0; i < n-2; i++)
03994 {
03995
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
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
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
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
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
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
04070 A(n-1, n-1) = 1.0/A(n-1, n-1);
04071
04072
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
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
04123
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
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
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
04354
04355
04356
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
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
04391
04392
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
04413
04414
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
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
04463
04464
04465 }
04466
04467 #define SELDON_FILE_LAPACK_EIGENVALUES_CXX
04468 #endif
