Seldon user's guide

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