LapackDriver.cpp

Go to the documentation of this file.
00001 //Copyright (c) 2004-2005, Baris Sumengen
00002 //All rights reserved.
00003 //
00004 // CIMPL Matrix Performance Library
00005 //
00006 //Redistribution and use in source and binary
00007 //forms, with or without modification, are
00008 //permitted provided that the following
00009 //conditions are met:
00010 //
00011 //    * No commercial use is allowed. 
00012 //    This software can only be used
00013 //    for non-commercial purposes. This 
00014 //    distribution is mainly intended for
00015 //    academic research and teaching.
00016 //    * Redistributions of source code must
00017 //    retain the above copyright notice, this
00018 //    list of conditions and the following
00019 //    disclaimer.
00020 //    * Redistributions of binary form must
00021 //    mention the above copyright notice, this
00022 //    list of conditions and the following
00023 //    disclaimer in a clearly visible part 
00024 //    in associated product manual, 
00025 //    readme, and web site of the redistributed 
00026 //    software.
00027 //    * Redistributions in binary form must
00028 //    reproduce the above copyright notice,
00029 //    this list of conditions and the
00030 //    following disclaimer in the
00031 //    documentation and/or other materials
00032 //    provided with the distribution.
00033 //    * The name of Baris Sumengen may not be
00034 //    used to endorse or promote products
00035 //    derived from this software without
00036 //    specific prior written permission.
00037 //
00038 //THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT
00039 //HOLDERS AND CONTRIBUTORS "AS IS" AND ANY
00040 //EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT
00041 //NOT LIMITED TO, THE IMPLIED WARRANTIES OF
00042 //MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00043 //PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
00044 //CONTRIBUTORS BE LIABLE FOR ANY
00045 //DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00046 //EXEMPLARY, OR CONSEQUENTIAL DAMAGES
00047 //(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT
00048 //OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
00049 //DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
00050 //HOWEVER CAUSED AND ON ANY THEORY OF
00051 //LIABILITY, WHETHER IN CONTRACT, STRICT
00052 //LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
00053 //OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
00054 //OF THIS SOFTWARE, EVEN IF ADVISED OF THE
00055 //POSSIBILITY OF SUCH DAMAGE.
00056 
00057 
00058 
00059 #include "./Lapack.h"
00060 #include "mkl.h"
00061 
00062 
00063 
00064 //=========================================================
00065 // Linear equations solver
00066 //=========================================================
00067 
00068 // Gesv
00069 Matrix<float> Lapack::Gesv(Matrix<float>& A, Matrix<float>& B)
00070 {
00071         if(A.Columns() != A.Rows())
00072         {
00073                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00074                 Utility::RunTimeError("Matrix is not square!");
00075         }
00076         if(A.Columns() != B.Rows())
00077         {
00078                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00079                 Utility::RunTimeError("Matrix dimensions does not match!");
00080         }
00081         int info;
00082         Vector<int> ipiv(B.Rows());
00083         Matrix<float> LU = A.Clone();
00084         Matrix<float> X = B.Clone();
00085         int xrows = X.Rows();
00086         int xcols = X.Columns();
00087         int lurows = LU.Rows();
00088         SGESV(&xrows, &xcols, LU.Data(), &lurows, ipiv.Data(), X.Data(), &xrows, &info);
00089         if(info != 0)
00090         {
00091                 if(info < 0)
00092                 {
00093                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00094                         char s[500];
00095                         sprintf(s, "%i'th argument of SGESV is an illegal value!", -info);
00096                         Utility::RunTimeError(s);
00097                 }
00098                 else if(info > 0)  
00099                 {
00100                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00101                         char s[500];
00102                         sprintf(s, "U(%i,%i) is exactly zero. The LU factorization has been completed, but the factor U is exactly singular, so the solution could not be computed!", info, info);
00103                         Utility::RunTimeError(s);
00104                 }
00105         }
00106         return X;
00107 }
00108 
00109 Matrix<double> Lapack::Gesv(Matrix<double>& A, Matrix<double>& B)
00110 {
00111         if(A.Columns() != A.Rows())
00112         {
00113                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00114                 Utility::RunTimeError("Matrix is not square!");
00115         }
00116         if(A.Columns() != B.Rows())
00117         {
00118                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00119                 Utility::RunTimeError("Matrix dimensions does not match!");
00120         }
00121         int info;
00122         Vector<int> ipiv(B.Rows());
00123         Matrix<double> LU = A.Clone();
00124         Matrix<double> X = B.Clone();
00125         int xrows = X.Rows();
00126         int xcols = X.Columns();
00127         int lurows = LU.Rows();
00128         DGESV(&xrows, &xcols, LU.Data(), &lurows, ipiv.Data(), X.Data(), &xrows, &info);
00129         if(info != 0)
00130         {
00131                 if(info < 0)
00132                 {
00133                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00134                         char s[500];
00135                         sprintf(s, "%i'th argument of DGESV is an illegal value!", -info);
00136                         Utility::RunTimeError(s);
00137                 }
00138                 else if(info > 0)  
00139                 {
00140                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00141                         char s[500];
00142                         sprintf(s, "U(%i,%i) is exactly zero. The LU factorization has been completed, but the factor U is exactly singular, so the solution could not be computed!", info, info);
00143                         Utility::RunTimeError(s);
00144                 }
00145         }
00146         return X;
00147 }
00148 
00149 Matrix<ComplexFloat> Lapack::Gesv(Matrix<ComplexFloat>& A, Matrix<ComplexFloat>& B)
00150 {
00151         if(A.Columns() != A.Rows())
00152         {
00153                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00154                 Utility::RunTimeError("Matrix is not square!");
00155         }
00156         if(A.Columns() != B.Rows())
00157         {
00158                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00159                 Utility::RunTimeError("Matrix dimensions does not match!");
00160         }
00161         int info;
00162         Vector<int> ipiv(B.Rows());
00163         Matrix<ComplexFloat> LU = A.Clone();
00164         Matrix<ComplexFloat> X = B.Clone();
00165         int xrows = X.Rows();
00166         int xcols = X.Columns();
00167         int lurows = LU.Rows();
00168         CGESV(&xrows, &xcols, reinterpret_cast<MKL_Complex8*>(LU.Data()), &lurows, ipiv.Data(), reinterpret_cast<MKL_Complex8*>(X.Data()), &xrows, &info);
00169         if(info != 0)
00170         {
00171                 if(info < 0)
00172                 {
00173                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00174                         char s[500];
00175                         sprintf(s, "%i'th argument of CGESV is an illegal value!", -info);
00176                         Utility::RunTimeError(s);
00177                 }
00178                 else if(info > 0)  
00179                 {
00180                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00181                         char s[500];
00182                         sprintf(s, "U(%i,%i) is exactly zero. The LU factorization has been completed, but the factor U is exactly singular, so the solution could not be computed!", info, info);
00183                         Utility::RunTimeError(s);
00184                 }
00185         }
00186         return X;
00187 }
00188 
00189 Matrix<ComplexDouble> Lapack::Gesv(Matrix<ComplexDouble>& A, Matrix<ComplexDouble>& B)
00190 {
00191         if(A.Columns() != A.Rows())
00192         {
00193                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00194                 Utility::RunTimeError("Matrix is not square!");
00195         }
00196         if(A.Columns() != B.Rows())
00197         {
00198                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00199                 Utility::RunTimeError("Matrix dimensions does not match!");
00200         }
00201         int info;
00202         Vector<int> ipiv(B.Rows());
00203         Matrix<ComplexDouble> LU = A.Clone();
00204         Matrix<ComplexDouble> X = B.Clone();
00205         int xrows = X.Rows();
00206         int xcols = X.Columns();
00207         int lurows = LU.Rows();
00208         ZGESV(&xrows, &xcols, reinterpret_cast<MKL_Complex16*>(LU.Data()), &lurows, ipiv.Data(), reinterpret_cast<MKL_Complex16*>(X.Data()), &xrows, &info);
00209         if(info != 0)
00210         {
00211                 if(info < 0)
00212                 {
00213                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00214                         char s[500];
00215                         sprintf(s, "%i'th argument of ZGESV is an illegal value!", -info);
00216                         Utility::RunTimeError(s);
00217                 }
00218                 else if(info > 0)  
00219                 {
00220                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00221                         char s[500];
00222                         sprintf(s, "U(%i,%i) is exactly zero. The LU factorization has been completed, but the factor U is exactly singular, so the solution could not be computed!", info, info);
00223                         Utility::RunTimeError(s);
00224                 }
00225         }
00226         return X;
00227 }
00228 
00229 
00230 //Gesvx
00231 
00232 
00233 //Posv
00234 Matrix<float> Lapack::Posv(Matrix<float>& A, Matrix<float>& B)
00235 {
00236         if(A.Columns() != A.Rows())
00237         {
00238                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00239                 Utility::RunTimeError("Matrix is not square!");
00240         }
00241         if(A.Columns() != B.Rows())
00242         {
00243                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00244                 Utility::RunTimeError("Matrix dimensions does not match!");
00245         }
00246         int info;
00247         Matrix<float> LU = A.Clone();
00248         Matrix<float> X = B.Clone();
00249         int xrows = X.Rows();
00250         int xcols = X.Columns();
00251         int lurows = LU.Rows();
00252         SPOSV("U", &xrows, &xcols, LU.Data(), &lurows, X.Data(), &xrows, &info);
00253         if(info != 0)
00254         {
00255                 if(info < 0)
00256                 {
00257                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00258                         char s[500];
00259                         sprintf(s, "%i'th argument of SPOSV is an illegal value!", -info);
00260                         Utility::RunTimeError(s);
00261                 }
00262                 else if(info > 0)  
00263                 {
00264                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00265                         char s[500];
00266                         sprintf(s, "The leading minor of order %i (and hence the matrix A itself) is not positive definite, so the factorization could not be completed, and the solution has not been computed!", info);
00267                         Utility::RunTimeError(s);
00268                 }
00269         }
00270         return X;
00271 }
00272 
00273 Matrix<double> Lapack::Posv(Matrix<double>& A, Matrix<double>& B)
00274 {
00275         if(A.Columns() != A.Rows())
00276         {
00277                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00278                 Utility::RunTimeError("Matrix is not square!");
00279         }
00280         if(A.Columns() != B.Rows())
00281         {
00282                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00283                 Utility::RunTimeError("Matrix dimensions does not match!");
00284         }
00285         int info;
00286         Matrix<double> LU = A.Clone();
00287         Matrix<double> X = B.Clone();
00288         int xrows = X.Rows();
00289         int xcols = X.Columns();
00290         int lurows = LU.Rows();
00291         DPOSV("U", &xrows, &xcols, LU.Data(), &lurows, X.Data(), &xrows, &info);
00292         if(info != 0)
00293         {
00294                 if(info < 0)
00295                 {
00296                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00297                         char s[500];
00298                         sprintf(s, "%i'th argument of DPOSV is an illegal value!", -info);
00299                         Utility::RunTimeError(s);
00300                 }
00301                 else if(info > 0)  
00302                 {
00303                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00304                         char s[500];
00305                         sprintf(s, "The leading minor of order %i (and hence the matrix A itself) is not positive definite, so the factorization could not be completed, and the solution has not been computed!", info);
00306                         Utility::RunTimeError(s);
00307                 }
00308         }
00309         return X;
00310 }
00311 
00312 Matrix<ComplexFloat> Lapack::Posv(Matrix<ComplexFloat>& A, Matrix<ComplexFloat>& B)
00313 {
00314         if(A.Columns() != A.Rows())
00315         {
00316                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00317                 Utility::RunTimeError("Matrix is not square!");
00318         }
00319         if(A.Columns() != B.Rows())
00320         {
00321                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00322                 Utility::RunTimeError("Matrix dimensions does not match!");
00323         }
00324         int info;
00325         Matrix<ComplexFloat> LU = A.Clone();
00326         Matrix<ComplexFloat> X = B.Clone();
00327         int xrows = X.Rows();
00328         int xcols = X.Columns();
00329         int lurows = LU.Rows();
00330         CPOSV("U", &xrows, &xcols, reinterpret_cast<MKL_Complex8*>(LU.Data()), &lurows, reinterpret_cast<MKL_Complex8*>(X.Data()), &xrows, &info);
00331         if(info != 0)
00332         {
00333                 if(info < 0)
00334                 {
00335                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00336                         char s[500];
00337                         sprintf(s, "%i'th argument of CPOSV is an illegal value!", -info);
00338                         Utility::RunTimeError(s);
00339                 }
00340                 else if(info > 0)  
00341                 {
00342                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00343                         char s[500];
00344                         sprintf(s, "The leading minor of order %i (and hence the matrix A itself) is not positive definite, so the factorization could not be completed, and the solution has not been computed!", info);
00345                         Utility::RunTimeError(s);
00346                 }
00347         }
00348         return X;
00349 }
00350 
00351 Matrix<ComplexDouble> Lapack::Posv(Matrix<ComplexDouble>& A, Matrix<ComplexDouble>& B)
00352 {
00353         if(A.Columns() != A.Rows())
00354         {
00355                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00356                 Utility::RunTimeError("Matrix is not square!");
00357         }
00358         if(A.Columns() != B.Rows())
00359         {
00360                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00361                 Utility::RunTimeError("Matrix dimensions does not match!");
00362         }
00363         int info;
00364         Matrix<ComplexDouble> LU = A.Clone();
00365         Matrix<ComplexDouble> X = B.Clone();
00366         int xrows = X.Rows();
00367         int xcols = X.Columns();
00368         int lurows = LU.Rows();
00369         ZPOSV("U", &xrows, &xcols, reinterpret_cast<MKL_Complex16*>(LU.Data()), &lurows, reinterpret_cast<MKL_Complex16*>(X.Data()), &xrows, &info);
00370         if(info != 0)
00371         {
00372                 if(info < 0)
00373                 {
00374                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00375                         char s[500];
00376                         sprintf(s, "%i'th argument of ZPOSV is an illegal value!", -info);
00377                         Utility::RunTimeError(s);
00378                 }
00379                 else if(info > 0)  
00380                 {
00381                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00382                         char s[500];
00383                         sprintf(s, "The leading minor of order %i (and hence the matrix A itself) is not positive definite, so the factorization could not be completed, and the solution has not been computed!", info);
00384                         Utility::RunTimeError(s);
00385                 }
00386         }
00387         return X;
00388 }
00389 
00390 //Posvx
00391 
00392 
00393 //Sysv
00394 Matrix<float> Lapack::Sysv(Matrix<float>& A, Matrix<float>& B)
00395 {
00396         if(A.Columns() != A.Rows())
00397         {
00398                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00399                 Utility::RunTimeError("Matrix is not square!");
00400         }
00401         if(A.Columns() != B.Rows())
00402         {
00403                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00404                 Utility::RunTimeError("Matrix dimensions does not match!");
00405         }
00406         int info;
00407         Vector<int> ipiv(B.Rows());
00408         Matrix<float> LU = A.Clone();
00409         Matrix<float> X = B.Clone();
00410         int xrows = X.Rows();
00411         int xcols = X.Columns();
00412         int lurows = LU.Rows();
00413         int bs = 64*lurows;
00414         Vector<float> work(bs);
00415         SSYSV("U", &xrows, &xcols, LU.Data(), &lurows, ipiv.Data(), X.Data(), &xrows, work.Data(), &bs, &info);
00416 
00417         if(info != 0)
00418         {
00419                 if(info < 0)
00420                 {
00421                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00422                         char s[500];
00423                         sprintf(s, "%i'th argument of SSYSV is an illegal value!", -info);
00424                         Utility::RunTimeError(s);
00425                 }
00426                 else if(info > 0)  
00427                 {
00428                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00429                         char s[500];
00430                         sprintf(s, "d(%i,%i) is 0. The factorization has been completed, but D is exactly singular, so the solution could not be computed!", info, info);
00431                         Utility::RunTimeError(s);
00432                 }
00433         }
00434         return X;
00435 }
00436 
00437 Matrix<double> Lapack::Sysv(Matrix<double>& A, Matrix<double>& B)
00438 {
00439         if(A.Columns() != A.Rows())
00440         {
00441                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00442                 Utility::RunTimeError("Matrix is not square!");
00443         }
00444         if(A.Columns() != B.Rows())
00445         {
00446                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00447                 Utility::RunTimeError("Matrix dimensions does not match!");
00448         }
00449         int info;
00450         Vector<int> ipiv(B.Rows());
00451         Matrix<double> LU = A.Clone();
00452         Matrix<double> X = B.Clone();
00453         int xrows = X.Rows();
00454         int xcols = X.Columns();
00455         int lurows = LU.Rows();
00456         int bs = 64*lurows;
00457         Vector<double> work(bs);
00458         DSYSV("U", &xrows, &xcols, LU.Data(), &lurows, ipiv.Data(), X.Data(), &xrows, work.Data(), &bs, &info);
00459 
00460         if(info != 0)
00461         {
00462                 if(info < 0)
00463                 {
00464                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00465                         char s[500];
00466                         sprintf(s, "%i'th argument of DSYSV is an illegal value!", -info);
00467                         Utility::RunTimeError(s);
00468                 }
00469                 else if(info > 0)  
00470                 {
00471                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00472                         char s[500];
00473                         sprintf(s, "d(%i,%i) is 0. The factorization has been completed, but D is exactly singular, so the solution could not be computed!", info, info);
00474                         Utility::RunTimeError(s);
00475                 }
00476         }
00477         return X;
00478 }
00479 
00480 Matrix<ComplexFloat> Lapack::Sysv(Matrix<ComplexFloat>& A, Matrix<ComplexFloat>& B)
00481 {
00482         if(A.Columns() != A.Rows())
00483         {
00484                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00485                 Utility::RunTimeError("Matrix is not square!");
00486         }
00487         if(A.Columns() != B.Rows())
00488         {
00489                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00490                 Utility::RunTimeError("Matrix dimensions does not match!");
00491         }
00492         int info;
00493         Vector<int> ipiv(B.Rows());
00494         Matrix<ComplexFloat> LU = A.Clone();
00495         Matrix<ComplexFloat> X = B.Clone();
00496         int xrows = X.Rows();
00497         int xcols = X.Columns();
00498         int lurows = LU.Rows();
00499         int bs = 64*lurows;
00500         Vector<ComplexFloat> work(bs);
00501         CSYSV("U", &xrows, &xcols, reinterpret_cast<MKL_Complex8*>(LU.Data()), &lurows, ipiv.Data(), reinterpret_cast<MKL_Complex8*>(X.Data()), &xrows, reinterpret_cast<MKL_Complex8*>(work.Data()), &bs, &info);
00502 
00503         if(info != 0)
00504         {
00505                 if(info < 0)
00506                 {
00507                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00508                         char s[500];
00509                         sprintf(s, "%i'th argument of CSYSV is an illegal value!", -info);
00510                         Utility::RunTimeError(s);
00511                 }
00512                 else if(info > 0)  
00513                 {
00514                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00515                         char s[500];
00516                         sprintf(s, "d(%i,%i) is 0. The factorization has been completed, but D is exactly singular, so the solution could not be computed!", info, info);
00517                         Utility::RunTimeError(s);
00518                 }
00519         }
00520         return X;
00521 }
00522 
00523 Matrix<ComplexDouble> Lapack::Sysv(Matrix<ComplexDouble>& A, Matrix<ComplexDouble>& B)
00524 {
00525         if(A.Columns() != A.Rows())
00526         {
00527                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00528                 Utility::RunTimeError("Matrix is not square!");
00529         }
00530         if(A.Columns() != B.Rows())
00531         {
00532                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00533                 Utility::RunTimeError("Matrix dimensions does not match!");
00534         }
00535         int info;
00536         Vector<int> ipiv(B.Rows());
00537         Matrix<ComplexDouble> LU = A.Clone();
00538         Matrix<ComplexDouble> X = B.Clone();
00539         int xrows = X.Rows();
00540         int xcols = X.Columns();
00541         int lurows = LU.Rows();
00542         int bs = 64*lurows;
00543         Vector<ComplexDouble> work(bs);
00544         ZSYSV("U", &xrows, &xcols, reinterpret_cast<MKL_Complex16*>(LU.Data()), &lurows, ipiv.Data(), reinterpret_cast<MKL_Complex16*>(X.Data()), &xrows, reinterpret_cast<MKL_Complex16*>(work.Data()), &bs, &info);
00545 
00546         if(info != 0)
00547         {
00548                 if(info < 0)
00549                 {
00550                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00551                         char s[500];
00552                         sprintf(s, "%i'th argument of ZSYSV is an illegal value!", -info);
00553                         Utility::RunTimeError(s);
00554                 }
00555                 else if(info > 0)  
00556                 {
00557                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00558                         char s[500];
00559                         sprintf(s, "d(%i,%i) is 0. The factorization has been completed, but D is exactly singular, so the solution could not be computed!", info, info);
00560                         Utility::RunTimeError(s);
00561                 }
00562         }
00563         return X;
00564 }
00565 
00566 //Sysvx
00567 
00568 
00569 
00570 
00571 //Hesv
00572 Matrix<ComplexFloat> Lapack::Hesv(Matrix<ComplexFloat>& A, Matrix<ComplexFloat>& B)
00573 {
00574         if(A.Columns() != A.Rows())
00575         {
00576                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00577                 Utility::RunTimeError("Matrix is not square!");
00578         }
00579         if(A.Columns() != B.Rows())
00580         {
00581                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00582                 Utility::RunTimeError("Matrix dimensions does not match!");
00583         }
00584         int info;
00585         Vector<int> ipiv(B.Rows());
00586         Matrix<ComplexFloat> LU = A.Clone();
00587         Matrix<ComplexFloat> X = B.Clone();
00588         int xrows = X.Rows();
00589         int xcols = X.Columns();
00590         int lurows = LU.Rows();
00591         int bs = 64*lurows;
00592         Vector<ComplexFloat> work(bs);
00593         CHESV("U", &xrows, &xcols, reinterpret_cast<MKL_Complex8*>(LU.Data()), &lurows, ipiv.Data(), reinterpret_cast<MKL_Complex8*>(X.Data()), &xrows, reinterpret_cast<MKL_Complex8*>(work.Data()), &bs, &info);
00594 
00595         if(info != 0)
00596         {
00597                 if(info < 0)
00598                 {
00599                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00600                         char s[500];
00601                         sprintf(s, "%i'th argument of CHESV is an illegal value!", -info);
00602                         Utility::RunTimeError(s);
00603                 }
00604                 else if(info > 0)  
00605                 {
00606                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00607                         char s[500];
00608                         sprintf(s, "d(%i,%i) is 0. The factorization has been completed, but D is exactly singular, so the solution could not be computed!", info, info);
00609                         Utility::RunTimeError(s);
00610                 }
00611         }
00612         return X;
00613 }
00614 
00615 Matrix<ComplexDouble> Lapack::Hesv(Matrix<ComplexDouble>& A, Matrix<ComplexDouble>& B)
00616 {
00617         if(A.Columns() != A.Rows())
00618         {
00619                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00620                 Utility::RunTimeError("Matrix is not square!");
00621         }
00622         if(A.Columns() != B.Rows())
00623         {
00624                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00625                 Utility::RunTimeError("Matrix dimensions does not match!");
00626         }
00627         int info;
00628         Vector<int> ipiv(B.Rows());
00629         Matrix<ComplexDouble> LU = A.Clone();
00630         Matrix<ComplexDouble> X = B.Clone();
00631         int xrows = X.Rows();
00632         int xcols = X.Columns();
00633         int lurows = LU.Rows();
00634         int bs = 64*lurows;
00635         Vector<ComplexDouble> work(bs);
00636         ZHESV("U", &xrows, &xcols, reinterpret_cast<MKL_Complex16*>(LU.Data()), &lurows, ipiv.Data(), reinterpret_cast<MKL_Complex16*>(X.Data()), &xrows, reinterpret_cast<MKL_Complex16*>(work.Data()), &bs, &info);
00637 
00638         if(info != 0)
00639         {
00640                 if(info < 0)
00641                 {
00642                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00643                         char s[500];
00644                         sprintf(s, "%i'th argument of ZHESV is an illegal value!", -info);
00645                         Utility::RunTimeError(s);
00646                 }
00647                 else if(info > 0)  
00648                 {
00649                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00650                         char s[500];
00651                         sprintf(s, "d(%i,%i) is 0. The factorization has been completed, but D is exactly singular, so the solution could not be computed!", info, info);
00652                         Utility::RunTimeError(s);
00653                 }
00654         }
00655         return X;
00656 }
00657 
00658 
00659 //Hesvx
00660 
00661 
00662 
00663 
00664 
00665 
00666 
00667 
00668 
00669 
00670 
00671 
00672 
00673 
00674 
00675 //=========================================================
00676 // Linear Least Squares problems
00677 //=========================================================
00678 
00679 
00680 
00681 // Gels
00682 Matrix<float> Lapack::Gels(Matrix<float>& A, Matrix<float>& B)
00683 {
00684         if(A.Rows() != B.Rows())
00685         {
00686                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00687                 Utility::RunTimeError("Matrix dimensions does not match!");
00688         }
00689 
00690         int info;
00691         Matrix<float> F = A.Clone();
00692         Matrix<float> X;
00693         if(B.Rows() >= A.Columns())
00694         {
00695                 X = B.Clone();
00696         }
00697         else
00698         {
00699                 X = Matrix<float>(A.Columns(), B.Columns());
00700                 X.ReadFromMatrix(B);
00701         }
00702         int m = F.Rows();
00703         int nrhs = X.Columns();
00704         int n = F.Columns();
00705         int ldb = X.Rows();
00706 
00707         int min = (m>n)?n:m;
00708         int max = (m>n)?m:n;
00709         max = (max>nrhs)?max:nrhs;
00710         int bs = min+64*max;
00711         Vector<float> work(bs);
00712         SGELS("N", &m, &n, &nrhs, F.Data(), &m, X.Data(), &ldb, work.Data(), &bs, &info);
00713         if(info != 0)
00714         {
00715                 if(info < 0)
00716                 {
00717                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00718                         char s[500];
00719                         sprintf(s, "%i'th argument of SGELS is an illegal value!", -info);
00720                         Utility::RunTimeError(s);
00721                 }
00722                 else if(info > 0)  
00723                 {
00724                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00725                         Utility::RunTimeError("Something went wrong when running SGELS!");
00726                 }
00727         }
00728         if(B.Rows() > A.Columns())
00729         {
00730                 X = X(0, A.Columns()-1, 0, nrhs-1);
00731         }
00732         return X;
00733 }
00734 
00735 Matrix<double> Lapack::Gels(Matrix<double>& A, Matrix<double>& B)
00736 {
00737         if(A.Rows() != B.Rows())
00738         {
00739                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00740                 Utility::RunTimeError("Matrix dimensions does not match!");
00741         }
00742 
00743         int info;
00744         Matrix<double> F = A.Clone();
00745         Matrix<double> X;
00746         if(B.Rows() >= A.Columns())
00747         {
00748                 X = B.Clone();
00749         }
00750         else
00751         {
00752                 X = Matrix<double>(A.Columns(), B.Columns());
00753                 X.ReadFromMatrix(B);
00754         }
00755         int m = F.Rows();
00756         int nrhs = X.Columns();
00757         int n = F.Columns();
00758         int ldb = X.Rows();
00759 
00760         int min = (m>n)?n:m;
00761         int max = (m>n)?m:n;
00762         max = (max>nrhs)?max:nrhs;
00763         int bs = min+64*max;
00764         Vector<double> work(bs);
00765         DGELS("N", &m, &n, &nrhs, F.Data(), &m, X.Data(), &ldb, work.Data(), &bs, &info);
00766         if(info != 0)
00767         {
00768                 if(info < 0)
00769                 {
00770                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00771                         char s[500];
00772                         sprintf(s, "%i'th argument of DGELS is an illegal value!", -info);
00773                         Utility::RunTimeError(s);
00774                 }
00775                 else if(info > 0)  
00776                 {
00777                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00778                         Utility::RunTimeError("Something went wrong when running DGELS!");
00779                 }
00780         }
00781         if(B.Rows() > A.Columns())
00782         {
00783                 X = X(0, A.Columns()-1, 0, nrhs-1);
00784         }
00785         return X;
00786 }
00787 
00788 Matrix<ComplexFloat> Lapack::Gels(Matrix<ComplexFloat>& A, Matrix<ComplexFloat>& B)
00789 {
00790         if(A.Rows() != B.Rows())
00791         {
00792                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00793                 Utility::RunTimeError("Matrix dimensions does not match!");
00794         }
00795 
00796         int info;
00797         Matrix<ComplexFloat> F = A.Clone();
00798         Matrix<ComplexFloat> X;
00799         if(B.Rows() >= A.Columns())
00800         {
00801                 X = B.Clone();
00802         }
00803         else
00804         {
00805                 X = Matrix<ComplexFloat>(A.Columns(), B.Columns());
00806                 X.ReadFromMatrix(B);
00807         }
00808         int m = F.Rows();
00809         int nrhs = X.Columns();
00810         int n = F.Columns();
00811         int ldb = X.Rows();
00812 
00813         int min = (m>n)?n:m;
00814         int max = (m>n)?m:n;
00815         max = (max>nrhs)?max:nrhs;
00816         int bs = min+64*max;
00817         Vector<ComplexFloat> work(bs);
00818         CGELS("N", &m, &n, &nrhs, reinterpret_cast<MKL_Complex8*>(F.Data()), &m, reinterpret_cast<MKL_Complex8*>(X.Data()), &ldb, reinterpret_cast<MKL_Complex8*>(work.Data()), &bs, &info);
00819         if(info != 0)
00820         {
00821                 if(info < 0)
00822                 {
00823                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00824                         char s[500];
00825                         sprintf(s, "%i'th argument of CGELS is an illegal value!", -info);
00826                         Utility::RunTimeError(s);
00827                 }
00828                 else if(info > 0)  
00829                 {
00830                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00831                         Utility::RunTimeError("Something went wrong when running CGELS!");
00832                 }
00833         }
00834         if(B.Rows() > A.Columns())
00835         {
00836                 X = X(0, A.Columns()-1, 0, nrhs-1);
00837         }
00838         return X;
00839 }
00840 
00841 Matrix<ComplexDouble> Lapack::Gels(Matrix<ComplexDouble>& A, Matrix<ComplexDouble>& B)
00842 {
00843         if(A.Rows() != B.Rows())
00844         {
00845                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00846                 Utility::RunTimeError("Matrix dimensions does not match!");
00847         }
00848 
00849         int info;
00850         Matrix<ComplexDouble> F = A.Clone();
00851         Matrix<ComplexDouble> X;
00852         if(B.Rows() >= A.Columns())
00853         {
00854                 X = B.Clone();
00855         }
00856         else
00857         {
00858                 X = Matrix<ComplexDouble>(A.Columns(), B.Columns());
00859                 X.ReadFromMatrix(B);
00860         }
00861         int m = F.Rows();
00862         int nrhs = X.Columns();
00863         int n = F.Columns();
00864         int ldb = X.Rows();
00865 
00866         int min = (m>n)?n:m;
00867         int max = (m>n)?m:n;
00868         max = (max>nrhs)?max:nrhs;
00869         int bs = min+64*max;
00870         Vector<ComplexDouble> work(bs);
00871         ZGELS("N", &m, &n, &nrhs, reinterpret_cast<MKL_Complex16*>(F.Data()), &m, reinterpret_cast<MKL_Complex16*>(X.Data()), &ldb, reinterpret_cast<MKL_Complex16*>(work.Data()), &bs, &info);
00872         if(info != 0)
00873         {
00874                 if(info < 0)
00875                 {
00876                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00877                         char s[500];
00878                         sprintf(s, "%i'th argument of ZGELS is an illegal value!", -info);
00879                         Utility::RunTimeError(s);
00880                 }
00881                 else if(info > 0)  
00882                 {
00883                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00884                         Utility::RunTimeError("Something went wrong when running ZGELS!");
00885                 }
00886         }
00887 
00888         if(B.Rows() > A.Columns())
00889         {
00890                 X = X(0, A.Columns()-1, 0, nrhs-1);
00891         }
00892 
00893         return X;
00894 }
00895 
00896 
00897 // Gelsy
00898 Matrix<float> Lapack::Gelsy(Matrix<float>& A, Matrix<float>& B)
00899 {
00900         if(A.Rows() != B.Rows())
00901         {
00902                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00903                 Utility::RunTimeError("Matrix dimensions does not match!");
00904         }
00905 
00906         int info;
00907         int rank;
00908         Matrix<float> F = A.Clone();
00909         Matrix<float> X;
00910         if(B.Rows() >= A.Columns())
00911         {
00912                 X = B.Clone();
00913         }
00914         else
00915         {
00916                 X = Matrix<float>(A.Columns(), B.Columns());
00917                 X.ReadFromMatrix(B);
00918         }
00919         int m = F.Rows();
00920         int nrhs = X.Columns();
00921         int n = F.Columns();
00922         int ldb = X.Rows();
00923 
00924         int min = (m>n)?n:m;
00925         int max = (m>n)?m:n;
00926         int bs = min+2*n+64*(n+1);
00927         int temp = 2*min+64*nrhs;
00928         bs = (bs>temp)?bs:temp;
00929         Vector<float> work(bs);
00930         Vector<int> jpvt(n);
00931         float rcond = max*numeric_limits<float>::epsilon();
00932 
00933         SGELSY(&m, &n, &nrhs, F.Data(), &m, X.Data(), &ldb, jpvt.Data(), &rcond, &rank, work.Data(), &bs, &info);
00934         if(rank < min && info == 0)
00935         {
00936                         char s[500];
00937                         sprintf(s, "Matrix is rank deficient (rank is %i)!", rank);
00938                         Utility::Warning(s);
00939         }
00940         if(info != 0)
00941         {
00942                 if(info < 0)
00943                 {
00944                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00945                         char s[500];
00946                         sprintf(s, "%i'th argument of SGELSY is an illegal value!", -info);
00947                         Utility::RunTimeError(s);
00948                 }
00949                 else if(info > 0)  
00950                 {
00951                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00952                         Utility::RunTimeError("Something went wrong when running SGELSY!");
00953                 }
00954         }
00955         if(B.Rows() > A.Columns())
00956         {
00957                 X = X(0, A.Columns()-1, 0, nrhs-1);
00958         }
00959         return X;
00960 }
00961 
00962 
00963 
00964 Matrix<double> Lapack::Gelsy(Matrix<double>& A, Matrix<double>& B)
00965 {
00966         if(A.Rows() != B.Rows())
00967         {
00968                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00969                 Utility::RunTimeError("Matrix dimensions does not match!");
00970         }
00971 
00972         int info;
00973         int rank;
00974         Matrix<double> F = A.Clone();
00975         Matrix<double> X;
00976         if(B.Rows() >= A.Columns())
00977         {
00978                 X = B.Clone();
00979         }
00980         else
00981         {
00982                 X = Matrix<double>(A.Columns(), B.Columns());
00983                 X.ReadFromMatrix(B);
00984         }
00985         int m = F.Rows();
00986         int nrhs = X.Columns();
00987         int n = F.Columns();
00988         int ldb = X.Rows();
00989 
00990         int min = (m>n)?n:m;
00991         int max = (m>n)?m:n;
00992         int bs = min+2*n+64*(n+1);
00993         int temp = 2*min+64*nrhs;
00994         bs = (bs>temp)?bs:temp;
00995         Vector<double> work(bs);
00996         Vector<int> jpvt(n);
00997         double rcond = max*numeric_limits<double>::epsilon();
00998 
00999         DGELSY(&m, &n, &nrhs, F.Data(), &m, X.Data(), &ldb, jpvt.Data(), &rcond, &rank, work.Data(), &bs, &info);
01000         if(rank < min && info == 0)
01001         {
01002                         char s[500];
01003                         sprintf(s, "Matrix is rank deficient (rank is %i)!", rank);
01004                         Utility::Warning(s);
01005         }
01006         if(info != 0)
01007         {
01008                 if(info < 0)
01009                 {
01010                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
01011                         char s[500];
01012                         sprintf(s, "%i'th argument of DGELSY is an illegal value!", -info);
01013                         Utility::RunTimeError(s);
01014                 }
01015                 else if(info > 0)  
01016                 {
01017                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
01018                         Utility::RunTimeError("Something went wrong when running DGELSY!");
01019                 }
01020         }
01021         if(B.Rows() > A.Columns())
01022         {
01023                 X = X(0, A.Columns()-1, 0, nrhs-1);
01024         }
01025         return X;
01026 }
01027 
01028 Matrix<ComplexFloat> Lapack::Gelsy(Matrix<ComplexFloat>& A, Matrix<ComplexFloat>& B)
01029 {
01030         if(A.Rows() != B.Rows())
01031         {
01032                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
01033                 Utility::RunTimeError("Matrix dimensions does not match!");
01034         }
01035 
01036         int info;
01037         int rank;
01038         Matrix<ComplexFloat> F = A.Clone();
01039         Matrix<ComplexFloat> X;
01040         if(B.Rows() >= A.Columns())
01041         {
01042                 X = B.Clone();
01043         }
01044         else
01045         {
01046                 X = Matrix<ComplexFloat>(A.Columns(), B.Columns());
01047                 X.ReadFromMatrix(B);
01048         }
01049         int m = F.Rows();
01050         int nrhs = X.Columns();
01051         int n = F.Columns();
01052         int ldb = X.Rows();
01053 
01054         int min = (m>n)?n:m;
01055         int max = (m>n)?m:n;
01056         int bs = 2*min;
01057         int temp = 64*(n+1);
01058         bs = (bs>temp)?bs:temp;
01059         temp = min+64*min;
01060         bs = (bs>temp)?bs:temp;
01061         temp = min+64*nrhs;
01062         bs = (bs>temp)?bs:temp;
01063         bs += min;
01064         Vector<ComplexFloat> work(bs);
01065         Vector<float> rwork(2*n);
01066         Vector<int> jpvt(n);
01067         float rcond = max*numeric_limits<float>::epsilon();
01068 
01069         CGELSY(&m, &n, &nrhs, reinterpret_cast<MKL_Complex8*>(F.Data()), &m, reinterpret_cast<MKL_Complex8*>(X.Data()), &ldb, jpvt.Data(), &rcond, &rank, reinterpret_cast<MKL_Complex8*>(work.Data()), &bs, rwork.Data(), &info);
01070         if(rank < min && info == 0)
01071         {
01072                         char s[500];
01073                         sprintf(s, "Matrix is rank deficient (rank is %i)!", rank);
01074                         Utility::Warning(s);
01075         }
01076         if(info != 0)
01077         {
01078                 if(info < 0)
01079                 {
01080                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
01081                         char s[500];
01082                         sprintf(s, "%i'th argument of CGELSY is an illegal value!", -info);
01083                         Utility::RunTimeError(s);
01084                 }
01085                 else if(info > 0)  
01086                 {
01087                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
01088                         Utility::RunTimeError("Something went wrong when running CGELSY!");
01089                 }
01090         }
01091         if(B.Rows() > A.Columns())
01092         {
01093                 X = X(0, A.Columns()-1, 0, nrhs-1);
01094         }
01095         return X;
01096 }
01097 
01098 Matrix<ComplexDouble> Lapack::Gelsy(Matrix<ComplexDouble>& A, Matrix<ComplexDouble>& B)
01099 {
01100         if(A.Rows() != B.Rows())
01101         {
01102                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
01103                 Utility::RunTimeError("Matrix dimensions does not match!");
01104         }
01105 
01106         int info;
01107         int rank;
01108         Matrix<ComplexDouble> F = A.Clone();
01109         Matrix<ComplexDouble> X;
01110         if(B.Rows() >= A.Columns())
01111         {
01112                 X = B.Clone();
01113         }
01114         else
01115         {
01116                 X = Matrix<ComplexDouble>(A.Columns(), B.Columns());
01117                 X.ReadFromMatrix(B);
01118         }
01119         int m = F.Rows();
01120         int nrhs = X.Columns();
01121         int n = F.Columns();
01122         int ldb = X.Rows();
01123 
01124         int min = (m>n)?n:m;
01125         int max = (m>n)?m:n;
01126         int bs = 2*min;
01127         int temp = 64*(n+1);
01128         bs = (bs>temp)?bs:temp;
01129         temp = min+64*min;
01130         bs = (bs>temp)?bs:temp;
01131         temp = min+64*nrhs;
01132         bs = (bs>temp)?bs:temp;
01133         bs += min;
01134         Vector<ComplexDouble> work(bs);
01135         Vector<double> rwork(2*n);
01136         Vector<int> jpvt(n);
01137         double rcond = max*numeric_limits<double>::epsilon();
01138 
01139         ZGELSY(&m, &n, &nrhs, reinterpret_cast<MKL_Complex16*>(F.Data()), &m, reinterpret_cast<MKL_Complex16*>(X.Data()), &ldb, jpvt.Data(), &rcond, &rank, reinterpret_cast<MKL_Complex16*>(work.Data()), &bs, rwork.Data(), &info);
01140         if(rank < min && info == 0)
01141         {
01142                         char s[500];
01143                         sprintf(s, "Matrix is rank deficient (rank is %i)!", rank);
01144                         Utility::Warning(s);
01145         }
01146         if(info != 0)
01147         {
01148                 if(info < 0)
01149                 {
01150                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
01151                         char s[500];
01152                         sprintf(s, "%i'th argument of ZGELSY is an illegal value!", -info);
01153                         Utility::RunTimeError(s);
01154                 }
01155                 else if(info > 0)  
01156                 {
01157                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
01158                         Utility::RunTimeError("Something went wrong when running ZGELSY!");
01159                 }
01160         }
01161         if(B.Rows() > A.Columns())
01162         {
01163                 X = X(0, A.Columns()-1, 0, nrhs-1);
01164         }
01165         return X;
01166 }
01167 
01168 
01169 // Gelss
01170 Matrix<float> Lapack::Gelss(Matrix<float>& A, Matrix<float>& B)
01171 {
01172         if(A.Rows() != B.Rows())
01173         {
01174                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
01175                 Utility::RunTimeError("Matrix dimensions does not match!");
01176         }
01177 
01178         int info;
01179         int rank;
01180         Matrix<float> F = A.Clone();
01181         Matrix<float> X;
01182         if(B.Rows() >= A.Columns())
01183         {
01184                 X = B.Clone();
01185         }
01186         else
01187         {
01188                 X = Matrix<float>(A.Columns(), B.Columns());
01189                 X.ReadFromMatrix(B);
01190         }
01191         int m = F.Rows();
01192         int nrhs = X.Columns();
01193         int n = F.Columns();
01194         int ldb = X.Rows();
01195 
01196         int min = (m>n)?n:m;
01197         int max = (m>n)?m:n;
01198         int bs = max;
01199         int temp = 2*min;
01200         bs = (bs>temp)?bs:temp;
01201         bs = (bs>nrhs)?bs:nrhs;
01202         bs += 3*min;
01203         Vector<float> work(bs);
01204         Vector<float> s(min);
01205         float rcond = max*numeric_limits<float>::epsilon();
01206 
01207         SGELSS(&m, &n, &nrhs, F.Data(), &m, X.Data(), &ldb, s.Data(), &rcond, &rank, work.Data(), &bs, &info);
01208         if(rank < min && info == 0)
01209         {
01210                         char s[500];
01211                         sprintf(s, "Matrix is rank deficient (rank is %i)!", rank);
01212                         Utility::Warning(s);
01213         }
01214         if(info != 0)
01215         {
01216                 if(info < 0)
01217                 {
01218                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
01219                         char s[500];
01220                         sprintf(s, "%i'th argument of SGELSS is an illegal value!", -info);
01221                         Utility::RunTimeError(s);
01222                 }
01223                 else if(info > 0)  
01224                 {
01225                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
01226                         char s[500];
01227                         sprintf(s, "The algorithm for computing the SVD failed to converge; %i off-diagonal elements of an intermediate bidiagonal form did not converge to zero!", info);
01228                         Utility::RunTimeError(s);
01229                 }
01230         }
01231         if(B.Rows() > A.Columns())
01232         {
01233                 X = X(0, A.Columns()-1, 0, nrhs-1);
01234         }
01235         return X;
01236 }
01237 
01238 Matrix<double> Lapack::Gelss(Matrix<double>& A, Matrix<double>& B)
01239 {
01240         if(A.Rows() != B.Rows())
01241         {
01242                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
01243                 Utility::RunTimeError("Matrix dimensions does not match!");
01244         }
01245 
01246         int info;
01247         int rank;
01248         Matrix<double> F = A.Clone();
01249         Matrix<double> X;
01250         if(B.Rows() >= A.Columns())
01251         {
01252                 X = B.Clone();
01253         }
01254         else
01255         {
01256                 X = Matrix<double>(A.Columns(), B.Columns());
01257                 X.ReadFromMatrix(B);
01258         }
01259         int m = F.Rows();
01260         int nrhs = X.Columns();
01261         int n = F.Columns();
01262         int ldb = X.Rows();
01263 
01264         int min = (m>n)?n:m;
01265         int max = (m>n)?m:n;
01266         int bs = max;
01267         int temp = 2*min;
01268         bs = (bs>temp)?bs:temp;
01269         bs = (bs>nrhs)?bs:nrhs;
01270         bs += 3*min;
01271         Vector<double> work(bs);
01272         Vector<double> s(min);
01273         double rcond = max*numeric_limits<double>::epsilon();
01274 
01275         DGELSS(&m, &n, &nrhs, F.Data(), &m, X.Data(), &ldb, s.Data(), &rcond, &rank, work.Data(), &bs, &info);
01276         if(rank < min && info == 0)
01277         {
01278                         char s[500];
01279                         sprintf(s, "Matrix is rank deficient (rank is %i)!", rank);
01280                         Utility::Warning(s);
01281         }
01282         if(info != 0)
01283         {
01284                 if(info < 0)
01285                 {
01286                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
01287                         char s[500];
01288                         sprintf(s, "%i'th argument of DGELSS is an illegal value!", -info);
01289                         Utility::RunTimeError(s);
01290                 }
01291                 else if(info > 0)  
01292                 {
01293                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
01294                         char s[500];
01295                         sprintf(s, "The algorithm for computing the SVD failed to converge; %i off-diagonal elements of an intermediate bidiagonal form did not converge to zero!", info);
01296                         Utility::RunTimeError(s);
01297                 }
01298         }
01299         if(B.Rows() > A.Columns())
01300         {
01301                 X = X(0, A.Columns()-1, 0, nrhs-1);
01302         }
01303         return X;
01304 }
01305 
01306 Matrix<ComplexFloat> Lapack::Gelss(Matrix<ComplexFloat>& A, Matrix<ComplexFloat>& B)
01307 {
01308         if(A.Rows() != B.Rows())
01309         {
01310                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
01311                 Utility::RunTimeError("Matrix dimensions does not match!");
01312         }
01313 
01314         int info;
01315         int rank;
01316         Matrix<ComplexFloat> F = A.Clone();
01317         Matrix<ComplexFloat> X;
01318         if(B.Rows() >= A.Columns())
01319         {
01320                 X = B.Clone();
01321         }
01322         else
01323         {
01324                 X = Matrix<ComplexFloat>(A.Columns(), B.Columns());
01325                 X.ReadFromMatrix(B);
01326         }
01327         int m = F.Rows();
01328         int nrhs = X.Columns();
01329         int n = F.Columns();
01330         int ldb = X.Rows();
01331 
01332         int min = (m>n)?n:m;
01333         int max = (m>n)?m:n;
01334         int bs = (max>nrhs)?max:nrhs;
01335         bs += 2*min;
01336         Vector<ComplexFloat> work(bs);
01337         Vector<float> s(min);
01338         Vector<float> rwork(5*min);
01339         float rcond = max*numeric_limits<float>::epsilon();
01340 
01341         CGELSS(&m, &n, &nrhs, reinterpret_cast<MKL_Complex8*>(F.Data()), &m, reinterpret_cast<MKL_Complex8*>(X.Data()), &ldb, s.Data(), &rcond, &rank, reinterpret_cast<MKL_Complex8*>(work.Data()), &bs, rwork.Data(), &info);
01342         if(rank < min && info == 0)
01343         {
01344                         char s[500];
01345                         sprintf(s, "Matrix is rank deficient (rank is %i)!", rank);
01346                         Utility::Warning(s);
01347         }
01348         if(info != 0)
01349         {
01350                 if(info < 0)
01351                 {
01352                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
01353                         char s[500];
01354                         sprintf(s, "%i'th argument of CGELSS is an illegal value!", -info);
01355                         Utility::RunTimeError(s);
01356                 }
01357                 else if(info > 0)  
01358                 {
01359                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
01360                         char s[500];
01361                         sprintf(s, "The algorithm for computing the SVD failed to converge; %i off-diagonal elements of an intermediate bidiagonal form did not converge to zero!", info);
01362                         Utility::RunTimeError(s);
01363                 }
01364         }
01365         if(B.Rows() > A.Columns())
01366         {
01367                 X = X(0, A.Columns()-1, 0, nrhs-1);
01368         }
01369         return X;
01370 }
01371 
01372 Matrix<ComplexDouble> Lapack::Gelss(Matrix<ComplexDouble>& A, Matrix<ComplexDouble>& B)
01373 {
01374         if(A.Rows() != B.Rows())
01375         {
01376                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
01377                 Utility::RunTimeError("Matrix dimensions does not match!");
01378         }
01379 
01380         int info;
01381         int rank;
01382         Matrix<ComplexDouble> F = A.Clone();
01383         Matrix<ComplexDouble> X;
01384         if(B.Rows() >= A.Columns())
01385         {
01386                 X = B.Clone();
01387         }
01388         else
01389         {
01390                 X = Matrix<ComplexDouble>(A.Columns(), B.Columns());
01391                 X.ReadFromMatrix(B);
01392         }
01393         int m = F.Rows();
01394         int nrhs = X.Columns();
01395         int n = F.Columns();
01396         int ldb = X.Rows();
01397 
01398         int min = (m>n)?n:m;
01399         int max = (m>n)?m:n;
01400         int bs = (max>nrhs)?max:nrhs;
01401         bs += 2*min;
01402         Vector<ComplexDouble> work(bs);
01403         Vector<double> s(min);
01404         Vector<double> rwork(5*min);
01405         double rcond = max*numeric_limits<double>::epsilon();
01406 
01407         ZGELSS(&m, &n, &nrhs, reinterpret_cast<MKL_Complex16*>(F.Data()), &m, reinterpret_cast<MKL_Complex16*>(X.Data()), &ldb, s.Data(), &rcond, &rank, reinterpret_cast<MKL_Complex16*>(work.Data()), &bs, rwork.Data(), &info);
01408         if(rank < min && info == 0)
01409         {
01410                         char s[500];
01411                         sprintf(s, "Matrix is rank deficient (rank is %i)!", rank);
01412                         Utility::Warning(s);
01413         }
01414         if(info != 0)
01415         {
01416                 if(info < 0)
01417                 {
01418                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
01419                         char s[500];
01420                         sprintf(s, "%i'th argument of ZGELSS is an illegal value!", -info);
01421                         Utility::RunTimeError(s);
01422                 }
01423                 else if(info > 0)  
01424                 {
01425                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
01426                         char s[500];
01427                         sprintf(s, "The algorithm for computing the SVD failed to converge; %i off-diagonal elements of an intermediate bidiagonal form did not converge to zero!", info);
01428                         Utility::RunTimeError(s);
01429                 }
01430         }
01431         if(B.Rows() > A.Columns())
01432         {
01433                 X = X(0, A.Columns()-1, 0, nrhs-1);
01434         }
01435         return X;
01436 }
01437 
01438 // Gelsd
01439 Matrix<float> Lapack::Gelsd(Matrix<float>& A, Matrix<float>& B)
01440 {
01441         if(A.Rows() != B.Rows())
01442         {
01443                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
01444                 Utility::RunTimeError("Matrix dimensions does not match!");
01445         }
01446 
01447         int info;
01448         int rank;
01449         Matrix<float> F = A.Clone();
01450         Matrix<float> X;
01451         if(B.Rows() >= A.Columns())
01452         {
01453                 X = B.Clone();
01454         }
01455         else
01456         {
01457                 X = Matrix<float>(A.Columns(), B.Columns());
01458                 X.ReadFromMatrix(B);
01459         }
01460         int m = F.Rows();
01461         int nrhs = X.Columns();
01462         int n = F.Columns();
01463         int ldb = X.Rows();
01464 
01465         int min = (m>n)?n:m;
01466         int max = (m>n)?m:n;
01467         int smlsiz = 25;
01468         int nlvl = (int)(log((double)min/(double)(smlsiz+1))/log(2.0))+1;
01469         nlvl = (nlvl>0)?nlvl:0;
01470         int bs = 0;
01471         if(m>=n)
01472         {
01473                 bs = 12*n + 2*n*smlsiz + 8*n*nlvl + n*nrhs + (smlsiz+1)*(smlsiz+1);
01474         }
01475         else
01476         {
01477                 bs = 12*m + 2*m*smlsiz + 8*m*nlvl + m*nrhs + (smlsiz+1)*(smlsiz+1);
01478         }
01479         Vector<float> work(bs);
01480         int iwork_size = 3*min*nlvl + 11*min;
01481         Vector<int> iwork(iwork_size);
01482         Vector<float> s(min);
01483         float rcond = max*numeric_limits<float>::epsilon();
01484 
01485         SGELSD(&m, &n, &nrhs, F.Data(), &m, X.Data(), &ldb, s.Data(), &rcond, &rank, work.Data(), &bs, iwork.Data(), &info);
01486 
01487         if(rank < min && info == 0)
01488         {
01489                         char s[500];
01490                         sprintf(s, "Matrix is rank deficient (rank is %i)!", rank);
01491                         Utility::Warning(s);
01492         }
01493         if(info != 0)
01494         {
01495                 if(info < 0)
01496                 {
01497                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
01498                         char s[500];
01499                         sprintf(s, "%i'th argument of SGELSD is an illegal value!", -info);
01500                         Utility::RunTimeError(s);
01501                 }
01502                 else if(info > 0)  
01503                 {
01504                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
01505                         char s[500];
01506                         sprintf(s, "The algorithm for computing the SVD failed to converge; %i off-diagonal elements of an intermediate bidiagonal form did not converge to zero!", info);
01507                         Utility::RunTimeError(s);
01508                 }
01509         }
01510         if(B.Rows() > A.Columns())
01511         {
01512                 X = X(0, A.Columns()-1, 0, nrhs-1);
01513         }
01514         return X;
01515 }
01516 
01517 Matrix<double> Lapack::Gelsd(Matrix<double>& A, Matrix<double>& B)
01518 {
01519         if(A.Rows() != B.Rows())
01520         {
01521                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
01522                 Utility::RunTimeError("Matrix dimensions does not match!");
01523         }
01524 
01525         int info;
01526         int rank;
01527         Matrix<double> F = A.Clone();
01528         Matrix<double> X;
01529         if(B.Rows() >= A.Columns())
01530         {
01531                 X = B.Clone();
01532         }
01533         else
01534         {
01535                 X = Matrix<double>(A.Columns(), B.Columns());
01536                 X.ReadFromMatrix(B);
01537         }
01538         int m = F.Rows();
01539         int nrhs = X.Columns();
01540         int n = F.Columns();
01541         int ldb = X.Rows();
01542 
01543         int min = (m>n)?n:m;
01544         int max = (m>n)?m:n;
01545         int smlsiz = 25;
01546         int nlvl = (int)(log((double)min/(double)(smlsiz+1))/log(2.0))+1;
01547         nlvl = (nlvl>0)?nlvl:0;
01548         int bs = 0;
01549         if(m>=n)
01550         {
01551                 bs = 12*n + 2*n*smlsiz + 8*n*nlvl + n*nrhs + (smlsiz+1)*(smlsiz+1);
01552         }
01553         else
01554         {
01555                 bs = 12*m + 2*m*smlsiz + 8*m*nlvl + m*nrhs + (smlsiz+1)*(smlsiz+1);
01556         }
01557         Vector<double> work(bs);
01558         int iwork_size = 3*min*nlvl + 11*min;
01559         Vector<int> iwork(iwork_size);
01560         Vector<double> s(min);
01561         double rcond = max*numeric_limits<double>::epsilon();
01562 
01563         DGELSD(&m, &n, &nrhs, F.Data(), &m, X.Data(), &ldb, s.Data(), &rcond, &rank, work.Data(), &bs, iwork.Data(), &info);
01564 
01565         if(rank < min && info == 0)
01566         {
01567                         char s[500];
01568                         sprintf(s, "Matrix is rank deficient (rank is %i)!", rank);
01569                         Utility::Warning(s);
01570         }
01571         if(info != 0)
01572         {
01573                 if(info < 0)
01574                 {
01575                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
01576                         char s[500];
01577                         sprintf(s, "%i'th argument of DGELSD is an illegal value!", -info);
01578                         Utility::RunTimeError(s);
01579                 }
01580                 else if(info > 0)  
01581                 {
01582                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
01583                         char s[500];
01584                         sprintf(s, "The algorithm for computing the SVD failed to converge; %i off-diagonal elements of an intermediate bidiagonal form did not converge to zero!", info);
01585                         Utility::RunTimeError(s);
01586                 }
01587         }
01588         if(B.Rows() > A.Columns())
01589         {
01590                 X = X(0, A.Columns()-1, 0, nrhs-1);
01591         }
01592         return X;
01593 }
01594 
01595 Matrix<ComplexFloat> Lapack::Gelsd(Matrix<ComplexFloat>& A, Matrix<ComplexFloat>& B)
01596 {
01597         if(A.Rows() != B.Rows())
01598         {
01599                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
01600                 Utility::RunTimeError("Matrix dimensions does not match!");
01601         }
01602 
01603         int info;
01604         int rank;
01605         Matrix<ComplexFloat> F = A.Clone();
01606         Matrix<ComplexFloat> X;
01607         if(B.Rows() >= A.Columns())
01608         {
01609                 X = B.Clone();
01610         }
01611         else
01612         {
01613                 X = Matrix<ComplexFloat>(A.Columns(), B.Columns());
01614                 X.ReadFromMatrix(B);
01615         }
01616         int m = F.Rows();
01617         int nrhs = X.Columns();
01618         int n = F.Columns();
01619         int ldb = X.Rows();
01620 
01621         int min = (m>n)?n:m;
01622         int max = (m>n)?m:n;
01623         int smlsiz = 25;
01624         int nlvl = (int)(log((double)min/(double)(smlsiz+1))/log(2.0))+1;
01625         nlvl = (nlvl>0)?nlvl:0;
01626         int bs = -1;
01627 
01628         Vector<ComplexFloat> work(1);
01629         int lrwork = 0;
01630         if(m>=n)
01631         {
01632                 lrwork = 10*n + 2*n*smlsiz + 8*n*nlvl + 3*smlsiz*nrhs + (smlsiz+1)*(smlsiz+1);
01633         }
01634         else
01635         {
01636                 lrwork = 10*m + 2*m*smlsiz + 8*m*nlvl + 3*smlsiz*nrhs + (smlsiz+1)*(smlsiz+1);
01637         }
01638         Vector<float> rwork(lrwork);
01639         int iwork_size = 3*min*nlvl + 11*min;
01640         Vector<int> iwork(iwork_size);
01641         Vector<float> s(min);
01642         float rcond = max*numeric_limits<float>::epsilon();
01643 
01644         CGELSD(&m, &n, &nrhs, reinterpret_cast<MKL_Complex8*>(F.Data()), &m, reinterpret_cast<MKL_Complex8*>(X.Data()), &ldb, s.Data(), &rcond, &rank, reinterpret_cast<MKL_Complex8*>(work.Data()), &bs, rwork.Data(), iwork.Data(), &info);
01645         bs = (int)real(work[0]);
01646         work = Vector<ComplexFloat>(bs);
01647         CGELSD(&m, &n, &nrhs, reinterpret_cast<MKL_Complex8*>(F.Data()), &m, reinterpret_cast<MKL_Complex8*>(X.Data()), &ldb, s.Data(), &rcond, &rank, reinterpret_cast<MKL_Complex8*>(work.Data()), &bs, rwork.Data(), iwork.Data(), &info);
01648 
01649         if(rank < min && info == 0)
01650         {
01651                         char s[500];
01652                         sprintf(s, "Matrix is rank deficient (rank is %i)!", rank);
01653                         Utility::Warning(s);
01654         }
01655         if(info != 0)
01656         {
01657                 if(info < 0)
01658                 {
01659                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
01660                         char s[500];
01661                         sprintf(s, "%i'th argument of CGELSD is an illegal value!", -info);
01662                         Utility::RunTimeError(s);
01663                 }
01664                 else if(info > 0)  
01665                 {
01666                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
01667                         char s[500];
01668                         sprintf(s, "The algorithm for computing the SVD failed to converge; %i off-diagonal elements of an intermediate bidiagonal form did not converge to zero!", info);
01669                         Utility::RunTimeError(s);
01670                 }
01671         }
01672         if(B.Rows() > A.Columns())
01673         {
01674                 X = X(0, A.Columns()-1, 0, nrhs-1);
01675         }
01676         return X;
01677 }
01678 
01679 Matrix<ComplexDouble> Lapack::Gelsd(Matrix<ComplexDouble>& A, Matrix<ComplexDouble>& B)
01680 {
01681         if(A.Rows() != B.Rows())
01682         {
01683                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
01684                 Utility::RunTimeError("Matrix dimensions does not match!");
01685         }
01686 
01687         int info;
01688         int rank;
01689         Matrix<ComplexDouble> F = A.Clone();
01690         Matrix<ComplexDouble> X;
01691         if(B.Rows() >= A.Columns())
01692         {
01693                 X = B.Clone();
01694         }
01695         else
01696         {
01697                 X = Matrix<ComplexDouble>(A.Columns(), B.Columns());
01698                 X.ReadFromMatrix(B);
01699         }
01700         int m = F.Rows();
01701         int nrhs = X.Columns();
01702         int n = F.Columns();
01703         int ldb = X.Rows();
01704 
01705         int min = (m>n)?n:m;
01706         int max = (m>n)?m:n;
01707         int smlsiz = 25;
01708         int nlvl = (int)(log((double)min/(double)(smlsiz+1))/log(2.0))+1;
01709         nlvl = (nlvl>0)?nlvl:0;
01710         int bs = -1;
01711 
01712         Vector<ComplexDouble> work(1);
01713         int lrwork = 0;
01714         if(m>=n)
01715         {
01716                 lrwork = 10*n + 2*n*smlsiz + 8*n*nlvl + 3*smlsiz*nrhs + (smlsiz+1)*(smlsiz+1);
01717         }
01718         else
01719         {
01720                 lrwork = 10*m + 2*m*smlsiz + 8*m*nlvl + 3*smlsiz*nrhs + (smlsiz+1)*(smlsiz+1);
01721         }
01722         Vector<double> rwork(lrwork);
01723         int iwork_size = 3*min*nlvl + 11*min;
01724         Vector<int> iwork(iwork_size);
01725         Vector<double> s(min);
01726         double rcond = max*numeric_limits<double>::epsilon();
01727         
01728         ZGELSD(&m, &n, &nrhs, reinterpret_cast<MKL_Complex16*>(F.Data()), &m, reinterpret_cast<MKL_Complex16*>(X.Data()), &ldb, s.Data(), &rcond, &rank, reinterpret_cast<MKL_Complex16*>(work.Data()), &bs, rwork.Data(), iwork.Data(), &info);
01729         bs = (int)real(work[0]);
01730         work = Vector<ComplexDouble>(bs);
01731         ZGELSD(&m, &n, &nrhs, reinterpret_cast<MKL_Complex16*>(F.Data()), &m, reinterpret_cast<MKL_Complex16*>(X.Data()), &ldb, s.Data(), &rcond, &rank, reinterpret_cast<MKL_Complex16*>(work.Data()), &bs, rwork.Data(), iwork.Data(), &info);
01732 
01733         if(rank < min && info == 0)
01734         {
01735                         char s[500];
01736                         sprintf(s, "Matrix is rank deficient (rank is %i)!", rank);
01737                         Utility::Warning(s);
01738         }
01739         if(info != 0)
01740         {
01741                 if(info < 0)
01742                 {
01743                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
01744                         char s[500];
01745                         sprintf(s, "%i'th argument of ZGELSD is an illegal value!", -info);
01746                         Utility::RunTimeError(s);
01747                 }
01748                 else if(info > 0)  
01749                 {
01750                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
01751                         char s[500];
01752                         sprintf(s, "The algorithm for computing the SVD failed to converge; %i off-diagonal elements of an intermediate bidiagonal form did not converge to zero!", info);
01753                         Utility::RunTimeError(s);
01754                 }
01755         }
01756         if(B.Rows() > A.Columns())
01757         {
01758                 X = X(0, A.Columns()-1, 0, nrhs-1);
01759         }
01760         return X;
01761 }
01762 
01763 
01764 
01765 
01766 
01767 
01768 //=========================================================
01769 // Generalized Linear Least Squares problems
01770 //=========================================================
01771 
01772 
01773 // Gglse
01774 Vector<float> Lapack::Gglse(Matrix<float>& A, Matrix<float>& B, Vector<float>& c, Vector<float>& d)
01775 {
01776         if(A.Columns() != B.Columns())
01777         {
01778                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
01779                 Utility::RunTimeError("Matrix dimensions does not match!");
01780         }
01781         if(A.Rows() != c.Length() || B.Rows() != d.Length())
01782         {
01783                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
01784                 Utility::RunTimeError("Matrix/Vector dimensions does not match!");
01785         }
01786 
01787         if(B.Rows() > A.Columns() || (B.Rows() + A.Rows()) < A.Columns())
01788         {
01789                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
01790                 Utility::RunTimeError("Matrix/Vector dimensions have incompatibilities!");
01791         }
01792 
01793         int info;
01794         Matrix<float> A2 = A.Clone();
01795         Matrix<float> B2 = B.Clone();
01796         Vector<float> c2 = c.Clone();
01797         Vector<float> d2 = d.Clone();
01798 
01799         int m = A2.Rows();
01800         int n = A2.Columns();
01801         int p = B2.Rows();
01802         Vector<float> x(n);
01803 
01804         int min = (m>n)?n:m;
01805         int max = (m>n)?m:n;
01806         int nb = 64;
01807         int bs = p + min + max*nb;
01808         Vector<float> work(bs);
01809 
01810         SGGLSE(&m, &n, &p, A2.Data(), &m, B2.Data(), &p, c2.Data(), d2.Data(), x.Data(), work.Data(), &bs, &info);
01811 
01812         if(info != 0)
01813         {
01814                 if(info < 0)
01815                 {
01816                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
01817                         char s[500];
01818                         sprintf(s, "%i'th argument of SGGLSE is an illegal value!", -info);
01819                         Utility::RunTimeError(s);
01820                 }
01821                 else if(info > 0)  
01822                 {
01823                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
01824                         char s[500];
01825                         sprintf(s, "Problem experienced when running SGGLSE!");
01826                         Utility::RunTimeError(s);
01827                 }
01828         }
01829         return x;
01830 }
01831 
01832 
01833 
01834 Vector<double> Lapack::Gglse(Matrix<double>& A, Matrix<double>& B, Vector<double>& c, Vector<double>& d)
01835 {
01836         if(A.Columns() != B.Columns())
01837         {
01838                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
01839                 Utility::RunTimeError("Matrix dimensions does not match!");
01840         }
01841         if(A.Rows() != c.Length() || B.Rows() != d.Length())
01842         {
01843                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
01844                 Utility::RunTimeError("Matrix/Vector dimensions does not match!");
01845         }
01846 
01847         if(B.Rows() > A.Columns() || (B.Rows() + A.Rows()) < A.Columns())
01848         {
01849                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
01850                 Utility::RunTimeError("Matrix/Vector dimensions have incompatibilities!");
01851         }
01852 
01853         int info;
01854         Matrix<double> A2 = A.Clone();
01855         Matrix<double> B2 = B.Clone();
01856         Vector<double> c2 = c.Clone();
01857         Vector<double> d2 = d.Clone();
01858 
01859         int m = A2.Rows();
01860         int n = A2.Columns();
01861         int p = B2.Rows();
01862         Vector<double> x(n);
01863 
01864         int min = (m>n)?n:m;
01865         int max = (m>n)?m:n;
01866         int nb = 64;
01867         int bs = p + min + max*nb;
01868         Vector<double> work(bs);
01869 
01870         DGGLSE(&m, &n, &p, A2.Data(), &m, B2.Data(), &p, c2.Data(), d2.Data(), x.Data(), work.Data(), &bs, &info);
01871 
01872         if(info != 0)
01873         {
01874                 if(info < 0)
01875                 {
01876                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
01877                         char s[500];
01878                         sprintf(s, "%i'th argument of DGGLSE is an illegal value!", -info);
01879                         Utility::RunTimeError(s);
01880                 }
01881                 else if(info > 0)  
01882                 {
01883                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
01884                         char s[500];
01885                         sprintf(s, "Problem experienced when running DGGLSE!");
01886                         Utility::RunTimeError(s);
01887                 }
01888         }
01889         return x;
01890 }
01891 
01892 
01893 Vector<ComplexFloat> Lapack::Gglse(Matrix<ComplexFloat>& A, Matrix<ComplexFloat>& B, Vector<ComplexFloat>& c, Vector<ComplexFloat>& d)
01894 {
01895         if(A.Columns() != B.Columns())
01896         {
01897                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
01898                 Utility::RunTimeError("Matrix dimensions does not match!");
01899         }
01900         if(A.Rows() != c.Length() || B.Rows() != d.Length())
01901         {
01902                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
01903                 Utility::RunTimeError("Matrix/Vector dimensions does not match!");
01904         }
01905 
01906         if(B.Rows() > A.Columns() || (B.Rows() + A.Rows()) < A.Columns())
01907         {
01908                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
01909                 Utility::RunTimeError("Matrix/Vector dimensions have incompatibilities!");
01910         }
01911 
01912         int info;
01913         Matrix<ComplexFloat> A2 = A.Clone();
01914         Matrix<ComplexFloat> B2 = B.Clone();
01915         Vector<ComplexFloat> c2 = c.Clone();
01916         Vector<ComplexFloat> d2 = d.Clone();
01917 
01918         int m = A2.Rows();
01919         int n = A2.Columns();
01920         int p = B2.Rows();
01921         Vector<ComplexFloat> x(n);
01922 
01923         int min = (m>n)?n:m;
01924         int max = (m>n)?m:n;
01925         int nb = 64;
01926         int bs = p + min + max*nb;
01927         Vector<ComplexFloat> work(bs);
01928 
01929         CGGLSE(&m, &n, &p, reinterpret_cast<MKL_Complex8*>(A2.Data()), &m, reinterpret_cast<MKL_Complex8*>(B2.Data()), &p, reinterpret_cast<MKL_Complex8*>(c2.Data()), reinterpret_cast<MKL_Complex8*>(d2.Data()), reinterpret_cast<MKL_Complex8*>(x.Data()), reinterpret_cast<MKL_Complex8*>(work.Data()), &bs, &info);
01930 
01931         if(info != 0)
01932         {
01933                 if(info < 0)
01934                 {
01935                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
01936                         char s[500];
01937                         sprintf(s, "%i'th argument of CGGLSE is an illegal value!", -info);
01938                         Utility::RunTimeError(s);
01939                 }
01940                 else if(info > 0)  
01941                 {
01942                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
01943                         char s[500];
01944                         sprintf(s, "Problem experienced when running CGGLSE!");
01945                         Utility::RunTimeError(s);
01946                 }
01947         }
01948         return x;
01949 }
01950 
01951 
01952 Vector<ComplexDouble> Lapack::Gglse(Matrix<ComplexDouble>& A, Matrix<ComplexDouble>& B, Vector<ComplexDouble>& c, Vector<ComplexDouble>& d)
01953 {
01954         if(A.Columns() != B.Columns())
01955         {
01956                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
01957                 Utility::RunTimeError("Matrix dimensions does not match!");
01958         }
01959         if(A.Rows() != c.Length() || B.Rows() != d.Length())
01960         {
01961                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
01962                 Utility::RunTimeError("Matrix/Vector dimensions does not match!");
01963         }
01964 
01965         if(B.Rows() > A.Columns() || (B.Rows() + A.Rows()) < A.Columns())
01966         {
01967                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
01968                 Utility::RunTimeError("Matrix/Vector dimensions have incompatibilities!");
01969         }
01970 
01971         int info;
01972         Matrix<ComplexDouble> A2 = A.Clone();
01973         Matrix<ComplexDouble> B2 = B.Clone();
01974         Vector<ComplexDouble> c2 = c.Clone();
01975         Vector<ComplexDouble> d2 = d.Clone();
01976 
01977         int m = A2.Rows();
01978         int n = A2.Columns();
01979         int p = B2.Rows();
01980         Vector<ComplexDouble> x(n);
01981 
01982         int min = (m>n)?n:m;
01983         int max = (m>n)?m:n;
01984         int nb = 64;
01985         int bs = p + min + max*nb;
01986         Vector<ComplexDouble> work(bs);
01987 
01988         ZGGLSE(&m, &n, &p, reinterpret_cast<MKL_Complex16*>(A2.Data()), &m, reinterpret_cast<MKL_Complex16*>(B2.Data()), &p, reinterpret_cast<MKL_Complex16*>(c2.Data()), reinterpret_cast<MKL_Complex16*>(d2.Data()), reinterpret_cast<MKL_Complex16*>(x.Data()), reinterpret_cast<MKL_Complex16*>(work.Data()), &bs, &info);
01989 
01990         if(info != 0)
01991         {
01992                 if(info < 0)
01993                 {
01994                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
01995                         char s[500];
01996                         sprintf(s, "%i'th argument of ZGGLSE is an illegal value!", -info);
01997                         Utility::RunTimeError(s);
01998                 }
01999                 else if(info > 0)  
02000                 {
02001                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
02002                         char s[500];
02003                         sprintf(s, "Problem experienced when running ZGGLSE!");
02004                         Utility::RunTimeError(s);
02005                 }
02006         }
02007         return x;
02008 }
02009 
02010 
02011 
02012 
02013 
02014 // Ggglm
02015 void Lapack::Ggglm(Matrix<float>& A, Matrix<float>& B, Vector<float>& d, Vector<float>& x, Vector<float>& y)
02016 {
02017         if(A.Rows() != B.Rows() || B.Columns() < (A.Rows() - A.Columns()))
02018         {
02019                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
02020                 Utility::RunTimeError("Matrix dimensions does not match!");
02021         }
02022 
02023         if(d.Length() != A.Rows())
02024         {
02025                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
02026                 Utility::RunTimeError("Matrix/Vector dimensions have incompatibilities!");
02027         }
02028 
02029         int info;
02030         Matrix<float> A2 = A.Clone();
02031         Matrix<float> B2 = B.Clone();
02032         Vector<float> d2 = d.Clone();
02033 
02034         int m = A2.Columns();
02035         int n = A2.Rows();
02036         int p = B2.Columns();
02037         x = Vector<float>(m);
02038         y = Vector<float>(p);
02039 
02040         int min = (n>p)?p:n;
02041         int max = (n>p)?n:p;
02042         int nb = 64;
02043         int bs = m + min + max*nb;
02044         Vector<float> work(bs);
02045 
02046         SGGGLM(&n, &m, &p, A2.Data(), &n, B2.Data(), &n, d2.Data(), x.Data(), y.Data(), work.Data(), &bs, &info);
02047 
02048         if(info != 0)
02049         {
02050                 if(info < 0)
02051                 {
02052                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
02053                         char s[500];
02054                         sprintf(s, "%i'th argument of SGGGLM is an illegal value!", -info);
02055                         Utility::RunTimeError(s);
02056                 }
02057                 else if(info > 0)  
02058                 {
02059                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
02060                         char s[500];
02061                         sprintf(s, "Problem experienced when running SGGGLM!");
02062                         Utility::RunTimeError(s);
02063                 }
02064         }
02065 
02066 }
02067 
02068 
02069 void Lapack::Ggglm(Matrix<double>& A, Matrix<double>& B, Vector<double>& d, Vector<double>& x, Vector<double>& y)
02070 {
02071         if(A.Rows() != B.Rows() || B.Columns() < (A.Rows() - A.Columns()))
02072         {
02073                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
02074                 Utility::RunTimeError("Matrix dimensions does not match!");
02075         }
02076 
02077         if(d.Length() != A.Rows())
02078         {
02079                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
02080                 Utility::RunTimeError("Matrix/Vector dimensions have incompatibilities!");
02081         }
02082 
02083         int info;
02084         Matrix<double> A2 = A.Clone();
02085         Matrix<double> B2 = B.Clone();
02086         Vector<double> d2 = d.Clone();
02087 
02088         int m = A2.Columns();
02089         int n = A2.Rows();
02090         int p = B2.Columns();
02091         x = Vector<double>(m);
02092         y = Vector<double>(p);
02093 
02094         int min = (n>p)?p:n;
02095         int max = (n>p)?n:p;
02096         int nb = 64;
02097         int bs = m + min + max*nb;
02098         Vector<double> work(bs);
02099 
02100         DGGGLM(&n, &m, &p, A2.Data(), &n, B2.Data(), &n, d2.Data(), x.Data(), y.Data(), work.Data(), &bs, &info);
02101 
02102         if(info != 0)
02103         {
02104                 if(info < 0)
02105                 {
02106                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
02107                         char s[500];
02108                         sprintf(s, "%i'th argument of DGGGLM is an illegal value!", -info);
02109                         Utility::RunTimeError(s);
02110                 }
02111                 else if(info > 0)  
02112                 {
02113                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
02114                         char s[500];
02115                         sprintf(s, "Problem experienced when running DGGGLM!");
02116                         Utility::RunTimeError(s);
02117                 }
02118         }
02119 }
02120 
02121 
02122 void Lapack::Ggglm(Matrix<ComplexFloat>& A, Matrix<ComplexFloat>& B, Vector<ComplexFloat>& d, Vector<ComplexFloat>& x, Vector<ComplexFloat>& y)
02123 {
02124         if(A.Rows() != B.Rows() || B.Columns() < (A.Rows() - A.Columns()))
02125         {
02126                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
02127                 Utility::RunTimeError("Matrix dimensions does not match!");
02128         }
02129 
02130         if(d.Length() != A.Rows())
02131         {
02132                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
02133                 Utility::RunTimeError("Matrix/Vector dimensions have incompatibilities!");
02134         }
02135 
02136         int info;
02137         Matrix<ComplexFloat> A2 = A.Clone();
02138         Matrix<ComplexFloat> B2 = B.Clone();
02139         Vector<ComplexFloat> d2 = d.Clone();
02140 
02141         int m = A2.Columns();
02142         int n = A2.Rows();
02143         int p = B2.Columns();
02144         x = Vector<ComplexFloat>(m);
02145         y = Vector<ComplexFloat>(p);
02146 
02147         int min = (n>p)?p:n;
02148         int max = (n>p)?n:p;
02149         int nb = 64;
02150         int bs = m + min + max*nb;
02151         Vector<ComplexFloat> work(bs);
02152 
02153         CGGGLM(&n, &m, &p, reinterpret_cast<MKL_Complex8*>(A2.Data()), &n, reinterpret_cast<MKL_Complex8*>(B2.Data()), &n, reinterpret_cast<MKL_Complex8*>(d2.Data()), reinterpret_cast<MKL_Complex8*>(x.Data()), reinterpret_cast<MKL_Complex8*>(y.Data()), reinterpret_cast<MKL_Complex8*>(work.Data()), &bs, &info);
02154 
02155         if(info != 0)
02156         {
02157                 if(info < 0)
02158                 {
02159                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
02160                         char s[500];
02161                         sprintf(s, "%i'th argument of CGGGLM is an illegal value!", -info);
02162                         Utility::RunTimeError(s);
02163                 }
02164                 else if(info > 0)  
02165                 {
02166                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
02167                         char s[500];
02168                         sprintf(s, "Problem experienced when running CGGGLM!");
02169                         Utility::RunTimeError(s);
02170                 }
02171         }
02172 }
02173 
02174 
02175 void Lapack::Ggglm(Matrix<ComplexDouble>& A, Matrix<ComplexDouble>& B, Vector<ComplexDouble>& d, Vector<ComplexDouble>& x, Vector<ComplexDouble>& y)
02176 {
02177         if(A.Rows() != B.Rows() || B.Columns() < (A.Rows() - A.Columns()))
02178         {
02179                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
02180                 Utility::RunTimeError("Matrix dimensions does not match!");
02181         }
02182 
02183         if(d.Length() != A.Rows())
02184         {
02185                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
02186                 Utility::RunTimeError("Matrix/Vector dimensions have incompatibilities!");
02187         }
02188 
02189         int info;
02190         Matrix<ComplexDouble> A2 = A.Clone();
02191         Matrix<ComplexDouble> B2 = B.Clone();
02192         Vector<ComplexDouble> d2 = d.Clone();
02193 
02194         int m = A2.Columns();
02195         int n = A2.Rows();
02196         int p = B2.Columns();
02197         x = Vector<ComplexDouble>(m);
02198         y = Vector<ComplexDouble>(p);
02199 
02200         int min = (n>p)?p:n;
02201         int max = (n>p)?n:p;
02202         int nb = 64;
02203         int bs = m + min + max*nb;
02204         Vector<ComplexDouble> work(bs);
02205 
02206         ZGGGLM(&n, &m, &p, reinterpret_cast<MKL_Complex16*>(A2.Data()), &n, reinterpret_cast<MKL_Complex16*>(B2.Data()), &n, reinterpret_cast<MKL_Complex16*>(d2.Data()), reinterpret_cast<MKL_Complex16*>(x.Data()), reinterpret_cast<MKL_Complex16*>(y.Data()), reinterpret_cast<MKL_Complex16*>(work.Data()), &bs, &info);
02207 
02208         if(info != 0)
02209         {
02210                 if(info < 0)
02211                 {
02212                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
02213                         char s[500];
02214                         sprintf(s, "%i'th argument of ZGGGLM is an illegal value!", -info);
02215                         Utility::RunTimeError(s);
02216                 }
02217                 else if(info > 0)  
02218                 {
02219                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
02220                         char s[500];
02221                         sprintf(s, "Problem experienced when running ZGGGLM!");
02222                         Utility::RunTimeError(s);
02223                 }
02224         }
02225 }
02226 
02227 
02228 
02229 
02230 
02231 
02232 //=========================================================
02233 // Symmetric eigenproblem
02234 //=========================================================
02235 
02236 
02237 
02238 // Syevd or Heevd
02239 Vector<float> Lapack::Syevd(Matrix<float>& A)
02240 {
02241         if(A.Rows() != A.Columns())
02242         {
02243                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
02244                 Utility::RunTimeError("Matrix is not square!");
02245         }
02246         
02247         int info;
02248         Matrix<float> F = A.Clone();
02249         int n = F.Rows();
02250         Vector<float> w(n);
02251 
02252         int bs = 2*n+1;
02253         Vector<float> work(bs);
02254 
02255         int iwork_size = 1;
02256         Vector<int> iwork(iwork_size);
02257 
02258         SSYEVD("N", "U", &n, F.Data(), &n, w.Data(), work.Data(), &bs, iwork.Data(), &iwork_size, &info);
02259 
02260         if(info != 0)
02261         {
02262                 if(info < 0)
02263                 {
02264                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
02265                         char s[500];
02266                         sprintf(s, "%i'th argument of SSYEVD is an illegal value!", -info);
02267                         Utility::RunTimeError(s);
02268                 }
02269                 else if(info > 0)  
02270                 {
02271                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
02272                         char s[500];
02273                         sprintf(s, "The algorithm failed to converge; %i elements of an intermediate tridiagonal form did not converge to zero!", info);
02274                         Utility::RunTimeError(s);
02275                 }
02276         }
02277         return w;
02278 }
02279 
02280 
02281 Vector<float> Lapack::Syevd(Matrix<float>& A, Matrix<float>& EV)
02282 {
02283         if(A.Rows() != A.Columns())
02284         {
02285                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
02286                 Utility::RunTimeError("Matrix is not square!");
02287         }
02288         
02289         int info;
02290         EV = A.Clone();
02291         int n = EV.Rows();
02292         Vector<float> w(n);
02293 
02294         int k = 1;
02295         if(n>1)
02296         {
02297                 k = (int)(log((double)(n-1))/log(2.0))+1;
02298         }
02299         int bs = 3*n*n+(5+2*k)*n+1;
02300         Vector<float> work(bs);
02301 
02302         int iwork_size = 5*n+3;
02303         Vector<int> iwork(iwork_size);
02304 
02305         SSYEVD("V", "U", &n, EV.Data(), &n, w.Data(), work.Data(), &bs, iwork.Data(), &iwork_size, &info);
02306         
02307         if(info != 0)
02308         {
02309                 if(info < 0)
02310                 {
02311                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
02312                         char s[500];
02313                         sprintf(s, "%i'th argument of SSYEVD is an illegal value!", -info);
02314                         Utility::RunTimeError(s);
02315                 }
02316                 else if(info > 0)  
02317                 {
02318                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
02319                         char s[500];
02320                         sprintf(s, "The algorithm failed to converge; %i elements of an intermediate tridiagonal form did not converge to zero!", info);
02321                         Utility::RunTimeError(s);
02322                 }
02323         }
02324         return w;
02325 }
02326 
02327 
02328 Vector<double> Lapack::Syevd(Matrix<double>& A)
02329 {
02330         if(A.Rows() != A.Columns())
02331         {
02332                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
02333                 Utility::RunTimeError("Matrix is not square!");
02334         }
02335         
02336         int info;
02337         Matrix<double> F = A.Clone();
02338         int n = F.Rows();
02339         Vector<double> w(n);
02340 
02341         int bs = 2*n+1;
02342         Vector<double> work(bs);
02343 
02344         int iwork_size = 1;
02345         Vector<int> iwork(iwork_size);
02346 
02347         DSYEVD("N", "U", &n, F.Data(), &n, w.Data(), work.Data(), &bs, iwork.Data(), &iwork_size, &info);
02348 
02349         if(info != 0)
02350         {
02351                 if(info < 0)
02352                 {
02353                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
02354                         char s[500];
02355                         sprintf(s, "%i'th argument of DSYEVD is an illegal value!", -info);
02356                         Utility::RunTimeError(s);
02357                 }
02358                 else if(info > 0)  
02359                 {
02360                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
02361                         char s[500];
02362                         sprintf(s, "The algorithm failed to converge; %i elements of an intermediate tridiagonal form did not converge to zero!", info);
02363                         Utility::RunTimeError(s);
02364                 }
02365         }
02366         return w;
02367 }
02368 
02369 Vector<double> Lapack::Syevd(Matrix<double>& A, Matrix<double>& EV)
02370 {
02371         if(A.Rows() != A.Columns())
02372         {
02373                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
02374                 Utility::RunTimeError("Matrix is not square!");
02375         }
02376         
02377         int info;
02378         EV = A.Clone();
02379         int n = EV.Rows();
02380         Vector<double> w(n);
02381 
02382         int k = 1;
02383         if(n>1)
02384         {
02385                 k = (int)(log((double)(n-1))/log(2.0))+1;
02386         }
02387         int bs = 3*n*n+(5+2*k)*n+1;
02388         Vector<double> work(bs);
02389 
02390         int iwork_size = 5*n+3;
02391         Vector<int> iwork(iwork_size);
02392 
02393         DSYEVD("V", "U", &n, EV.Data(), &n, w.Data(), work.Data(), &bs, iwork.Data(), &iwork_size, &info);
02394         
02395         if(info != 0)
02396         {
02397                 if(info < 0)
02398                 {
02399                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
02400                         char s[500];
02401                         sprintf(s, "%i'th argument of DSYEVD is an illegal value!", -info);
02402                         Utility::RunTimeError(s);
02403                 }
02404                 else if(info > 0)  
02405                 {
02406                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
02407                         char s[500];
02408                         sprintf(s, "The algorithm failed to converge; %i elements of an intermediate tridiagonal form did not converge to zero!", info);
02409                         Utility::RunTimeError(s);
02410                 }
02411         }
02412         return w;
02413 }
02414 
02415 Vector<float> Lapack::Heevd(Matrix<ComplexFloat>& A)
02416 {
02417         if(A.Rows() != A.Columns())
02418         {
02419                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
02420                 Utility::RunTimeError("Matrix is not square!");
02421         }
02422         
02423         int info;
02424         Matrix<ComplexFloat> F = A.Clone();
02425         int n = F.Rows();
02426         Vector<float> w(n);
02427 
02428         int bs = n+1;
02429         Vector<ComplexFloat> work(bs);
02430 
02431         int lrwork = n;
02432         Vector<float> rwork(lrwork);
02433 
02434         int iwork_size = 1;
02435         Vector<int> iwork(iwork_size);
02436 
02437         CHEEVD("N", "U", &n, reinterpret_cast<MKL_Complex8*>(F.Data()), &n, w.Data(), reinterpret_cast<MKL_Complex8*>(work.Data()), &bs, rwork.Data(), &lrwork, iwork.Data(), &iwork_size, &info);
02438 
02439         if(info != 0)
02440         {
02441                 if(info < 0)
02442                 {
02443                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
02444                         char s[500];
02445                         sprintf(s, "%i'th argument of CHEEVD is an illegal value!", -info);
02446                         Utility::RunTimeError(s);
02447                 }
02448                 else if(info > 0)  
02449                 {
02450                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
02451                         char s[500];
02452                         sprintf(s, "The algorithm failed to converge; %i elements of an intermediate tridiagonal form did not converge to zero!", info);
02453                         Utility::RunTimeError(s);
02454                 }
02455         }
02456         return w;
02457 }
02458 
02459 Vector<float> Lapack::Heevd(Matrix<ComplexFloat>& A, Matrix<ComplexFloat>& EV)
02460 {
02461         if(A.Rows() != A.Columns())
02462         {
02463                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
02464                 Utility::RunTimeError("Matrix is not square!");
02465         }
02466         
02467         int info;
02468         EV = A.Clone();
02469         int n = EV.Rows();
02470         Vector<float> w(n);
02471 
02472         int k = 1;
02473         if(n>1)
02474         {
02475                 k = (int)(log((double)(n-1))/log(2.0))+1;
02476         }
02477         int bs = n*n+2*n;
02478         Vector<ComplexFloat> work(bs);
02479 
02480         int lrwork = 3*n*n+(4+2*k)*n+1;
02481         Vector<float> rwork(lrwork);
02482 
02483         int iwork_size = 5*n+3;
02484         Vector<int> iwork(iwork_size);
02485 
02486         CHEEVD("V", "U", &n, reinterpret_cast<MKL_Complex8*>(EV.Data()), &n, w.Data(), reinterpret_cast<MKL_Complex8*>(work.Data()), &bs, rwork.Data(), &lrwork, iwork.Data(), &iwork_size, &info);
02487         
02488         if(info != 0)
02489         {
02490                 if(info < 0)
02491                 {
02492                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
02493                         char s[500];
02494                         sprintf(s, "%i'th argument of CHEEVD is an illegal value!", -info);
02495                         Utility::RunTimeError(s);
02496                 }
02497                 else if(info > 0)  
02498                 {
02499                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
02500                         char s[500];
02501                         sprintf(s, "The algorithm failed to converge; %i elements of an intermediate tridiagonal form did not converge to zero!", info);
02502                         Utility::RunTimeError(s);
02503                 }
02504         }
02505         return w;
02506 }
02507 
02508 Vector<double> Lapack::Heevd(Matrix<ComplexDouble>& A)
02509 {
02510         if(A.Rows() != A.Columns())
02511         {
02512                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
02513                 Utility::RunTimeError("Matrix is not square!");
02514         }
02515         
02516         int info;
02517         Matrix<ComplexDouble> F = A.Clone();
02518         int n = F.Rows();
02519         Vector<double> w(n);
02520 
02521         int bs = n+1;
02522         Vector<ComplexDouble> work(bs);
02523 
02524         int lrwork = n;
02525         Vector<double> rwork(lrwork);
02526 
02527         int iwork_size = 1;
02528         Vector<int> iwork(iwork_size);
02529 
02530         ZHEEVD("N", "U", &n, reinterpret_cast<MKL_Complex16*>(F.Data()), &n, w.Data(), reinterpret_cast<MKL_Complex16*>(work.Data()), &bs, rwork.Data(), &lrwork, iwork.Data(), &iwork_size, &info);
02531 
02532         if(info != 0)
02533         {
02534                 if(info < 0)
02535                 {
02536                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
02537                         char s[500];
02538                         sprintf(s, "%i'th argument of ZHEEVD is an illegal value!", -info);
02539                         Utility::RunTimeError(s);
02540                 }
02541                 else if(info > 0)  
02542                 {
02543                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
02544                         char s[500];
02545                         sprintf(s, "The algorithm failed to converge; %i elements of an intermediate tridiagonal form did not converge to zero!", info);
02546                         Utility::RunTimeError(s);
02547                 }
02548         }
02549         return w;
02550 }
02551 
02552 Vector<double> Lapack::Heevd(Matrix<ComplexDouble>& A, Matrix<ComplexDouble>& EV)
02553 {
02554         if(A.Rows() != A.Columns())
02555         {
02556                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
02557                 Utility::RunTimeError("Matrix is not square!");
02558         }
02559         
02560         int info;
02561         EV = A.Clone();
02562         int n = EV.Rows();
02563         Vector<double> w(n);
02564 
02565         int k = 1;
02566         if(n>1)
02567         {
02568                 k = (int)(log((double)(n-1))/log(2.0))+1;
02569         }
02570         int bs = n*n+2*n;
02571         Vector<ComplexDouble> work(bs);
02572 
02573         int lrwork = 3*n*n+(4+2*k)*n+1;
02574         Vector<double> rwork(lrwork);
02575 
02576         int iwork_size = 5*n+3;
02577         Vector<int> iwork(iwork_size);
02578 
02579         ZHEEVD("V", "U", &n, reinterpret_cast<MKL_Complex16*>(EV.Data()), &n, w.Data(), reinterpret_cast<MKL_Complex16*>(work.Data()), &bs, rwork.Data(), &lrwork, iwork.Data(), &iwork_size, &info);
02580         
02581         if(info != 0)
02582         {
02583                 if(info < 0)
02584                 {
02585                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
02586                         char s[500];
02587                         sprintf(s, "%i'th argument of ZHEEVD is an illegal value!", -info);
02588                         Utility::RunTimeError(s);
02589                 }
02590                 else if(info > 0)  
02591                 {
02592                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
02593                         char s[500];
02594                         sprintf(s, "The algorithm failed to converge; %i elements of an intermediate tridiagonal form did not converge to zero!", info);
02595                         Utility::RunTimeError(s);
02596                 }
02597         }
02598         return w;
02599 }
02600 
02601 
02602 
02603 
02604 
02605 
02606 
02607 
02608 
02609 
02610 // Syevr
02611 Vector<float> Lapack::Syevr(Matrix<float>& A)
02612 {
02613         if(A.Rows() != A.Columns())
02614         {
02615                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
02616                 Utility::RunTimeError("Matrix is not square!");
02617         }
02618         
02619         int info;
02620         float abstol = SLAMCH("S");
02621 
02622         int m;
02623         float vl, vu;
02624         int il, iu;
02625         
02626         Matrix<float> F = A.Clone();
02627         int n = F.Rows();
02628         Vector<float> w(n);
02629 
02630         Vector<int> isuppz(2*n);
02631 
02632         int bs = (64+6)*n;
02633         Vector<float> work(bs);
02634 
02635         int iwork_size = 10*n;
02636         Vector<int> iwork(iwork_size);
02637 
02638         SSYEVR("N", "A", "U", &n, F.Data(), &n, &vl, &vu, &il, &iu, &abstol, &m,  w.Data(), 0, &n, isuppz.Data(), work.Data(), &bs, iwork.Data(), &iwork_size, &info);
02639         
02640         if(info != 0)
02641         {
02642                 if(info < 0)
02643                 {
02644                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
02645                         char s[500];
02646                         sprintf(s, "%i'th argument of SSYEVR is an illegal value!", -info);
02647                         Utility::RunTimeError(s);
02648                 }
02649                 else if(info > 0)  
02650                 {
02651                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
02652                         char s[500];
02653                         sprintf(s, "Internal error occured when running SSYEVR");
02654                         Utility::RunTimeError(s);
02655                 }
02656         }
02657         return w;
02658 }
02659 
02660 
02661 Vector<float> Lapack::Syevr(Matrix<float>& A, int il, int iu)
02662 {
02663         if(A.Rows() != A.Columns())
02664         {
02665                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
02666                 Utility::RunTimeError("Matrix is not square!");
02667         }
02668 
02669         if(il > iu || il > A.Rows() || iu > A.Rows())
02670         {
02671                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
02672                 Utility::RunTimeError("The choice of parameters (il and iu) for selecting a subset of eigenvalues is not correct!");
02673         }
02674         
02675         int info;
02676         float abstol = SLAMCH("S");
02677 
02678         int m;
02679         float vl, vu;
02680         
02681         Matrix<float> F = A.Clone();
02682         int n = F.Rows();
02683         Vector<float> w(iu-il+1);
02684 
02685         Vector<int> isuppz(2*n);
02686 
02687         int bs = (64+6)*n;
02688         Vector<float> work(bs);
02689 
02690         int iwork_size = 10*n;
02691         Vector<int> iwork(iwork_size);
02692 
02693         SSYEVR("N", "I", "U", &n, F.Data(), &n, &vl, &vu, &il, &iu, &abstol, &m,  w.Data(), 0, &n, isuppz.Data(), work.Data(), &bs, iwork.Data(), &iwork_size, &info);
02694         
02695         if(info != 0)
02696         {
02697                 if(info < 0)
02698                 {
02699                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
02700                         char s[500];
02701                         sprintf(s, "%i'th argument of SSYEVR is an illegal value!", -info);
02702                         Utility::RunTimeError(s);
02703                 }
02704                 else if(info > 0)  
02705                 {
02706                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
02707                         char s[500];
02708                         sprintf(s, "Internal error occured when running SSYEVR");
02709                         Utility::RunTimeError(s);
02710                 }
02711         }
02712 
02713         return w;
02714 }
02715 
02716 Vector<float> Lapack::Syevr(Matrix<float>& A, Matrix<float>& EV)
02717 {
02718         if(A.Rows() != A.Columns())
02719         {
02720                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
02721                 Utility::RunTimeError("Matrix is not square!");
02722         }
02723         
02724         int info;
02725         float abstol = SLAMCH("S");
02726 
02727         int m;
02728         float vl, vu;
02729         int il, iu;
02730         
02731         Matrix<float> F = A.Clone();
02732         int n = F.Rows();
02733         Vector<float> w(n);
02734         EV = Matrix<float>(n,n);
02735 
02736         Vector<int> isuppz(2*n);
02737 
02738         int bs = (64+6)*n;
02739         Vector<float> work(bs);
02740 
02741         int iwork_size = 10*n;
02742         Vector<int> iwork(iwork_size);
02743 
02744         SSYEVR("V", "A", "U", &n, F.Data(), &n, &vl, &vu, &il, &iu, &abstol, &m,  w.Data(), EV.Data(), &n, isuppz.Data(), work.Data(), &bs, iwork.Data(), &iwork_size, &info);
02745         
02746         if(info != 0)
02747         {
02748                 if(info < 0)
02749                 {
02750                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
02751                         char s[500];
02752                         sprintf(s, "%i'th argument of SSYEVR is an illegal value!", -info);
02753                         Utility::RunTimeError(s);
02754                 }
02755                 else if(info > 0)  
02756                 {
02757                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
02758                         char s[500];
02759                         sprintf(s, "Internal error occured when running SSYEVR");
02760                         Utility::RunTimeError(s);
02761                 }
02762         }
02763         return w;
02764 }
02765 
02766 Vector<float> Lapack::Syevr(Matrix<float>& A, Matrix<float>& EV, int il, int iu)
02767 {
02768         if(A.Rows() != A.Columns())
02769         {
02770                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
02771                 Utility::RunTimeError("Matrix is not square!");
02772         }
02773 
02774         if(il > iu || il > A.Rows() || iu > A.Rows())
02775         {
02776                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
02777                 Utility::RunTimeError("The choice of parameters (il and iu) for selecting a subset of eigenvalues is not correct!");
02778         }
02779         
02780         int info;
02781         float abstol = SLAMCH("S");
02782 
02783         int m;
02784         float vl, vu;
02785         
02786         Matrix<float> F = A.Clone();
02787         int n = F.Rows();
02788         Vector<float> w(iu-il+1);
02789         EV = Matrix<float>(n,iu-il+1);
02790 
02791         Vector<int> isuppz(2*n);
02792 
02793         int bs = (64+6)*n;
02794         Vector<float> work(bs);
02795 
02796         int iwork_size = 10*n;
02797         Vector<int> iwork(iwork_size);
02798 
02799         SSYEVR("V", "I", "U", &n, F.Data(), &n, &vl, &vu, &il, &iu, &abstol, &m,  w.Data(), EV.Data(), &n, isuppz.Data(), work.Data(), &bs, iwork.Data(), &iwork_size, &info);
02800         
02801         if(info != 0)
02802         {
02803                 if(info < 0)
02804                 {
02805                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
02806                         char s[500];
02807                         sprintf(s, "%i'th argument of SSYEVR is an illegal value!", -info);
02808                         Utility::RunTimeError(s);
02809                 }
02810                 else if(info > 0)  
02811                 {
02812                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
02813                         char s[500];
02814                         sprintf(s, "Internal error occured when running SSYEVR");
02815                         Utility::RunTimeError(s);
02816                 }
02817         }
02818 
02819         return w;
02820 }
02821 
02822 
02823 
02824 
02825 
02827 
02828 Vector<double> Lapack::Syevr(Matrix<double>& A)
02829 {
02830         if(A.Rows() != A.Columns())
02831         {
02832                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
02833                 Utility::RunTimeError("Matrix is not square!");
02834         }
02835         
02836         int info;
02837         double abstol = SLAMCH("S");
02838 
02839         int m;
02840         double vl, vu;
02841         int il, iu;
02842         
02843         Matrix<double> F = A.Clone();
02844         int n = F.Rows();
02845         Vector<double> w(n);
02846 
02847         Vector<int> isuppz(2*n);
02848 
02849         int bs = (64+6)*n;
02850         Vector<double> work(bs);
02851 
02852         int iwork_size = 10*n;
02853         Vector<int> iwork(iwork_size);
02854 
02855         DSYEVR("N", "A", "U", &n, F.Data(), &n, &vl, &vu, &il, &iu, &abstol, &m,  w.Data(), 0, &n, isuppz.Data(), work.Data(), &bs, iwork.Data(), &iwork_size, &info);
02856         
02857         if(info != 0)
02858         {
02859                 if(info < 0)
02860                 {
02861                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
02862                         char s[500];
02863                         sprintf(s, "%i'th argument of DSYEVR is an illegal value!", -info);
02864                         Utility::RunTimeError(s);
02865                 }
02866                 else if(info > 0)  
02867                 {
02868                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
02869                         char s[500];
02870                         sprintf(s, "Internal error occured when running DSYEVR");
02871                         Utility::RunTimeError(s);
02872                 }
02873         }
02874         return w;
02875 }
02876 
02877 
02878 Vector<double> Lapack::Syevr(Matrix<double>& A, int il, int iu)
02879 {
02880         if(A.Rows() != A.Columns())
02881         {
02882                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
02883                 Utility::RunTimeError("Matrix is not square!");
02884         }
02885 
02886         if(il > iu || il > A.Rows() || iu > A.Rows())
02887         {
02888                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
02889                 Utility::RunTimeError("The choice of parameters (il and iu) for selecting a subset of eigenvalues is not correct!");
02890         }
02891         
02892         int info;
02893         double abstol = SLAMCH("S");
02894 
02895         int m;
02896         double vl, vu;
02897         
02898         Matrix<double> F = A.Clone();
02899         int n = F.Rows();
02900         Vector<double> w(iu-il+1);
02901 
02902         Vector<int> isuppz(2*n);
02903 
02904         int bs = (64+6)*n;
02905         Vector<double> work(bs);
02906 
02907         int iwork_size = 10*n;
02908         Vector<int> iwork(iwork_size);
02909 
02910         DSYEVR("N", "I", "U", &n, F.Data(), &n, &vl, &vu, &il, &iu, &abstol, &m,  w.Data(), 0, &n, isuppz.Data(), work.Data(), &bs, iwork.Data(), &iwork_size, &info);
02911         
02912         if(info != 0)
02913         {
02914                 if(info < 0)
02915                 {
02916                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
02917                         char s[500];
02918                         sprintf(s, "%i'th argument of DSYEVR is an illegal value!", -info);
02919                         Utility::RunTimeError(s);
02920                 }
02921                 else if(info > 0)  
02922                 {
02923                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
02924                         char s[500];
02925                         sprintf(s, "Internal error occured when running DSYEVR");
02926                         Utility::RunTimeError(s);
02927                 }
02928         }
02929 
02930         return w;
02931 }
02932 
02933 Vector<double> Lapack::Syevr(Matrix<double>& A, Matrix<double>& EV)
02934 {
02935         if(A.Rows() != A.Columns())
02936         {
02937                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
02938                 Utility::RunTimeError("Matrix is not square!");
02939         }
02940         
02941         int info;
02942         double abstol = SLAMCH("S");
02943 
02944         int m;
02945         double vl, vu;
02946         int il, iu;
02947         
02948         Matrix<double> F = A.Clone();
02949         int n = F.Rows();
02950         Vector<double> w(n);
02951         EV = Matrix<double>(n,n);
02952 
02953         Vector<int> isuppz(2*n);
02954 
02955         int bs = (64+6)*n;
02956         Vector<double> work(bs);
02957 
02958         int iwork_size = 10*n;
02959         Vector<int> iwork(iwork_size);
02960 
02961         DSYEVR("V", "A", "U", &n, F.Data(), &n, &vl, &vu, &il, &iu, &abstol, &m,  w.Data(), EV.Data(), &n, isuppz.Data(), work.Data(), &bs, iwork.Data(), &iwork_size, &info);
02962         
02963         if(info != 0)
02964         {
02965                 if(info < 0)
02966                 {
02967                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
02968                         char s[500];
02969                         sprintf(s, "%i'th argument of DSYEVR is an illegal value!", -info);
02970                         Utility::RunTimeError(s);
02971                 }
02972                 else if(info > 0)  
02973                 {
02974                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
02975                         char s[500];
02976                         sprintf(s, "Internal error occured when running DSYEVR");
02977                         Utility::RunTimeError(s);
02978                 }
02979         }
02980         return w;
02981 }
02982 
02983 Vector<double> Lapack::Syevr(Matrix<double>& A, Matrix<double>& EV, int il, int iu)
02984 {
02985         if(A.Rows() != A.Columns())
02986         {
02987                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
02988                 Utility::RunTimeError("Matrix is not square!");
02989         }
02990 
02991         if(il > iu || il > A.Rows() || iu > A.Rows())
02992         {
02993                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
02994                 Utility::RunTimeError("The choice of parameters (il and iu) for selecting a subset of eigenvalues is not correct!");
02995         }
02996         
02997         int info;
02998         double abstol = SLAMCH("S");
02999 
03000         int m;
03001         double vl, vu;
03002         
03003         Matrix<double> F = A.Clone();
03004         int n = F.Rows();
03005         Vector<double> w(iu-il+1);
03006         EV = Matrix<double>(n,iu-il+1);
03007 
03008         Vector<int> isuppz(2*n);
03009 
03010         int bs = (64+6)*n;
03011         Vector<double> work(bs);
03012 
03013         int iwork_size = 10*n;
03014         Vector<int> iwork(iwork_size);
03015 
03016         DSYEVR("V", "I", "U", &n, F.Data(), &n, &vl, &vu, &il, &iu, &abstol, &m,  w.Data(), EV.Data(), &n, isuppz.Data(), work.Data(), &bs, iwork.Data(), &iwork_size, &info);
03017         
03018         if(info != 0)
03019         {
03020                 if(info < 0)
03021                 {
03022                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
03023                         char s[500];
03024                         sprintf(s, "%i'th argument of DSYEVR is an illegal value!", -info);
03025                         Utility::RunTimeError(s);
03026                 }
03027                 else if(info > 0)  
03028                 {
03029                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
03030                         char s[500];
03031                         sprintf(s, "Internal error occured when running DSYEVR");
03032                         Utility::RunTimeError(s);
03033                 }
03034         }
03035 
03036         return w;
03037 }
03038 
03039 
03040 
03041 
03043 
03044 Vector<float> Lapack::Heevr(Matrix<ComplexFloat>& A)
03045 {
03046         if(A.Rows() != A.Columns())
03047         {
03048                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
03049                 Utility::RunTimeError("Matrix is not square!");
03050         }
03051         
03052         int info;
03053         float abstol = SLAMCH("S");
03054 
03055         int m;
03056         float vl, vu;
03057         int il, iu;
03058         
03059         Matrix<ComplexFloat> F = A.Clone();
03060         int n = F.Rows();
03061         Vector<float> w(n);
03062 
03063         Vector<int> isuppz(2*n);
03064 
03065         int bs = (64+1)*n;
03066         Vector<ComplexFloat> work(bs);
03067 
03068         int rwork_size = 24*n;
03069         Vector<float> rwork(rwork_size);
03070 
03071         int iwork_size = 10*n;
03072         Vector<int> iwork(iwork_size);
03073 
03074         CHEEVR("N", "A", "U", &n, reinterpret_cast<MKL_Complex8*>(F.Data()), &n, &vl, &vu, &il, &iu, &abstol, &m,  w.Data(), 0, &n, isuppz.Data(), reinterpret_cast<MKL_Complex8*>(work.Data()), &bs, rwork.Data(), &rwork_size, iwork.Data(), &iwork_size, &info);
03075         
03076         if(info != 0)
03077         {
03078                 if(info < 0)
03079                 {
03080                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
03081                         char s[500];
03082                         sprintf(s, "%i'th argument of CHEEVR is an illegal value!", -info);
03083                         Utility::RunTimeError(s);
03084                 }
03085                 else if(info > 0)  
03086                 {
03087                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
03088                         char s[500];
03089                         sprintf(s, "Internal error occured when running CHEEVR");
03090                         Utility::RunTimeError(s);
03091                 }
03092         }
03093         return w;
03094 }
03095 
03096 
03097 Vector<float> Lapack::Heevr(Matrix<ComplexFloat>& A, int il, int iu)
03098 {
03099         if(A.Rows() != A.Columns())
03100         {
03101                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
03102                 Utility::RunTimeError("Matrix is not square!");
03103         }
03104 
03105         if(il > iu || il > A.Rows() || iu > A.Rows())
03106         {
03107                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
03108                 Utility::RunTimeError("The choice of parameters (il and iu) for selecting a subset of eigenvalues is not correct!");
03109         }
03110         
03111         int info;
03112         float abstol = SLAMCH("S");
03113 
03114         int m;
03115         float vl, vu;
03116         
03117         Matrix<ComplexFloat> F = A.Clone();
03118         int n = F.Rows();
03119         Vector<float> w(iu-il+1);
03120 
03121         Vector<int> isuppz(2*n);
03122 
03123         int bs = (64+6)*n;
03124         Vector<ComplexFloat> work(bs);
03125 
03126         int rwork_size = 24*n;
03127         Vector<float> rwork(rwork_size);
03128         int iwork_size = 10*n;
03129         Vector<int> iwork(iwork_size);
03130 
03131         CHEEVR("N", "I", "U", &n, reinterpret_cast<MKL_Complex8*>(F.Data()), &n, &vl, &vu, &il, &iu, &abstol, &m,  w.Data(), 0, &n, isuppz.Data(), reinterpret_cast<MKL_Complex8*>(work.Data()), &bs, rwork.Data(), &rwork_size, iwork.Data(), &iwork_size, &info);
03132         
03133         if(info != 0)
03134         {
03135                 if(info < 0)
03136                 {
03137                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
03138                         char s[500];
03139                         sprintf(s, "%i'th argument of CHEEVR is an illegal value!", -info);
03140                         Utility::RunTimeError(s);
03141                 }
03142                 else if(info > 0)  
03143                 {
03144                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
03145                         char s[500];
03146                         sprintf(s, "Internal error occured when running CHEEVR");
03147                         Utility::RunTimeError(s);
03148                 }
03149         }
03150 
03151         return w;
03152 }
03153 
03154 Vector<float> Lapack::Heevr(Matrix<ComplexFloat>& A, Matrix<ComplexFloat>& EV)
03155 {
03156         if(A.Rows() != A.Columns())
03157         {
03158                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
03159                 Utility::RunTimeError("Matrix is not square!");
03160         }
03161         
03162         int info;
03163         float abstol = SLAMCH("S");
03164 
03165         int m;
03166         float vl, vu;
03167         int il, iu;
03168         
03169         Matrix<ComplexFloat> F = A.Clone();
03170         int n = F.Rows();
03171         Vector<float> w(n);
03172         EV = Matrix<ComplexFloat>(n,n);
03173 
03174         Vector<int> isuppz(2*n);
03175 
03176         int bs = (64+6)*n;
03177         Vector<ComplexFloat> work(bs);
03178 
03179         int rwork_size = 24*n;
03180         Vector<float> rwork(rwork_size);
03181         int iwork_size = 10*n;
03182         Vector<int> iwork(iwork_size);
03183 
03184         CHEEVR("V", "A", "U", &n, reinterpret_cast<MKL_Complex8*>(F.Data()), &n, &vl, &vu, &il, &iu, &abstol, &m,  w.Data(), reinterpret_cast<MKL_Complex8*>(EV.Data()), &n, isuppz.Data(), reinterpret_cast<MKL_Complex8*>(work.Data()), &bs, rwork.Data(), &rwork_size, iwork.Data(), &iwork_size, &info);
03185         
03186         if(info != 0)
03187         {
03188                 if(info < 0)
03189                 {
03190                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
03191                         char s[500];
03192                         sprintf(s, "%i'th argument of CHEEVR is an illegal value!", -info);
03193                         Utility::RunTimeError(s);
03194                 }
03195                 else if(info > 0)  
03196                 {
03197                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
03198                         char s[500];
03199                         sprintf(s, "Internal error occured when running CHEEVR");
03200                         Utility::RunTimeError(s);
03201                 }
03202         }
03203         return w;
03204 }
03205 
03206 Vector<float> Lapack::Heevr(Matrix<ComplexFloat>& A, Matrix<ComplexFloat>& EV, int il, int iu)
03207 {
03208         if(A.Rows() != A.Columns())
03209         {
03210                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
03211                 Utility::RunTimeError("Matrix is not square!");
03212         }
03213 
03214         if(il > iu || il > A.Rows() || iu > A.Rows())
03215         {
03216                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
03217                 Utility::RunTimeError("The choice of parameters (il and iu) for selecting a subset of eigenvalues is not correct!");
03218         }
03219         
03220         int info;
03221         float abstol = SLAMCH("S");
03222 
03223         int m;
03224         float vl, vu;
03225         
03226         Matrix<ComplexFloat> F = A.Clone();
03227         int n = F.Rows();
03228         Vector<float> w(iu-il+1);
03229         EV = Matrix<ComplexFloat>(n,iu-il+1);
03230 
03231         Vector<int> isuppz(2*n);
03232 
03233         int bs = (64+6)*n;
03234         Vector<ComplexFloat> work(bs);
03235 
03236         int rwork_size = 24*n;
03237         Vector<float> rwork(rwork_size);
03238         int iwork_size = 10*n;
03239         Vector<int> iwork(iwork_size);
03240 
03241         CHEEVR("V", "I", "U", &n, reinterpret_cast<MKL_Complex8*>(F.Data()), &n, &vl, &vu, &il, &iu, &abstol, &m,  w.Data(), reinterpret_cast<MKL_Complex8*>(EV.Data()), &n, isuppz.Data(), reinterpret_cast<MKL_Complex8*>(work.Data()), &bs, rwork.Data(), &rwork_size, iwork.Data(), &iwork_size, &info);
03242         
03243         if(info != 0)
03244         {
03245                 if(info < 0)
03246                 {
03247                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
03248                         char s[500];
03249                         sprintf(s, "%i'th argument of CHEEVR is an illegal value!", -info);
03250                         Utility::RunTimeError(s);
03251                 }
03252                 else if(info > 0)  
03253                 {
03254                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
03255                         char s[500];
03256                         sprintf(s, "Internal error occured when running CHEEVR");
03257                         Utility::RunTimeError(s);
03258                 }
03259         }
03260 
03261         return w;
03262 }
03263 
03264 
03265 
03267 
03268 Vector<double> Lapack::Heevr(Matrix<ComplexDouble>& A)
03269 {
03270         if(A.Rows() != A.Columns())
03271         {
03272                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
03273                 Utility::RunTimeError("Matrix is not square!");
03274         }
03275         
03276         int info;
03277         double abstol = SLAMCH("S");
03278 
03279         int m;
03280         double vl, vu;
03281         int il, iu;
03282         
03283         Matrix<ComplexDouble> F = A.Clone();
03284         int n = F.Rows();
03285         Vector<double> w(n);
03286 
03287         Vector<int> isuppz(2*n);
03288 
03289         int bs = (64+1)*n;
03290         Vector<ComplexDouble> work(bs);
03291 
03292         int rwork_size = 24*n;
03293         Vector<double> rwork(rwork_size);
03294 
03295         int iwork_size = 10*n;
03296         Vector<int> iwork(iwork_size);
03297 
03298         ZHEEVR("N", "A", "U", &n, reinterpret_cast<MKL_Complex16*>(F.Data()), &n, &vl, &vu, &il, &iu, &abstol, &m,  w.Data(), 0, &n, isuppz.Data(), reinterpret_cast<MKL_Complex16*>(work.Data()), &bs, rwork.Data(), &rwork_size, iwork.Data(), &iwork_size, &info);
03299         
03300         if(info != 0)
03301         {
03302                 if(info < 0)
03303                 {
03304                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
03305                         char s[500];
03306                         sprintf(s, "%i'th argument of ZHEEVR is an illegal value!", -info);
03307                         Utility::RunTimeError(s);
03308                 }
03309                 else if(info > 0)  
03310                 {
03311                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
03312                         char s[500];
03313                         sprintf(s, "Internal error occured when running ZHEEVR");
03314                         Utility::RunTimeError(s);
03315                 }
03316         }
03317         return w;
03318 }
03319 
03320 
03321 Vector<double> Lapack::Heevr(Matrix<ComplexDouble>& A, int il, int iu)
03322 {
03323         if(A.Rows() != A.Columns())
03324         {
03325                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
03326                 Utility::RunTimeError("Matrix is not square!");
03327         }
03328 
03329         if(il > iu || il > A.Rows() || iu > A.Rows())
03330         {
03331                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
03332                 Utility::RunTimeError("The choice of parameters (il and iu) for selecting a subset of eigenvalues is not correct!");
03333         }
03334         
03335         int info;
03336         double abstol = SLAMCH("S");
03337 
03338         int m;
03339         double vl, vu;
03340         
03341         Matrix<ComplexDouble> F = A.Clone();
03342         int n = F.Rows();
03343         Vector<double> w(iu-il+1);
03344 
03345         Vector<int> isuppz(2*n);
03346 
03347         int bs = (64+6)*n;
03348         Vector<ComplexDouble> work(bs);
03349 
03350         int rwork_size = 24*n;
03351         Vector<double> rwork(rwork_size);
03352         int iwork_size = 10*n;
03353         Vector<int> iwork(iwork_size);
03354 
03355         ZHEEVR("N", "I", "U", &n, reinterpret_cast<MKL_Complex16*>(F.Data()), &n, &vl, &vu, &il, &iu, &abstol, &m,  w.Data(), 0, &n, isuppz.Data(), reinterpret_cast<MKL_Complex16*>(work.Data()), &bs, rwork.Data(), &rwork_size, iwork.Data(), &iwork_size, &info);
03356         
03357         if(info != 0)
03358         {
03359                 if(info < 0)
03360                 {
03361                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
03362                         char s[500];
03363                         sprintf(s, "%i'th argument of ZHEEVR is an illegal value!", -info);
03364                         Utility::RunTimeError(s);
03365                 }
03366                 else if(info > 0)  
03367                 {
03368                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
03369                         char s[500];
03370                         sprintf(s, "Internal error occured when running ZHEEVR");
03371                         Utility::RunTimeError(s);
03372                 }
03373         }
03374 
03375         return w;
03376 }
03377 
03378 Vector<double> Lapack::Heevr(Matrix<ComplexDouble>& A, Matrix<ComplexDouble>& EV)
03379 {
03380         if(A.Rows() != A.Columns())
03381         {
03382                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
03383                 Utility::RunTimeError("Matrix is not square!");
03384         }
03385         
03386         int info;
03387         double abstol = SLAMCH("S");
03388 
03389         int m;
03390         double vl, vu;
03391         int il, iu;
03392         
03393         Matrix<ComplexDouble> F = A.Clone();
03394         int n = F.Rows();
03395         Vector<double> w(n);
03396         EV = Matrix<ComplexDouble>(n,n);
03397 
03398         Vector<int> isuppz(2*n);
03399 
03400         int bs = (64+6)*n;
03401         Vector<ComplexDouble> work(bs);
03402 
03403         int rwork_size = 24*n;
03404         Vector<double> rwork(rwork_size);
03405         int iwork_size = 10*n;
03406         Vector<int> iwork(iwork_size);
03407 
03408         ZHEEVR("V", "A", "U", &n, reinterpret_cast<MKL_Complex16*>(F.Data()), &n, &vl, &vu, &il, &iu, &abstol, &m,  w.Data(), reinterpret_cast<MKL_Complex16*>(EV.Data()), &n, isuppz.Data(), reinterpret_cast<MKL_Complex16*>(work.Data()), &bs, rwork.Data(), &rwork_size, iwork.Data(), &iwork_size, &info);
03409         
03410         if(info != 0)
03411         {
03412                 if(info < 0)
03413                 {
03414                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
03415                         char s[500];
03416                         sprintf(s, "%i'th argument of ZHEEVR is an illegal value!", -info);
03417                         Utility::RunTimeError(s);
03418                 }
03419                 else if(info > 0)  
03420                 {
03421                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
03422                         char s[500];
03423                         sprintf(s, "Internal error occured when running ZHEEVR");
03424                         Utility::RunTimeError(s);
03425                 }
03426         }
03427         return w;
03428 }
03429 
03430 Vector<double> Lapack::Heevr(Matrix<ComplexDouble>& A, Matrix<ComplexDouble>& EV, int il, int iu)
03431 {
03432         if(A.Rows() != A.Columns())
03433         {
03434                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
03435                 Utility::RunTimeError("Matrix is not square!");
03436         }
03437 
03438         if(il > iu || il > A.Rows() || iu > A.Rows())
03439         {
03440                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
03441                 Utility::RunTimeError("The choice of parameters (il and iu) for selecting a subset of eigenvalues is not correct!");
03442         }
03443         
03444         int info;
03445         double abstol = SLAMCH("S");
03446 
03447         int m;
03448         double vl, vu;
03449         
03450         Matrix<ComplexDouble> F = A.Clone();
03451         int n = F.Rows();
03452         Vector<double> w(iu-il+1);
03453         EV = Matrix<ComplexDouble>(n,iu-il+1);
03454 
03455         Vector<int> isuppz(2*n);
03456 
03457         int bs = (64+6)*n;
03458         Vector<ComplexDouble> work(bs);
03459 
03460         int rwork_size = 24*n;
03461         Vector<double> rwork(rwork_size);
03462         int iwork_size = 10*n;
03463         Vector<int> iwork(iwork_size);
03464 
03465         ZHEEVR("V", "I", "U", &n, reinterpret_cast<MKL_Complex16*>(F.Data()), &n, &vl, &vu, &il, &iu, &abstol, &m,  w.Data(), reinterpret_cast<MKL_Complex16*>(EV.Data()), &n, isuppz.Data(), reinterpret_cast<MKL_Complex16*>(work.Data()), &bs, rwork.Data(), &rwork_size, iwork.Data(), &iwork_size, &info);
03466         
03467         if(info != 0)
03468         {
03469                 if(info < 0)
03470                 {
03471                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
03472                         char s[500];
03473                         sprintf(s, "%i'th argument of ZHEEVR is an illegal value!", -info);
03474                         Utility::RunTimeError(s);
03475                 }
03476                 else if(info > 0)  
03477                 {
03478                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
03479                         char s[500];
03480                         sprintf(s, "Internal error occured when running ZHEEVR");
03481                         Utility::RunTimeError(s);
03482                 }
03483         }
03484 
03485         return w;
03486 }
03487 
03488 
03489 
03490 
03491 
03492 
03493 
03494 
03495 
03496 
03497 
03498 //=========================================================
03499 // Non-Symmetric eigenvalue problem
03500 //=========================================================
03501 
03502 
03503 // Gees
03504 Vector<ComplexFloat> Lapack::Gees(Matrix<float>& A, Matrix<float>& T)
03505 {
03506         if(A.Rows() != A.Columns())
03507         {
03508                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
03509                 Utility::RunTimeError("Matrix is not square!");
03510         }
03511         
03512         int info;
03513         T = A.Clone();
03514         int n = T.Rows();
03515         Vector<float> wr(n);
03516         Vector<float> wi(n);
03517 
03518         int bs = 3*n;
03519         Vector<float> work(bs);
03520         int sdim;
03521         
03522         SGEES("N", "N", 0, &n, T.Data(), &n, &sdim, wr.Data(), wi.Data(), 0, &n, work.Data(), &bs, 0, &info);
03523 
03524         Vector<ComplexFloat> w(n);
03525         for(int i=0;i<n;i++)
03526         {
03527                 w[i] = ComplexFloat(wr[i],wi[i]);
03528         }
03529 
03530         if(info != 0)
03531         {
03532                 if(info < 0)
03533                 {
03534                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
03535                         char s[500];
03536                         sprintf(s, "%i'th argument of SGEES is an illegal value!", -info);
03537                         Utility::RunTimeError(s);
03538                 }
03539                 else if(info > 0)  
03540                 {
03541                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
03542                         char s[500];
03543                         sprintf(s, "The QR algorithm failed to compute all the eigenvalues!");
03544                         Utility::RunTimeError(s);
03545                 }
03546         }
03547         return w;
03548 }
03549 
03550 Vector<ComplexFloat> Lapack::Gees(Matrix<float>& A, Matrix<float>& T, Matrix<float>& Z)
03551 {
03552         if(A.Rows() != A.Columns())
03553         {
03554                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
03555                 Utility::RunTimeError("Matrix is not square!");
03556         }
03557         
03558         int info;
03559         T = A.Clone();
03560         int n = T.Rows();
03561         Vector<float> wr(n);
03562         Vector<float> wi(n);
03563 
03564         Z = Matrix<float>(n,n);
03565 
03566         int bs = 3*n;
03567         Vector<float> work(bs);
03568         int sdim;
03569 
03570         SGEES("V", "N", 0, &n, T.Data(), &n, &sdim, wr.Data(), wi.Data(), Z.Data(), &n, work.Data(), &bs, 0, &info);
03571 
03572         Vector<ComplexFloat> w(n);
03573         for(int i=0;i<n;i++)
03574         {
03575                 w[i] = ComplexFloat(wr[i],wi[i]);
03576         }
03577 
03578         if(info != 0)
03579         {
03580                 if(info < 0)
03581                 {
03582                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
03583                         char s[500];
03584                         sprintf(s, "%i'th argument of SGEES is an illegal value!", -info);
03585                         Utility::RunTimeError(s);
03586                 }
03587                 else if(info > 0)  
03588                 {
03589                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
03590                         char s[500];
03591                         sprintf(s, "The QR algorithm failed to compute all the eigenvalues!");
03592                         Utility::RunTimeError(s);
03593                 }
03594         }
03595         return w;
03596 }
03597 
03598 Vector<ComplexDouble> Lapack::Gees(Matrix<double>& A, Matrix<double>& T)
03599 {
03600         if(A.Rows() != A.Columns())
03601         {
03602                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
03603                 Utility::RunTimeError("Matrix is not square!");
03604         }
03605         
03606         int info;
03607         T = A.Clone();
03608         int n = T.Rows();
03609         Vector<double> wr(n);
03610         Vector<double> wi(n);
03611 
03612         int bs = 3*n;
03613         Vector<double> work(bs);
03614         int sdim;
03615         
03616         DGEES("N", "N", 0, &n, T.Data(), &n, &sdim, wr.Data(), wi.Data(), 0, &n, work.Data(), &bs, 0, &info);
03617 
03618         Vector<ComplexDouble> w(n);
03619         for(int i=0;i<n;i++)
03620         {
03621                 w[i] = ComplexDouble(wr[i],wi[i]);
03622         }
03623 
03624         if(info != 0)
03625         {
03626                 if(info < 0)
03627                 {
03628                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
03629                         char s[500];
03630                         sprintf(s, "%i'th argument of DGEES is an illegal value!", -info);
03631                         Utility::RunTimeError(s);
03632                 }
03633                 else if(info > 0)  
03634                 {
03635                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
03636                         char s[500];
03637                         sprintf(s, "The QR algorithm failed to compute all the eigenvalues!");
03638                         Utility::RunTimeError(s);
03639                 }
03640         }
03641         return w;
03642 }
03643 
03644 Vector<ComplexDouble> Lapack::Gees(Matrix<double>& A, Matrix<double>& T, Matrix<double>& Z)
03645 {
03646         if(A.Rows() != A.Columns())
03647         {
03648                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
03649                 Utility::RunTimeError("Matrix is not square!");
03650         }
03651         
03652         int info;
03653         T = A.Clone();
03654         int n = T.Rows();
03655         Vector<double> wr(n);
03656         Vector<double> wi(n);
03657 
03658         Z = Matrix<double>(n,n);
03659 
03660         int bs = 3*n;
03661         Vector<double> work(bs);
03662         int sdim;
03663 
03664         DGEES("V", "N", 0, &n, T.Data(), &n, &sdim, wr.Data(), wi.Data(), Z.Data(), &n, work.Data(), &bs, 0, &info);
03665 
03666         Vector<ComplexDouble> w(n);
03667         for(int i=0;i<n;i++)
03668         {
03669                 w[i] = ComplexDouble(wr[i],wi[i]);
03670         }
03671 
03672         if(info != 0)
03673         {
03674                 if(info < 0)
03675                 {
03676                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
03677                         char s[500];
03678                         sprintf(s, "%i'th argument of DGEES is an illegal value!", -info);
03679                         Utility::RunTimeError(s);
03680                 }
03681                 else if(info > 0)  
03682                 {
03683                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
03684                         char s[500];
03685                         sprintf(s, "The QR algorithm failed to compute all the eigenvalues!");
03686                         Utility::RunTimeError(s);
03687                 }
03688         }
03689         return w;
03690 }
03691 
03692 
03693 
03694 Vector<ComplexFloat> Lapack::Gees(Matrix<ComplexFloat>& A, Matrix<ComplexFloat>& T)
03695 {
03696         if(A.Rows() != A.Columns())
03697         {
03698                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
03699                 Utility::RunTimeError("Matrix is not square!");
03700         }
03701         
03702         int info;
03703         T = A.Clone();
03704         int n = T.Rows();
03705         Vector<ComplexFloat> w(n);
03706 
03707         int bs = 2*n;
03708         Vector<ComplexFloat> work(bs);
03709         Vector<float> rwork(n);
03710         int sdim;
03711 
03712         CGEES("N", "N", 0, &n, reinterpret_cast<MKL_Complex8*>(T.Data()), &n, &sdim, reinterpret_cast<MKL_Complex8*>(w.Data()), 0, &n, reinterpret_cast<MKL_Complex8*>(work.Data()), &bs, rwork.Data(), 0, &info);
03713 
03714         if(info != 0)
03715         {
03716                 if(info < 0)
03717                 {
03718                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
03719                         char s[500];
03720                         sprintf(s, "%i'th argument of CGEES is an illegal value!", -info);
03721                         Utility::RunTimeError(s);
03722                 }
03723                 else if(info > 0)  
03724                 {
03725                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
03726                         char s[500];
03727                         sprintf(s, "The QR algorithm failed to compute all the eigenvalues!");
03728                         Utility::RunTimeError(s);
03729                 }
03730         }
03731         return w;
03732 }
03733 
03734 
03735 Vector<ComplexFloat> Lapack::Gees(Matrix<ComplexFloat>& A, Matrix<ComplexFloat>& T, Matrix<ComplexFloat>& Z)
03736 {
03737         if(A.Rows() != A.Columns())
03738         {
03739                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
03740                 Utility::RunTimeError("Matrix is not square!");
03741         }
03742         
03743         int info;
03744         T = A.Clone();
03745         int n = T.Rows();
03746         Vector<ComplexFloat> w(n);
03747         Z = Matrix<ComplexFloat>(n,n);
03748 
03749         int bs = 2*n;
03750         Vector<ComplexFloat> work(bs);
03751         Vector<float> rwork(n);
03752         int sdim;
03753 
03754         CGEES("V", "N", 0, &n, reinterpret_cast<MKL_Complex8*>(T.Data()), &n, &sdim, reinterpret_cast<MKL_Complex8*>(w.Data()), reinterpret_cast<MKL_Complex8*>(Z.Data()), &n, reinterpret_cast<MKL_Complex8*>(work.Data()), &bs, rwork.Data(), 0, &info);
03755 
03756         if(info != 0)
03757         {
03758                 if(info < 0)
03759                 {
03760                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
03761                         char s[500];
03762                         sprintf(s, "%i'th argument of CGEES is an illegal value!", -info);
03763                         Utility::RunTimeError(s);
03764                 }
03765                 else if(info > 0)  
03766                 {
03767                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
03768                         char s[500];
03769                         sprintf(s, "The QR algorithm failed to compute all the eigenvalues!");
03770                         Utility::RunTimeError(s);
03771                 }
03772         }
03773         return w;
03774 }
03775 
03776 
03777 
03778 Vector<ComplexDouble> Lapack::Gees(Matrix<ComplexDouble>& A, Matrix<ComplexDouble>& T)
03779 {
03780         if(A.Rows() != A.Columns())
03781         {
03782                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
03783                 Utility::RunTimeError("Matrix is not square!");
03784         }
03785         
03786         int info;
03787         T = A.Clone();
03788         int n = T.Rows();
03789         Vector<ComplexDouble> w(n);
03790 
03791         int bs = 2*n;
03792         Vector<ComplexDouble> work(bs);
03793         Vector<double> rwork(n);
03794         int sdim;
03795 
03796         ZGEES("N", "N", 0, &n, reinterpret_cast<MKL_Complex16*>(T.Data()), &n, &sdim, reinterpret_cast<MKL_Complex16*>(w.Data()), 0, &n, reinterpret_cast<MKL_Complex16*>(work.Data()), &bs, rwork.Data(), 0, &info);
03797 
03798         if(info != 0)
03799         {
03800                 if(info < 0)
03801                 {
03802                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
03803                         char s[500];
03804                         sprintf(s, "%i'th argument of ZGEES is an illegal value!", -info);
03805                         Utility::RunTimeError(s);
03806                 }
03807                 else if(info > 0)  
03808                 {
03809                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
03810                         char s[500];
03811                         sprintf(s, "The QR algorithm failed to compute all the eigenvalues!");
03812                         Utility::RunTimeError(s);
03813                 }
03814         }
03815         return w;
03816 }
03817 
03818 
03819 Vector<ComplexDouble> Lapack::Gees(Matrix<ComplexDouble>& A, Matrix<ComplexDouble>& T, Matrix<ComplexDouble>& Z)
03820 {
03821         if(A.Rows() != A.Columns())
03822         {
03823                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
03824                 Utility::RunTimeError("Matrix is not square!");
03825         }
03826         
03827         int info;
03828         T = A.Clone();
03829         int n = T.Rows();
03830         Vector<ComplexDouble> w(n);
03831         Z = Matrix<ComplexDouble>(n,n);
03832 
03833         int bs = 2*n;
03834         Vector<ComplexDouble> work(bs);
03835         Vector<double> rwork(n);
03836         int sdim;
03837 
03838         ZGEES("V", "N", 0, &n, reinterpret_cast<MKL_Complex16*>(T.Data()), &n, &sdim, reinterpret_cast<MKL_Complex16*>(w.Data()), reinterpret_cast<MKL_Complex16*>(Z.Data()), &n, reinterpret_cast<MKL_Complex16*>(work.Data()), &bs, rwork.Data(), 0, &info);
03839 
03840         if(info != 0)
03841         {
03842                 if(info < 0)
03843                 {
03844                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
03845                         char s[500];
03846                         sprintf(s, "%i'th argument of ZGEES is an illegal value!", -info);
03847                         Utility::RunTimeError(s);
03848                 }
03849                 else if(info > 0)  
03850                 {
03851                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
03852                         char s[500];
03853                         sprintf(s, "The QR algorithm failed to compute all the eigenvalues!");
03854                         Utility::RunTimeError(s);
03855                 }
03856         }
03857         return w;
03858 }
03859 
03860 
03861 
03862 
03863 
03864 
03865 // Geev
03866 
03867 Vector<ComplexFloat> Lapack::Geev(Matrix<float>& A)
03868 {
03869         if(A.Rows() != A.Columns())
03870         {
03871                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
03872                 Utility::RunTimeError("Matrix is not square!");
03873         }
03874         
03875         int info;
03876         Matrix<float> F = A.Clone();
03877         int n = F.Rows();
03878         Vector<float> wr(n);
03879         Vector<float> wi(n);
03880 
03881         int bs = 4*n;
03882         Vector<float> work(bs);
03883 
03884         SGEEV("N", "N", &n, F.Data(), &n, wr.Data(), wi.Data(), 0, &n, 0, &n, work.Data(), &bs, &info);
03885 
03886         Vector<ComplexFloat> w(n);
03887         for(int i=0;i<n;i++)
03888         {
03889                 w[i] = ComplexFloat(wr[i],wi[i]);
03890         }
03891 
03892         if(info != 0)
03893         {
03894                 if(info < 0)
03895                 {
03896                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
03897                         char s[500];
03898                         sprintf(s, "%i'th argument of SGEEV is an illegal value!", -info);
03899                         Utility::RunTimeError(s);
03900                 }
03901                 else if(info > 0)  
03902                 {
03903                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
03904                         char s[500];
03905                         sprintf(s, "The QR algorithm failed to compute all the eigenvalues, and no eigenvectors have been computed!");
03906                         Utility::RunTimeError(s);
03907                 }
03908         }
03909         return w;
03910 }
03911 
03912 Vector<ComplexFloat> Lapack::Geev(Matrix<float>& A, Matrix<ComplexFloat>& Vl, Matrix<ComplexFloat>& Vr)
03913 {
03914         if(A.Rows() != A.Columns())
03915         {
03916                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
03917                 Utility::RunTimeError("Matrix is not square!");
03918         }
03919         
03920         int info;
03921         Matrix<float> F = A.Clone();
03922         int n = F.Rows();
03923         Vector<float> wr(n);
03924         Vector<float> wi(n);
03925 
03926         Vl = Matrix<ComplexFloat>(n,n);
03927         Vr = Matrix<ComplexFloat>(n,n);
03928         Matrix<float> vl(n,n);
03929         Matrix<float> vr(n,n);
03930 
03931         int bs = 4*n;
03932         Vector<float> work(bs);
03933 
03934         SGEEV("V", "V", &n, F.Data(), &n, wr.Data(), wi.Data(), vl.Data(), &n, vr.Data(), &n, work.Data(), &bs, &info);
03935 
03936         Vector<ComplexFloat> w(n);
03937         bool conj = false;
03938         for(int i=0;i<n;i++)
03939         {
03940                 w[i] = ComplexFloat(wr[i],wi[i]);
03941                 if(fabs(wi[i]) < numeric_limits<float>::epsilon())
03942                 {
03943                         for(int j=0;j<n;j++)
03944                         {
03945                                 Vl.ElemNC(j,i) = ComplexFloat(vl.ElemNC(j,i),0);
03946                                 Vr.ElemNC(j,i) = ComplexFloat(vr.ElemNC(j,i),0);
03947                         }
03948                 }
03949                 else
03950                 {
03951                         for(int j=0;j<n;j++)
03952                         {
03953                                 if(conj)
03954                                 {
03955                                         Vl.ElemNC(j,i) = ComplexFloat(vl.ElemNC(j,i-1),-vl.ElemNC(j,i));
03956                                         Vr.ElemNC(j,i) = ComplexFloat(vr.ElemNC(j,i-1),-vr.ElemNC(j,i));
03957                                 }
03958                                 else
03959                                 {
03960                                         Vl.ElemNC(j,i) = ComplexFloat(vl.ElemNC(j,i),vl.ElemNC(j,i+1));
03961                                         Vr.ElemNC(j,i) = ComplexFloat(vr.ElemNC(j,i),vr.ElemNC(j,i+1));
03962                                 }
03963                         }
03964                         conj = !conj;
03965                 }
03966         }
03967 
03968         if(info != 0)
03969         {
03970                 if(info < 0)
03971                 {
03972                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
03973                         char s[500];
03974                         sprintf(s, "%i'th argument of SGEEV is an illegal value!", -info);
03975                         Utility::RunTimeError(s);
03976                 }
03977                 else if(info > 0)  
03978                 {
03979                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
03980                         char s[500];
03981                         sprintf(s, "The QR algorithm failed to compute all the eigenvalues, and no eigenvectors have been computed!");
03982                         Utility::RunTimeError(s);
03983                 }
03984         }
03985         return w;
03986 }
03987 
03988 Vector<ComplexDouble> Lapack::Geev(Matrix<double>& A)
03989 {
03990         if(A.Rows() != A.Columns())
03991         {
03992                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
03993                 Utility::RunTimeError("Matrix is not square!");
03994         }
03995         
03996         int info;
03997         Matrix<double> F = A.Clone();
03998         int n = F.Rows();
03999         Vector<double> wr(n);
04000         Vector<double> wi(n);
04001 
04002         int bs = 4*n;
04003         Vector<double> work(bs);
04004 
04005         DGEEV("N", "N", &n, F.Data(), &n, wr.Data(), wi.Data(), 0, &n, 0, &n, work.Data(), &bs, &info);
04006 
04007         Vector<ComplexDouble> w(n);
04008         for(int i=0;i<n;i++)
04009         {
04010                 w[i] = ComplexDouble(wr[i],wi[i]);
04011         }
04012 
04013         if(info != 0)
04014         {
04015                 if(info < 0)
04016                 {
04017                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
04018                         char s[500];
04019                         sprintf(s, "%i'th argument of DGEEV is an illegal value!", -info);
04020                         Utility::RunTimeError(s);
04021                 }
04022                 else if(info > 0)  
04023                 {
04024                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
04025                         char s[500];
04026                         sprintf(s, "The QR algorithm failed to compute all the eigenvalues, and no eigenvectors have been computed!");
04027                         Utility::RunTimeError(s);
04028                 }
04029         }
04030         return w;
04031 }
04032 
04033 Vector<ComplexDouble> Lapack::Geev(Matrix<double>& A, Matrix<ComplexDouble>& Vl, Matrix<ComplexDouble>& Vr)
04034 {
04035         if(A.Rows() != A.Columns())
04036         {
04037                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
04038                 Utility::RunTimeError("Matrix is not square!");
04039         }
04040         
04041         int info;
04042         Matrix<double> F = A.Clone();
04043         int n = F.Rows();
04044         Vector<double> wr(n);
04045         Vector<double> wi(n);
04046 
04047         Vl = Matrix<ComplexDouble>(n,n);
04048         Vr = Matrix<ComplexDouble>(n,n);
04049         Matrix<double> vl(n,n);
04050         Matrix<double> vr(n,n);
04051 
04052         int bs = 4*n;
04053         Vector<double> work(bs);
04054 
04055         DGEEV("V", "V", &n, F.Data(), &n, wr.Data(), wi.Data(), vl.Data(), &n, vr.Data(), &n, work.Data(), &bs, &info);
04056 
04057         Vector<ComplexDouble> w(n);
04058         bool conj = false;
04059         for(int i=0;i<n;i++)
04060         {
04061                 w[i] = ComplexDouble(wr[i],wi[i]);
04062                 if(fabs(wi[i]) < numeric_limits<double>::epsilon())
04063                 {
04064                         for(int j=0;j<n;j++)
04065                         {
04066                                 Vl.ElemNC(j,i) = ComplexDouble(vl.ElemNC(j,i),0);
04067                                 Vr.ElemNC(j,i) = ComplexDouble(vr.ElemNC(j,i),0);
04068                         }
04069                 }
04070                 else
04071                 {
04072                         for(int j=0;j<n;j++)
04073                         {
04074                                 if(conj)
04075                                 {
04076                                         Vl.ElemNC(j,i) = ComplexDouble(vl.ElemNC(j,i-1),-vl.ElemNC(j,i));
04077                                         Vr.ElemNC(j,i) = ComplexDouble(vr.ElemNC(j,i-1),-vr.ElemNC(j,i));
04078                                 }
04079                                 else
04080                                 {
04081                                         Vl.ElemNC(j,i) = ComplexDouble(vl.ElemNC(j,i),vl.ElemNC(j,i+1));
04082                                         Vr.ElemNC(j,i) = ComplexDouble(vr.ElemNC(j,i),vr.ElemNC(j,i+1));
04083                                 }
04084                         }
04085                         conj = !conj;
04086                 }
04087         }
04088 
04089         if(info != 0)
04090         {
04091                 if(info < 0)
04092                 {
04093                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
04094                         char s[500];
04095                         sprintf(s, "%i'th argument of DGEEV is an illegal value!", -info);
04096                         Utility::RunTimeError(s);
04097                 }
04098                 else if(info > 0)  
04099                 {
04100                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
04101                         char s[500];
04102                         sprintf(s, "The QR algorithm failed to compute all the eigenvalues, and no eigenvectors have been computed!");
04103                         Utility::RunTimeError(s);
04104                 }
04105         }
04106         return w;
04107 }
04108 
04109 Vector<ComplexFloat> Lapack::Geev(Matrix<ComplexFloat>& A)
04110 {
04111         if(A.Rows() != A.Columns())
04112         {
04113                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
04114                 Utility::RunTimeError("Matrix is not square!");
04115         }
04116         
04117         int info;
04118         Matrix<ComplexFloat> F = A.Clone();
04119         int n = F.Rows();
04120         Vector<ComplexFloat> w(n);
04121 
04122         int bs = 4*n;
04123         Vector<ComplexFloat> work(bs);
04124         
04125         Vector<float> rwork(2*n);
04126 
04127         CGEEV("N", "N", &n, reinterpret_cast<MKL_Complex8*>(F.Data()), &n, reinterpret_cast<MKL_Complex8*>(w.Data()), 0, &n, 0, &n, reinterpret_cast<MKL_Complex8*>(work.Data()), &bs, rwork.Data(), &info);
04128 
04129         if(info != 0)
04130         {
04131                 if(info < 0)
04132                 {
04133                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
04134                         char s[500];
04135                         sprintf(s, "%i'th argument of CGEEV is an illegal value!", -info);
04136                         Utility::RunTimeError(s);
04137                 }
04138                 else if(info > 0)  
04139                 {
04140                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
04141                         char s[500];
04142                         sprintf(s, "The QR algorithm failed to compute all the eigenvalues, and no eigenvectors have been computed!");
04143                         Utility::RunTimeError(s);
04144                 }
04145         }
04146         return w;
04147 }
04148 
04149 Vector<ComplexFloat> Lapack::Geev(Matrix<ComplexFloat>& A, Matrix<ComplexFloat>& Vl, Matrix<ComplexFloat>& Vr)
04150 {
04151         if(A.Rows() != A.Columns())
04152         {
04153                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
04154                 Utility::RunTimeError("Matrix is not square!");
04155         }
04156         
04157         int info;
04158         Matrix<ComplexFloat> F = A.Clone();
04159         int n = F.Rows();
04160         Vector<ComplexFloat> w(n);
04161 
04162         Vl = Matrix<ComplexFloat>(n,n);
04163         Vr = Matrix<ComplexFloat>(n,n);
04164 
04165         int bs = 4*n;
04166         Vector<ComplexFloat> work(bs);
04167 
04168         Vector<float> rwork(2*n);
04169 
04170         CGEEV("V", "V", &n, reinterpret_cast<MKL_Complex8*>(F.Data()), &n, reinterpret_cast<MKL_Complex8*>(w.Data()), reinterpret_cast<MKL_Complex8*>(Vl.Data()), &n, reinterpret_cast<MKL_Complex8*>(Vr.Data()), &n, reinterpret_cast<MKL_Complex8*>(work.Data()), &bs, rwork.Data(), &info);
04171 
04172         if(info != 0)
04173         {
04174                 if(info < 0)
04175                 {
04176                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
04177                         char s[500];
04178                         sprintf(s, "%i'th argument of CGEEV is an illegal value!", -info);
04179                         Utility::RunTimeError(s);
04180                 }
04181                 else if(info > 0)  
04182                 {
04183                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
04184                         char s[500];
04185                         sprintf(s, "The QR algorithm failed to compute all the eigenvalues, and no eigenvectors have been computed!");
04186                         Utility::RunTimeError(s);
04187                 }
04188         }
04189         return w;
04190 }
04191 
04192 Vector<ComplexDouble> Lapack::Geev(Matrix<ComplexDouble>& A)
04193 {
04194         if(A.Rows() != A.Columns())
04195         {
04196                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
04197                 Utility::RunTimeError("Matrix is not square!");
04198         }
04199         
04200         int info;
04201         Matrix<ComplexDouble> F = A.Clone();
04202         int n = F.Rows();
04203         Vector<ComplexDouble> w(n);
04204 
04205         int bs = 4*n;
04206         Vector<ComplexDouble> work(bs);
04207         
04208         Vector<double> rwork(2*n);
04209 
04210         ZGEEV("N", "N", &n, reinterpret_cast<MKL_Complex16*>(F.Data()), &n, reinterpret_cast<MKL_Complex16*>(w.Data()), 0, &n, 0, &n, reinterpret_cast<MKL_Complex16*>(work.Data()), &bs, rwork.Data(), &info);
04211 
04212         if(info != 0)
04213         {
04214                 if(info < 0)
04215                 {
04216                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
04217                         char s[500];
04218                         sprintf(s, "%i'th argument of ZGEEV is an illegal value!", -info);
04219                         Utility::RunTimeError(s);
04220                 }
04221                 else if(info > 0)  
04222                 {
04223                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
04224                         char s[500];
04225                         sprintf(s, "The QR algorithm failed to compute all the eigenvalues, and no eigenvectors have been computed!");
04226                         Utility::RunTimeError(s);
04227                 }
04228         }
04229         return w;
04230 }
04231 
04232 Vector<ComplexDouble> Lapack::Geev(Matrix<ComplexDouble>& A, Matrix<ComplexDouble>& Vl, Matrix<ComplexDouble>& Vr)
04233 {
04234         if(A.Rows() != A.Columns())
04235         {
04236                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
04237                 Utility::RunTimeError("Matrix is not square!");
04238         }
04239         
04240         int info;
04241         Matrix<ComplexDouble> F = A.Clone();
04242         int n = F.Rows();
04243         Vector<ComplexDouble> w(n);
04244 
04245         Vl = Matrix<ComplexDouble>(n,n);
04246         Vr = Matrix<ComplexDouble>(n,n);
04247 
04248         int bs = 4*n;
04249         Vector<ComplexDouble> work(bs);
04250 
04251         Vector<double> rwork(2*n);
04252 
04253         ZGEEV("V", "V", &n, reinterpret_cast<MKL_Complex16*>(F.Data()), &n, reinterpret_cast<MKL_Complex16*>(w.Data()), reinterpret_cast<MKL_Complex16*>(Vl.Data()), &n, reinterpret_cast<MKL_Complex16*>(Vr.Data()), &n, reinterpret_cast<MKL_Complex16*>(work.Data()), &bs, rwork.Data(), &info);
04254 
04255         if(info != 0)
04256         {
04257                 if(info < 0)
04258                 {
04259                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
04260                         char s[500];
04261                         sprintf(s, "%i'th argument of ZGEEV is an illegal value!", -info);
04262                         Utility::RunTimeError(s);
04263                 }
04264                 else if(info > 0)  
04265                 {
04266                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
04267                         char s[500];
04268                         sprintf(s, "The QR algorithm failed to compute all the eigenvalues, and no eigenvectors have been computed!");
04269                         Utility::RunTimeError(s);
04270                 }
04271         }
04272         return w;
04273 }
04274 
04275 
04276 
04277 
04278 
04279 
04280 //=========================================================
04281 // Singular Value Decomposition (SVD)
04282 //=========================================================
04283 
04284 
04285 
04286 
04287 // Gesdd
04288 Vector<float> Lapack::Gesdd(Matrix<float>& A)
04289 {
04290         int info;
04291         Matrix<float> F = A.Clone();
04292         int m = F.Rows();
04293         int n = F.Columns();
04294         int min = (m>n)?n:m;
04295         int max = (m>n)?m:n;
04296 
04297         Vector<float> s(min);
04298 
04299         int temp = (max>6*min) ? max : 6*min;
04300         int bs = 3*min + temp;
04301         Vector<float> work(bs);
04302         Vector<int> iwork(8*min);
04303         SGESDD("N", &m, &n, F.Data(), &m, s.Data(), 0, &m, 0, &n, work.Data(), &bs, iwork.Data(), &info);
04304         if(info != 0)
04305         {
04306                 if(info < 0)
04307                 {
04308                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
04309                         char s[500];
04310                         sprintf(s, "%i'th argument of SGESDD is an illegal value!", -info);
04311                         Utility::RunTimeError(s);
04312                 }
04313                 else if(info > 0)  
04314                 {
04315                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
04316                         char s[500];
04317                         sprintf(s, "_bdsdc did not converge, updating process failed!");
04318                         Utility::RunTimeError(s);
04319                 }
04320         }
04321         return s;
04322 }
04323 
04324 Vector<float> Lapack::Gesdd(Matrix<float>& A, Matrix<float>& U, Matrix<float>& V)
04325 {
04326         int info;
04327         Matrix<float> F = A.Clone();
04328         int m = F.Rows();
04329         int n = F.Columns();
04330         int min = (m>n)?n:m;
04331         int max = (m>n)?m:n;
04332         U = Matrix<float>(m,m);
04333         V = Matrix<float>(n,n);
04334 
04335         Vector<float> s(min);
04336 
04337         int temp = (max > 4*min*min+4*min) ? max : 4*min*min+4*min;
04338         int bs = 3*min*min + temp;
04339         Vector<float> work(bs);
04340         Vector<int> iwork(8*min);
04341         SGESDD("A", &m, &n, F.Data(), &m, s.Data(), U.Data(), &m, V.Data(), &n, work.Data(), &bs, iwork.Data(), &info);
04342         if(info != 0)
04343         {
04344                 if(info < 0)
04345                 {
04346                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
04347                         char s[500];
04348                         sprintf(s, "%i'th argument of SGESDD is an illegal value!", -info);
04349                         Utility::RunTimeError(s);
04350                 }
04351                 else if(info > 0)  
04352                 {
04353                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
04354                         char s[500];
04355                         sprintf(s, "_bdsdc did not converge, updating process failed!");
04356                         Utility::RunTimeError(s);
04357                 }
04358         }
04359         
04360         V.Transpose();
04361                 
04362         return s;
04363 }
04364 
04365 Vector<double> Lapack::Gesdd(Matrix<double>& A)
04366 {
04367         int info;
04368         Matrix<double> F = A.Clone();
04369         int m = F.Rows();
04370         int n = F.Columns();
04371         int min = (m>n)?n:m;
04372         int max = (m>n)?m:n;
04373 
04374         Vector<double> s(min);
04375 
04376         int temp = (max>6*min) ? max : 6*min;
04377         int bs = 3*min + temp;
04378         Vector<double> work(bs);
04379         Vector<int> iwork(8*min);
04380         DGESDD("N", &m, &n, F.Data(), &m, s.Data(), 0, &m, 0, &n, work.Data(), &bs, iwork.Data(), &info);
04381         if(info != 0)
04382         {
04383                 if(info < 0)
04384                 {
04385                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
04386                         char s[500];
04387                         sprintf(s, "%i'th argument of DGESDD is an illegal value!", -info);
04388                         Utility::RunTimeError(s);
04389                 }
04390                 else if(info > 0)  
04391                 {
04392                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
04393                         char s[500];
04394                         sprintf(s, "_bdsdc did not converge, updating process failed!");
04395                         Utility::RunTimeError(s);
04396                 }
04397         }
04398         return s;
04399 }
04400 
04401 Vector<double> Lapack::Gesdd(Matrix<double>& A, Matrix<double>& U, Matrix<double>& V)
04402 {
04403         int info;
04404         Matrix<double> F = A.Clone();
04405         int m = F.Rows();
04406         int n = F.Columns();
04407         int min = (m>n)?n:m;
04408         int max = (m>n)?m:n;
04409         U = Matrix<double>(m,m);
04410         V = Matrix<double>(n,n);
04411 
04412         Vector<double> s(min);
04413 
04414         int temp = (max > 4*min*min+4*min) ? max : 4*min*min+4*min;
04415         int bs = 3*min*min + temp;
04416         Vector<double> work(bs);
04417         Vector<int> iwork(8*min);
04418         DGESDD("A", &m, &n, F.Data(), &m, s.Data(), U.Data(), &m, V.Data(), &n, work.Data(), &bs, iwork.Data(), &info);
04419         if(info != 0)
04420         {
04421                 if(info < 0)
04422                 {
04423                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
04424                         char s[500];
04425                         sprintf(s, "%i'th argument of DGESDD is an illegal value!", -info);
04426                         Utility::RunTimeError(s);
04427                 }
04428                 else if(info > 0)  
04429                 {
04430                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
04431                         char s[500];
04432                         sprintf(s, "_bdsdc did not converge, updating process failed!");
04433                         Utility::RunTimeError(s);
04434                 }
04435         }
04436         
04437         V.Transpose();
04438                 
04439         return s;
04440 }
04441 
04442 Vector<float> Lapack::Gesdd(Matrix<ComplexFloat>& A)
04443 {
04444         int info;
04445         Matrix<ComplexFloat> F = A.Clone();
04446         int m = F.Rows();
04447         int n = F.Columns();
04448         int min = (m>n)?n:m;
04449         int max = (m>n)?m:n;
04450 
04451         Vector<float> s(min);
04452 
04453         int bs = 2*min + max;
04454         Vector<ComplexFloat> work(bs);
04455         Vector<int> iwork(8*min);
04456         Vector<float> rwork(5*min);
04457 
04458         CGESDD("N", &m, &n, reinterpret_cast<MKL_Complex8*>(F.Data()), &m, s.Data(), 0, &m, 0, &n, reinterpret_cast<MKL_Complex8*>(work.Data()), &bs, rwork.Data(), iwork.Data(), &info);
04459         if(info != 0)
04460         {
04461                 if(info < 0)
04462                 {
04463                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
04464                         char s[500];
04465                         sprintf(s, "%i'th argument of CGESDD is an illegal value!", -info);
04466                         Utility::RunTimeError(s);
04467                 }
04468                 else if(info > 0)  
04469                 {
04470                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
04471                         char s[500];
04472                         sprintf(s, "_bdsdc did not converge, updating process failed!");
04473                         Utility::RunTimeError(s);
04474                 }
04475         }
04476         return s;
04477 }
04478 
04479 Vector<float> Lapack::Gesdd(Matrix<ComplexFloat>& A, Matrix<ComplexFloat>& U, Matrix<ComplexFloat>& V)
04480 {
04481         int info;
04482         Matrix<ComplexFloat> F = A.Clone();
04483         int m = F.Rows();
04484         int n = F.Columns();
04485         int min = (m>n)?n:m;
04486         int max = (m>n)?m:n;
04487         U = Matrix<ComplexFloat>(m,m);
04488         V = Matrix<ComplexFloat>(n,n);
04489 
04490         Vector<float> s(min);
04491 
04492         int bs = min*min + max + 2*min;
04493         Vector<ComplexFloat> work(bs);
04494         Vector<int> iwork(8*min);
04495         Vector<float> rwork(5*min*min + 7*min);
04496         CGESDD("A", &m, &n, reinterpret_cast<MKL_Complex8*>(F.Data()), &m, s.Data(), reinterpret_cast<MKL_Complex8*>(U.Data()), &m, reinterpret_cast<MKL_Complex8*>(V.Data()), &n, reinterpret_cast<MKL_Complex8*>(work.Data()), &bs, rwork.Data(), iwork.Data(), &info);
04497         if(info != 0)
04498         {
04499                 if(info < 0)
04500                 {
04501                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
04502                         char s[500];
04503                         sprintf(s, "%i'th argument of CGESDD is an illegal value!", -info);
04504                         Utility::RunTimeError(s);
04505                 }
04506                 else if(info > 0)  
04507                 {
04508                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
04509                         char s[500];
04510                         sprintf(s, "_bdsdc did not converge, updating process failed!");
04511                         Utility::RunTimeError(s);
04512                 }
04513         }
04514         
04515         V.Transpose();
04516         for(int i=0; i<V.Length(); i++)
04517         {
04518                 V.ElemNC(i) = conj(V.ElemNC(i));
04519         }
04520                 
04521         return s;
04522 }
04523 
04524 Vector<double> Lapack::Gesdd(Matrix<ComplexDouble>& A)
04525 {
04526         int info;
04527         Matrix<ComplexDouble> F = A.Clone();
04528         int m = F.Rows();
04529         int n = F.Columns();
04530         int min = (m>n)?n:m;
04531         int max = (m>n)?m:n;
04532 
04533         Vector<double> s(min);
04534 
04535         int bs = 2*min + max;
04536         Vector<ComplexDouble> work(bs);
04537         Vector<int> iwork(8*min);
04538         Vector<double> rwork(5*min);
04539 
04540         ZGESDD("N", &m, &n, reinterpret_cast<MKL_Complex16*>(F.Data()), &m, s.Data(), 0, &m, 0, &n, reinterpret_cast<MKL_Complex16*>(work.Data()), &bs, rwork.Data(), iwork.Data(), &info);
04541         if(info != 0)
04542         {
04543                 if(info < 0)
04544                 {
04545                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
04546                         char s[500];
04547                         sprintf(s, "%i'th argument of ZGESDD is an illegal value!", -info);
04548                         Utility::RunTimeError(s);
04549                 }
04550                 else if(info > 0)  
04551                 {
04552                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
04553                         char s[500];
04554                         sprintf(s, "_bdsdc did not converge, updating process failed!");
04555                         Utility::RunTimeError(s);
04556                 }
04557         }
04558         return s;
04559 }
04560 
04561 Vector<double> Lapack::Gesdd(Matrix<ComplexDouble>& A, Matrix<ComplexDouble>& U, Matrix<ComplexDouble>& V)
04562 {
04563         int info;
04564         Matrix<ComplexDouble> F = A.Clone();
04565         int m = F.Rows();
04566         int n = F.Columns();
04567         int min = (m>n)?n:m;
04568         int max = (m>n)?m:n;
04569         U = Matrix<ComplexDouble>(m,m);
04570         V = Matrix<ComplexDouble>(n,n);
04571 
04572         Vector<double> s(min);
04573 
04574         int bs = min*min + max + 2*min;
04575         Vector<ComplexDouble> work(bs);
04576         Vector<int> iwork(8*min);
04577         Vector<double> rwork(5*min*min + 7*min);
04578         ZGESDD("A", &m, &n, reinterpret_cast<MKL_Complex16*>(F.Data()), &m, s.Data(), reinterpret_cast<MKL_Complex16*>(U.Data()), &m, reinterpret_cast<MKL_Complex16*>(V.Data()), &n, reinterpret_cast<MKL_Complex16*>(work.Data()), &bs, rwork.Data(), iwork.Data(), &info);
04579         if(info != 0)
04580         {
04581                 if(info < 0)
04582                 {
04583                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
04584                         char s[500];
04585                         sprintf(s, "%i'th argument of ZGESDD is an illegal value!", -info);
04586                         Utility::RunTimeError(s);
04587                 }
04588                 else if(info > 0)  
04589                 {
04590                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
04591                         char s[500];
04592                         sprintf(s, "_bdsdc did not converge, updating process failed!");
04593                         Utility::RunTimeError(s);
04594                 }
04595         }
04596         
04597         V.Transpose();
04598         for(int i=0; i<V.Length(); i++)
04599         {
04600                 V.ElemNC(i) = conj(V.ElemNC(i));
04601         }
04602 
04603         return s;
04604 }
04605 
04606 
04607 
04608 
04609 
04610 
04611 
04612 
04613 
04614 
04615 
04616 
04617 //=========================================================
04618 // Generalized Symmetric Definite eigenproblem
04619 //=========================================================
04620 
04621 
04622         // Sygv
04623         
04624         
04625         
04626         
04627         // Sygvx
04628         
04629         
04630         
04631         
04632         
04633         // Sygvd,Hegvd
04634 Vector<float> Lapack::Sygvd(Matrix<float>& A, Matrix<float>& B, int ptype)
04635 {
04636         if(ptype < 1 || ptype > 3 )
04637         {
04638                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
04639                 Utility::RunTimeError("Incorrext problem type. ptype can be 1, 2, or 3!");
04640         }
04641         
04642         if(A.Rows() != A.Columns() || B.Rows() != B.Columns())
04643         {
04644                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
04645                 Utility::RunTimeError("Matrix is not square!");
04646         }
04647         if(A.Rows() != B.Rows())
04648         {
04649                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
04650                 Utility::RunTimeError("Matrix dimensions does not match!");
04651         }
04652         
04653         int info;
04654         Matrix<float> A2 = A.Clone();
04655         Matrix<float> B2 = B.Clone();
04656         int n = A2.Rows();
04657         Vector<float> w(n);
04658 
04659         int bs = 2*n+1;
04660         Vector<float> work(bs);
04661 
04662         int iwork_size = 1;
04663         Vector<int> iwork(iwork_size);
04664 
04665         SSYGVD(&ptype, "N", "U", &n, A2.Data(), &n, B2.Data(), &n, w.Data(), work.Data(), &bs, iwork.Data(), &iwork_size, &info);
04666         if(info != 0)
04667         {
04668                 if(info < 0)
04669                 {
04670                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
04671                         char s[500];
04672                         sprintf(s, "%i'th argument of SSYGVD is an illegal value!", -info);
04673                         Utility::RunTimeError(s);
04674                 }
04675                 else if(info > 0)  
04676                 {
04677                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
04678                         cerr << "spotrf/dpotrf and ssyev/dsyev returned an error code: ";
04679                         if(info <= n)
04680                         {
04681                                 cerr << "ssyev/dsyev failed to converge, and " << info <<  " off-diagonal elements of an intermediate tridiagonal did not converge to zero!";
04682                         }
04683                         else if(info > n && info <=  2*n)
04684                         {
04685                                 cerr << "The leading minor of order " << info-n << " of B is not positive-definite. The factorization of B could not be completed and no eigenvalues or eigenvectors were computed!";
04686                         }
04687                         Utility::RunTimeError("");
04688                 }
04689         }
04690         return w;
04691 }
04692 
04693 Vector<float> Lapack::Sygvd(Matrix<float>& A, Matrix<float>& B, Matrix<float>& EV, int ptype)
04694 {
04695         if(ptype < 1 || ptype > 3 )
04696         {
04697                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
04698                 Utility::RunTimeError("Incorrext problem type. ptype can be 1, 2, or 3!");
04699         }
04700         
04701         if(A.Rows() != A.Columns() || B.Rows() != B.Columns())
04702         {
04703                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
04704                 Utility::RunTimeError("Matrix is not square!");
04705         }
04706         if(A.Rows() != B.Rows())
04707         {
04708                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
04709                 Utility::RunTimeError("Matrix dimensions does not match!");
04710         }
04711         
04712         int info;
04713         EV = A.Clone();
04714         Matrix<float> B2 = B.Clone();
04715         int n = B2.Rows();
04716         Vector<float> w(n);
04717 
04718         int bs = 2*n*n + 6*n +1;
04719         Vector<float> work(bs);
04720 
04721         int iwork_size = 5*n + 3;
04722         Vector<int> iwork(iwork_size);
04723 
04724         SSYGVD(&ptype, "V", "U", &n, EV.Data(), &n, B2.Data(), &n, w.Data(), work.Data(), &bs, iwork.Data(), &iwork_size, &info);
04725         if(info != 0)
04726         {
04727                 if(info < 0)
04728                 {
04729                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
04730                         char s[500];
04731                         sprintf(s, "%i'th argument of SSYGVD is an illegal value!", -info);
04732                         Utility::RunTimeError(s);
04733                 }
04734                 else if(info > 0)  
04735                 {
04736                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
04737                         cerr << "spotrf/dpotrf and ssyev/dsyev returned an error code: ";
04738                         if(info <= n)
04739                         {
04740                                 cerr << "ssyev/dsyev failed to converge, and " << info <<  " off-diagonal elements of an intermediate tridiagonal did not converge to zero!";
04741                         }
04742                         else if(info > n && info <=  2*n)
04743                         {
04744                                 cerr << "The leading minor of order " << info-n << " of B is not positive-definite. The factorization of B could not be completed and no eigenvalues or eigenvectors were computed!";
04745                         }
04746                         Utility::RunTimeError("");
04747                 }
04748         }
04749         return w;
04750 }
04751 
04752 Vector<double> Lapack::Sygvd(Matrix<double>& A, Matrix<double>& B, int ptype)
04753 {
04754         if(ptype < 1 || ptype > 3 )
04755         {
04756                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
04757                 Utility::RunTimeError("Incorrext problem type. ptype can be 1, 2, or 3!");
04758         }
04759         
04760         if(A.Rows() != A.Columns() || B.Rows() != B.Columns())
04761         {
04762                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
04763                 Utility::RunTimeError("Matrix is not square!");
04764         }
04765         if(A.Rows() != B.Rows())
04766         {
04767                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
04768                 Utility::RunTimeError("Matrix dimensions does not match!");
04769         }
04770         
04771         int info;
04772         Matrix<double> A2 = A.Clone();
04773         Matrix<double> B2 = B.Clone();
04774         int n = A2.Rows();
04775         Vector<double> w(n);
04776 
04777         int bs = 2*n+1;
04778         Vector<double> work(bs);
04779 
04780         int iwork_size = 1;
04781         Vector<int> iwork(iwork_size);
04782 
04783         DSYGVD(&ptype, "N", "U", &n, A2.Data(), &n, B2.Data(), &n, w.Data(), work.Data(), &bs, iwork.Data(), &iwork_size, &info);
04784         if(info != 0)
04785         {
04786                 if(info < 0)
04787                 {
04788                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
04789                         char s[500];
04790                         sprintf(s, "%i'th argument of DSYGVD is an illegal value!", -info);
04791                         Utility::RunTimeError(s);
04792                 }
04793                 else if(info > 0)  
04794                 {
04795                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
04796                         cerr << "spotrf/dpotrf and ssyev/dsyev returned an error code: ";
04797                         if(info <= n)
04798                         {
04799                                 cerr << "ssyev/dsyev failed to converge, and " << info <<  " off-diagonal elements of an intermediate tridiagonal did not converge to zero!";
04800                         }
04801                         else if(info > n && info <=  2*n)
04802                         {
04803                                 cerr << "The leading minor of order " << info-n << " of B is not positive-definite. The factorization of B could not be completed and no eigenvalues or eigenvectors were computed!";
04804                         }
04805                         Utility::RunTimeError("");
04806                 }
04807         }
04808         return w;
04809 }
04810 
04811 Vector<double> Lapack::Sygvd(Matrix<double>& A, Matrix<double>& B, Matrix<double>& EV, int ptype)
04812 {
04813         if(ptype < 1 || ptype > 3 )
04814         {
04815                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
04816                 Utility::RunTimeError("Incorrext problem type. ptype can be 1, 2, or 3!");
04817         }
04818         
04819         if(A.Rows() != A.Columns() || B.Rows() != B.Columns())
04820         {
04821                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
04822                 Utility::RunTimeError("Matrix is not square!");
04823         }
04824         if(A.Rows() != B.Rows())
04825         {
04826                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
04827                 Utility::RunTimeError("Matrix dimensions does not match!");
04828         }
04829         
04830         int info;
04831         EV = A.Clone();
04832         Matrix<double> B2 = B.Clone();
04833         int n = B2.Rows();
04834         Vector<double> w(n);
04835 
04836         int bs = 2*n*n + 6*n +1;
04837         Vector<double> work(bs);
04838 
04839         int iwork_size = 5*n + 3;
04840         Vector<int> iwork(iwork_size);
04841 
04842         DSYGVD(&ptype, "V", "U", &n, EV.Data(), &n, B2.Data(), &n, w.Data(), work.Data(), &bs, iwork.Data(), &iwork_size, &info);
04843         if(info != 0)
04844         {
04845                 if(info < 0)
04846                 {
04847                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
04848                         char s[500];
04849                         sprintf(s, "%i'th argument of DSYGVD is an illegal value!", -info);
04850                         Utility::RunTimeError(s);
04851                 }
04852                 else if(info > 0)  
04853                 {
04854                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
04855                         cerr << "spotrf/dpotrf and ssyev/dsyev returned an error code: ";
04856                         if(info <= n)
04857                         {
04858                                 cerr << "ssyev/dsyev failed to converge, and " << info <<  " off-diagonal elements of an intermediate tridiagonal did not converge to zero!";
04859                         }
04860                         else if(info > n && info <=  2*n)
04861                         {
04862                                 cerr << "The leading minor of order " << info-n << " of B is not positive-definite. The factorization of B could not be completed and no eigenvalues or eigenvectors were computed!";
04863                         }
04864                         Utility::RunTimeError("");
04865                 }
04866         }
04867         return w;
04868 }
04869 
04870 
04871 // Hegvd
04872 Vector<float> Lapack::Hegvd(Matrix<ComplexFloat>& A, Matrix<ComplexFloat>& B, int ptype)
04873 {
04874         if(ptype < 1 || ptype > 3 )
04875         {
04876                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
04877                 Utility::RunTimeError("Incorrext problem type. ptype can be 1, 2, or 3!");
04878         }
04879         
04880         if(A.Rows() != A.Columns() || B.Rows() != B.Columns())
04881         {
04882                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
04883                 Utility::RunTimeError("Matrix is not square!");
04884         }
04885         if(A.Rows() != B.Rows())
04886         {
04887                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
04888                 Utility::RunTimeError("Matrix dimensions does not match!");
04889         }
04890         
04891         int info;
04892         Matrix<ComplexFloat> A2 = A.Clone();
04893         Matrix<ComplexFloat> B2 = B.Clone();
04894         int n = A2.Rows();
04895         Vector<float> w(n);
04896 
04897         int bs = n+1;
04898         Vector<ComplexFloat> work(bs);
04899 
04900         int rwork_size = n;
04901         Vector<float> rwork(rwork_size);
04902 
04903         int iwork_size = 1;
04904         Vector<int> iwork(iwork_size);
04905 
04906         CHEGVD(&ptype, "N", "U", &n, reinterpret_cast<MKL_Complex8*>(A2.Data()), &n, reinterpret_cast<MKL_Complex8*>(B2.Data()), &n, w.Data(), reinterpret_cast<MKL_Complex8*>(work.Data()), &bs, rwork.Data(), &rwork_size, iwork.Data(), &iwork_size, &info);
04907         if(info != 0)
04908         {
04909                 if(info < 0)
04910                 {
04911                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
04912                         char s[500];
04913                         sprintf(s, "%i'th argument of CHEGVD is an illegal value!", -info);
04914                         Utility::RunTimeError(s);
04915                 }
04916                 else if(info > 0)  
04917                 {
04918                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
04919                         cerr << "cpotrf/zpotrf and cheev/zheev returned an error code: ";
04920                         if(info <= n)
04921                         {
04922                                 cerr << "cheev/zheev failed to converge, and " << info <<  " off-diagonal elements of an intermediate tridiagonal did not converge to zero!";
04923                         }
04924                         else if(info > n && info <=  2*n)
04925                         {
04926                                 cerr << "The leading minor of order " << info-n << " of B is not positive-definite. The factorization of B could not be completed and no eigenvalues or eigenvectors were computed!";
04927                         }
04928                         Utility::RunTimeError("");
04929                 }
04930         }
04931         return w;
04932 }
04933 
04934 Vector<float> Lapack::Hegvd(Matrix<ComplexFloat>& A, Matrix<ComplexFloat>& B, Matrix<ComplexFloat>& EV, int ptype)
04935 {
04936         if(ptype < 1 || ptype > 3 )
04937         {
04938                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
04939                 Utility::RunTimeError("Incorrext problem type. ptype can be 1, 2, or 3!");
04940         }
04941         
04942         if(A.Rows() != A.Columns() || B.Rows() != B.Columns())
04943         {
04944                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
04945                 Utility::RunTimeError("Matrix is not square!");
04946         }
04947         if(A.Rows() != B.Rows())
04948         {
04949                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
04950                 Utility::RunTimeError("Matrix dimensions does not match!");
04951         }
04952         
04953         int info;
04954         EV = A.Clone();
04955         Matrix<ComplexFloat> B2 = B.Clone();
04956         int n = B2.Rows();
04957         Vector<float> w(n);
04958 
04959         int bs = n*n+2*n;
04960         Vector<ComplexFloat> work(bs);
04961 
04962         int rwork_size = 2*n*n+5*n+1;
04963         Vector<float> rwork(rwork_size);
04964 
04965         int iwork_size = 5*n+3;
04966         Vector<int> iwork(iwork_size);
04967 
04968         CHEGVD(&ptype, "V", "U", &n, reinterpret_cast<MKL_Complex8*>(EV.Data()), &n, reinterpret_cast<MKL_Complex8*>(B2.Data()), &n, w.Data(), reinterpret_cast<MKL_Complex8*>(work.Data()), &bs, rwork.Data(), &rwork_size, iwork.Data(), &iwork_size, &info);
04969         if(info != 0)
04970         {
04971                 if(info < 0)
04972                 {
04973                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
04974                         char s[500];
04975                         sprintf(s, "%i'th argument of CHEGVD is an illegal value!", -info);
04976                         Utility::RunTimeError(s);
04977                 }
04978                 else if(info > 0)  
04979                 {
04980                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
04981                         cerr << "cpotrf/zpotrf and cheev/zheev returned an error code: ";
04982                         if(info <= n)
04983                         {
04984                                 cerr << "cheev/zheev failed to converge, and " << info <<  " off-diagonal elements of an intermediate tridiagonal did not converge to zero!";
04985                         }
04986                         else if(info > n && info <=  2*n)
04987                         {
04988                                 cerr << "The leading minor of order " << info-n << " of B is not positive-definite. The factorization of B could not be completed and no eigenvalues or eigenvectors were computed!";
04989                         }
04990                         Utility::RunTimeError("");
04991                 }
04992         }
04993         return w;
04994 }
04995 
04996 
04997 Vector<double> Lapack::Hegvd(Matrix<ComplexDouble>& A, Matrix<ComplexDouble>& B, int ptype)
04998 {
04999         if(ptype < 1 || ptype > 3 )
05000         {
05001                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
05002                 Utility::RunTimeError("Incorrext problem type. ptype can be 1, 2, or 3!");
05003         }
05004         
05005         if(A.Rows() != A.Columns() || B.Rows() != B.Columns())
05006         {
05007                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
05008                 Utility::RunTimeError("Matrix is not square!");
05009         }
05010         if(A.Rows() != B.Rows())
05011         {
05012                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
05013                 Utility::RunTimeError("Matrix dimensions does not match!");
05014         }
05015         
05016         int info;
05017         Matrix<ComplexDouble> A2 = A.Clone();
05018         Matrix<ComplexDouble> B2 = B.Clone();
05019         int n = A2.Rows();
05020         Vector<double> w(n);
05021 
05022         int bs = n+1;
05023         Vector<ComplexDouble> work(bs);
05024 
05025         int rwork_size = n;
05026         Vector<double> rwork(rwork_size);
05027 
05028         int iwork_size = 1;
05029         Vector<int> iwork(iwork_size);
05030 
05031         ZHEGVD(&ptype, "N", "U", &n, reinterpret_cast<MKL_Complex16*>(A2.Data()), &n, reinterpret_cast<MKL_Complex16*>(B2.Data()), &n, w.Data(), reinterpret_cast<MKL_Complex16*>(work.Data()), &bs, rwork.Data(), &rwork_size, iwork.Data(), &iwork_size, &info);
05032         if(info != 0)
05033         {
05034                 if(info < 0)
05035                 {
05036                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
05037                         char s[500];
05038                         sprintf(s, "%i'th argument of ZHEGVD is an illegal value!", -info);
05039                         Utility::RunTimeError(s);
05040                 }
05041                 else if(info > 0)  
05042                 {
05043                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
05044                         cerr << "cpotrf/zpotrf and cheev/zheev returned an error code: ";
05045                         if(info <= n)
05046                         {
05047                                 cerr << "cheev/zheev failed to converge, and " << info <<  " off-diagonal elements of an intermediate tridiagonal did not converge to zero!";
05048                         }
05049                         else if(info > n && info <=  2*n)
05050                         {
05051                                 cerr << "The leading minor of order " << info-n << " of B is not positive-definite. The factorization of B could not be completed and no eigenvalues or eigenvectors were computed!";
05052                         }
05053                         Utility::RunTimeError("");
05054                 }
05055         }
05056         return w;
05057 }
05058 
05059 Vector<double> Lapack::Hegvd(Matrix<ComplexDouble>& A, Matrix<ComplexDouble>& B, Matrix<ComplexDouble>& EV, int ptype)
05060 {
05061         if(ptype < 1 || ptype > 3 )
05062         {
05063                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
05064                 Utility::RunTimeError("Incorrext problem type. ptype can be 1, 2, or 3!");
05065         }
05066         
05067         if(A.Rows() != A.Columns() || B.Rows() != B.Columns())
05068         {
05069                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
05070                 Utility::RunTimeError("Matrix is not square!");
05071         }
05072         if(A.Rows() != B.Rows())
05073         {
05074                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
05075                 Utility::RunTimeError("Matrix dimensions does not match!");
05076         }
05077         
05078         int info;
05079         EV = A.Clone();
05080         Matrix<ComplexDouble> B2 = B.Clone();
05081         int n = B2.Rows();
05082         Vector<double> w(n);
05083 
05084         int bs = n*n+2*n;
05085         Vector<ComplexDouble> work(bs);
05086 
05087         int rwork_size = 2*n*n+5*n+1;
05088         Vector<double> rwork(rwork_size);
05089 
05090         int iwork_size = 5*n+3;
05091         Vector<int> iwork(iwork_size);
05092 
05093         ZHEGVD(&ptype, "V", "U", &n, reinterpret_cast<MKL_Complex16*>(EV.Data()), &n, reinterpret_cast<MKL_Complex16*>(B2.Data()), &n, w.Data(), reinterpret_cast<MKL_Complex16*>(work.Data()), &bs, rwork.Data(), &rwork_size, iwork.Data(), &iwork_size, &info);
05094         if(info != 0)
05095         {
05096                 if(info < 0)
05097                 {
05098                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
05099                         char s[500];
05100                         sprintf(s, "%i'th argument of ZHEGVD is an illegal value!", -info);
05101                         Utility::RunTimeError(s);
05102                 }
05103                 else if(info > 0)  
05104                 {
05105                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
05106                         cerr << "cpotrf/zpotrf and cheev/zheev returned an error code: ";
05107                         if(info <= n)
05108                         {
05109                                 cerr << "cheev/zheev failed to converge, and " << info <<  " off-diagonal elements of an intermediate tridiagonal did not converge to zero!";
05110                         }
05111                         else if(info > n && info <=  2*n)
05112                         {
05113                                 cerr << "The leading minor of order " << info-n << " of B is not positive-definite. The factorization of B could not be completed and no eigenvalues or eigenvectors were computed!";
05114                         }
05115                         Utility::RunTimeError("");
05116                 }
05117         }
05118         return w;
05119 }
05120 
05121 
05122 
05123 
05124 
05125 //=========================================================
05126 // Generalized Non-Symmetric eigenvalue problem
05127 //=========================================================
05128 
05129 
05130 // Gges
05131 void Lapack::Gges(Matrix<float>& A, Matrix<float>& B, Vector<ComplexFloat>& alpha, Vector<float>& beta, Matrix<float>& S, Matrix<float>& T)
05132 {
05133         if(A.Rows() != A.Columns() || B.Rows() != B.Columns())
05134         {
05135                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
05136                 Utility::RunTimeError("Matrix is not square!");
05137         }
05138         if(A.Rows() != B.Rows())
05139         {
05140                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
05141                 Utility::RunTimeError("Matrix dimensions does not match!");
05142         }
05143 
05144         
05145         int info;
05146         S = A.Clone();
05147         T = B.Clone();
05148         int n = A.Rows();
05149         Vector<float> wr(n);
05150         Vector<float> wi(n);
05151         beta = Vector<float>(n);
05152 
05153         int bs = 39*n+23;
05154         Vector<float> work(bs);
05155 
05156         int sdim;
05157 
05158         SGGES("N", "N", "N", 0, &n, S.Data(), &n, T.Data(), &n, &sdim, wr.Data(), wi.Data(), beta.Data(), 0, &n, 0, &n, work.Data(), &bs, 0, &info);
05159 
05160         alpha = Vector<ComplexFloat>(n);
05161         for(int i=0;i<n;i++)
05162         {
05163                 alpha[i] = ComplexFloat(wr[i],wi[i]);
05164         }
05165 
05166         if(info != 0)
05167         {
05168                 if(info < 0)
05169                 {
05170                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
05171                         char s[500];
05172                         sprintf(s, "%i'th argument of SGGES is an illegal value!", -info);
05173                         Utility::RunTimeError(s);
05174                 }
05175                 else if(info > 0)  
05176                 {
05177                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
05178                         char s[500];
05179                         sprintf(s, "Lapack error!");
05180                         Utility::RunTimeError(s);
05181                 }
05182         }
05183 }
05184 
05185 void Lapack::Gges(Matrix<float>& A, Matrix<float>& B, Vector<ComplexFloat>& alpha, Vector<float>& beta, Matrix<float>& S, Matrix<float>& T, Matrix<float>& Vl, Matrix<float>& Vr)
05186 {
05187         if(A.Rows() != A.Columns() || B.Rows() != B.Columns())
05188         {
05189                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
05190                 Utility::RunTimeError("Matrix is not square!");
05191         }
05192         if(A.Rows() != B.Rows())
05193         {
05194                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
05195                 Utility::RunTimeError("Matrix dimensions does not match!");
05196         }
05197 
05198         
05199         int info;
05200         S = A.Clone();
05201         T = B.Clone();
05202         int n = A.Rows();
05203         Vector<float> wr(n);
05204         Vector<float> wi(n);
05205         beta = Vector<float>(n);
05206 
05207         Vl = Matrix<float>(n,n);
05208         Vr = Matrix<float>(n,n);
05209 
05210         int bs = 39*n+23;
05211         Vector<float> work(bs);
05212 
05213         int sdim;
05214 
05215         SGGES("V", "V", "N", 0, &n, S.Data(), &n, T.Data(), &n, &sdim, wr.Data(), wi.Data(), beta.Data(), Vl.Data(), &n, Vr.Data(), &n, work.Data(), &bs, 0, &info);
05216 
05217         alpha = Vector<ComplexFloat>(n);
05218         for(int i=0;i<n;i++)
05219         {
05220                 alpha[i] = ComplexFloat(wr[i],wi[i]);
05221         }
05222 
05223         if(info != 0)
05224         {
05225                 if(info < 0)
05226                 {
05227                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
05228                         char s[500];
05229                         sprintf(s, "%i'th argument of SGGES is an illegal value!", -info);
05230                         Utility::RunTimeError(s);
05231                 }
05232                 else if(info > 0)  
05233                 {
05234                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
05235                         char s[500];
05236                         sprintf(s, "Lapack error!");
05237                         Utility::RunTimeError(s);
05238                 }
05239         }
05240 }
05241 
05242 void Lapack::Gges(Matrix<double>& A, Matrix<double>& B, Vector<ComplexDouble>& alpha, Vector<double>& beta, Matrix<double>& S, Matrix<double>& T)
05243 {
05244         if(A.Rows() != A.Columns() || B.Rows() != B.Columns())
05245         {
05246                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
05247                 Utility::RunTimeError("Matrix is not square!");
05248         }
05249         if(A.Rows() != B.Rows())
05250         {
05251                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
05252                 Utility::RunTimeError("Matrix dimensions does not match!");
05253         }
05254 
05255         
05256         int info;
05257         S = A.Clone();
05258         T = B.Clone();
05259         int n = A.Rows();
05260         Vector<double> wr(n);
05261         Vector<double> wi(n);
05262         beta = Vector<double>(n);
05263 
05264         int bs = 39*n+23;
05265         Vector<double> work(bs);
05266 
05267         int sdim;
05268 
05269         DGGES("N", "N", "N", 0, &n, S.Data(), &n, T.Data(), &n, &sdim, wr.Data(), wi.Data(), beta.Data(), 0, &n, 0, &n, work.Data(), &bs, 0, &info);
05270 
05271         alpha = Vector<ComplexDouble>(n);
05272         for(int i=0;i<n;i++)
05273         {
05274                 alpha[i] = ComplexDouble(wr[i],wi[i]);
05275         }
05276 
05277         if(info != 0)
05278         {
05279                 if(info < 0)
05280                 {
05281                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
05282                         char s[500];
05283                         sprintf(s, "%i'th argument of DGGES is an illegal value!", -info);
05284                         Utility::RunTimeError(s);
05285                 }
05286                 else if(info > 0)  
05287                 {
05288                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
05289                         char s[500];
05290                         sprintf(s, "Lapack error!");
05291                         Utility::RunTimeError(s);
05292                 }
05293         }
05294 }
05295 
05296 
05297 void Lapack::Gges(Matrix<double>& A, Matrix<double>& B, Vector<ComplexDouble>& alpha, Vector<double>& beta, Matrix<double>& S, Matrix<double>& T, Matrix<double>& Vl, Matrix<double>& Vr)
05298 {
05299         if(A.Rows() != A.Columns() || B.Rows() != B.Columns())
05300         {
05301                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
05302                 Utility::RunTimeError("Matrix is not square!");
05303         }
05304         if(A.Rows() != B.Rows())
05305         {
05306                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
05307                 Utility::RunTimeError("Matrix dimensions does not match!");
05308         }
05309 
05310         
05311         int info;
05312         S = A.Clone();
05313         T = B.Clone();
05314         int n = A.Rows();
05315         Vector<double> wr(n);
05316         Vector<double> wi(n);
05317         beta = Vector<double>(n);
05318 
05319         Vl = Matrix<double>(n,n);
05320         Vr = Matrix<double>(n,n);
05321 
05322         int bs = 39*n+23;
05323         Vector<double> work(bs);
05324 
05325         int sdim;
05326 
05327         DGGES("V", "V", "N", 0, &n, S.Data(), &n, T.Data(), &n, &sdim, wr.Data(), wi.Data(), beta.Data(), Vl.Data(), &n, Vr.Data(), &n, work.Data(), &bs, 0, &info);
05328 
05329         alpha = Vector<ComplexDouble>(n);
05330         for(int i=0;i<n;i++)
05331         {
05332                 alpha[i] = ComplexDouble(wr[i],wi[i]);
05333         }
05334 
05335         if(info != 0)
05336         {
05337                 if(info < 0)
05338                 {
05339                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
05340                         char s[500];
05341                         sprintf(s, "%i'th argument of DGGES is an illegal value!", -info);
05342                         Utility::RunTimeError(s);
05343                 }
05344                 else if(info > 0)  
05345                 {
05346                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
05347                         char s[500];
05348                         sprintf(s, "Lapack error!");
05349                         Utility::RunTimeError(s);
05350                 }
05351         }
05352 }
05353 
05354 
05355 
05356 void Lapack::Gges(Matrix<ComplexFloat>& A, Matrix<ComplexFloat>& B, Vector<ComplexFloat>& alpha, Vector<ComplexFloat>& beta, Matrix<ComplexFloat>& S, Matrix<ComplexFloat>& T)
05357 {
05358         if(A.Rows() != A.Columns() || B.Rows() != B.Columns())
05359         {
05360                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
05361                 Utility::RunTimeError("Matrix is not square!");
05362         }
05363         if(A.Rows() != B.Rows())
05364         {
05365                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
05366                 Utility::RunTimeError("Matrix dimensions does not match!");
05367         }
05368 
05369         
05370         int info;
05371         S = A.Clone();
05372         T = B.Clone();
05373         int n = A.Rows();
05374         alpha = Vector<ComplexFloat>(n);
05375         beta = Vector<ComplexFloat>(n);
05376 
05377         int bs = 33*n;
05378         Vector<ComplexFloat> work(bs);
05379 
05380         Vector<float> rwork(8*n);
05381         int sdim;
05382 
05383         CGGES("N", "N", "N", 0, &n, reinterpret_cast<MKL_Complex8*>(S.Data()), &n, reinterpret_cast<MKL_Complex8*>(T.Data()), &n, &sdim, reinterpret_cast<MKL_Complex8*>(alpha.Data()), reinterpret_cast<MKL_Complex8*>(beta.Data()), 0, &n, 0, &n, reinterpret_cast<MKL_Complex8*>(work.Data()), &bs, rwork.Data(), 0, &info);
05384 
05385 
05386         if(info != 0)
05387         {
05388                 if(info < 0)
05389                 {
05390                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
05391                         char s[500];
05392                         sprintf(s, "%i'th argument of CGGES is an illegal value!", -info);
05393                         Utility::RunTimeError(s);
05394                 }
05395                 else if(info > 0)  
05396                 {
05397                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
05398                         char s[500];
05399                         sprintf(s, "Lapack error!");
05400                         Utility::RunTimeError(s);
05401                 }
05402         }
05403 }
05404 
05405 void Lapack::Gges(Matrix<ComplexFloat>& A, Matrix<ComplexFloat>& B, Vector<ComplexFloat>& alpha, Vector<ComplexFloat>& beta, Matrix<ComplexFloat>& S, Matrix<ComplexFloat>& T, Matrix<ComplexFloat>& Vl, Matrix<ComplexFloat>& Vr)
05406 {
05407         if(A.Rows() != A.Columns() || B.Rows() != B.Columns())
05408         {
05409                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
05410                 Utility::RunTimeError("Matrix is not square!");
05411         }
05412         if(A.Rows() != B.Rows())
05413         {
05414                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
05415                 Utility::RunTimeError("Matrix dimensions does not match!");
05416         }
05417 
05418         
05419         int info;
05420         S = A.Clone();
05421         T = B.Clone();
05422         int n = A.Rows();
05423         alpha = Vector<ComplexFloat>(n);
05424         beta = Vector<ComplexFloat>(n);
05425 
05426         Vl = Matrix<ComplexFloat>(n,n);
05427         Vr = Matrix<ComplexFloat>(n,n);
05428 
05429         int bs = 33*n;
05430         Vector<ComplexFloat> work(bs);
05431 
05432         Vector<float> rwork(8*n);
05433         int sdim;
05434 
05435         CGGES("V", "V", "N", 0, &n, reinterpret_cast<MKL_Complex8*>(S.Data()), &n, reinterpret_cast<MKL_Complex8*>(T.Data()), &n, &sdim, reinterpret_cast<MKL_Complex8*>(alpha.Data()), reinterpret_cast<MKL_Complex8*>(beta.Data()), reinterpret_cast<MKL_Complex8*>(Vl.Data()), &n, reinterpret_cast<MKL_Complex8*>(Vr.Data()), &n, reinterpret_cast<MKL_Complex8*>(work.Data()), &bs, rwork.Data(), 0, &info);
05436 
05437         if(info != 0)
05438         {
05439                 if(info < 0)
05440                 {
05441                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
05442                         char s[500];
05443                         sprintf(s, "%i'th argument of CGGES is an illegal value!", -info);
05444                         Utility::RunTimeError(s);
05445                 }
05446                 else if(info > 0)  
05447                 {
05448                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
05449                         char s[500];
05450                         sprintf(s, "Lapack error!");
05451                         Utility::RunTimeError(s);
05452                 }
05453         }
05454 }
05455 
05456 void Lapack::Gges(Matrix<ComplexDouble>& A, Matrix<ComplexDouble>& B, Vector<ComplexDouble>& alpha, Vector<ComplexDouble>& beta, Matrix<ComplexDouble>& S, Matrix<ComplexDouble>& T)
05457 {
05458         if(A.Rows() != A.Columns() || B.Rows() != B.Columns())
05459         {
05460                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
05461                 Utility::RunTimeError("Matrix is not square!");
05462         }
05463         if(A.Rows() != B.Rows())
05464         {
05465                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
05466                 Utility::RunTimeError("Matrix dimensions does not match!");
05467         }
05468 
05469         
05470         int info;
05471         S = A.Clone();
05472         T = B.Clone();
05473         int n = A.Rows();
05474         alpha = Vector<ComplexDouble>(n);
05475         beta = Vector<ComplexDouble>(n);
05476 
05477         int bs = 33*n;
05478         Vector<ComplexDouble> work(bs);
05479 
05480         Vector<double> rwork(8*n);
05481         int sdim;
05482 
05483         ZGGES("N", "N", "N", 0, &n, reinterpret_cast<MKL_Complex16*>(S.Data()), &n, reinterpret_cast<MKL_Complex16*>(T.Data()), &n, &sdim, reinterpret_cast<MKL_Complex16*>(alpha.Data()), reinterpret_cast<MKL_Complex16*>(beta.Data()), 0, &n, 0, &n, reinterpret_cast<MKL_Complex16*>(work.Data()), &bs, rwork.Data(), 0, &info);
05484 
05485 
05486         if(info != 0)
05487         {
05488                 if(info < 0)
05489                 {
05490                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
05491                         char s[500];
05492                         sprintf(s, "%i'th argument of ZGGES is an illegal value!", -info);
05493                         Utility::RunTimeError(s);
05494                 }
05495                 else if(info > 0)  
05496                 {
05497                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
05498                         char s[500];
05499                         sprintf(s, "Lapack error!");
05500                         Utility::RunTimeError(s);
05501                 }
05502         }
05503 }
05504 
05505 void Lapack::Gges(Matrix<ComplexDouble>& A, Matrix<ComplexDouble>& B, Vector<ComplexDouble>& alpha, Vector<ComplexDouble>& beta, Matrix<ComplexDouble>& S, Matrix<ComplexDouble>& T, Matrix<ComplexDouble>& Vl, Matrix<ComplexDouble>& Vr)
05506 {
05507         if(A.Rows() != A.Columns() || B.Rows() != B.Columns())
05508         {
05509                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
05510                 Utility::RunTimeError("Matrix is not square!");
05511         }
05512         if(A.Rows() != B.Rows())
05513         {
05514                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
05515                 Utility::RunTimeError("Matrix dimensions does not match!");
05516         }
05517 
05518         
05519         int info;
05520         S = A.Clone();
05521         T = B.Clone();
05522         int n = A.Rows();
05523         alpha = Vector<ComplexDouble>(n);
05524         beta = Vector<ComplexDouble>(n);
05525 
05526         Vl = Matrix<ComplexDouble>(n,n);
05527         Vr = Matrix<ComplexDouble>(n,n);
05528 
05529         int bs = 33*n;
05530         Vector<ComplexDouble> work(bs);
05531 
05532         Vector<double> rwork(8*n);
05533         int sdim;
05534 
05535         ZGGES("V", "V", "N", 0, &n, reinterpret_cast<MKL_Complex16*>(S.Data()), &n, reinterpret_cast<MKL_Complex16*>(T.Data()), &n, &sdim, reinterpret_cast<MKL_Complex16*>(alpha.Data()), reinterpret_cast<MKL_Complex16*>(beta.Data()), reinterpret_cast<MKL_Complex16*>(Vl.Data()), &n, reinterpret_cast<MKL_Complex16*>(Vr.Data()), &n, reinterpret_cast<MKL_Complex16*>(work.Data()), &bs, rwork.Data(), 0, &info);
05536 
05537 
05538         if(info != 0)
05539         {
05540                 if(info < 0)
05541                 {
05542                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
05543                         char s[500];
05544                         sprintf(s, "%i'th argument of CGGES is an illegal value!", -info);
05545                         Utility::RunTimeError(s);
05546                 }
05547                 else if(info > 0)  
05548                 {
05549                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
05550                         char s[500];
05551                         sprintf(s, "Lapack error!");
05552                         Utility::RunTimeError(s);
05553                 }
05554         }
05555 }
05556 
05557 
05558 
05559 
05560 
05561 
05562 
05563 // Ggev
05564 
05565 void Lapack::Ggev(Matrix<float>& A, Matrix<float>& B, Vector<ComplexFloat>& alpha, Vector<float>& beta)
05566 {
05567         if(A.Rows() != A.Columns() || B.Rows() != B.Columns())
05568         {
05569                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
05570                 Utility::RunTimeError("Matrix is not square!");
05571         }
05572         if(A.Rows() != B.Rows())
05573         {
05574                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
05575                 Utility::RunTimeError("Matrix dimensions does not match!");
05576         }
05577 
05578         
05579         int info;
05580         Matrix<float> A2 = A.Clone();
05581         Matrix<float> B2 = B.Clone();
05582         int n = A2.Rows();
05583         Vector<float> wr(n);
05584         Vector<float> wi(n);
05585         beta = Vector<float>(n);
05586 
05587         int bs = 8*n+16;
05588         Vector<float> work(bs);
05589 
05590         SGGEV("N", "N", &n, A2.Data(), &n, B2.Data(), &n, wr.Data(), wi.Data(), beta.Data(), 0, &n, 0, &n, work.Data(), &bs, &info);
05591         
05592         alpha = Vector<ComplexFloat>(n);
05593         for(int i=0;i<n;i++)
05594         {
05595                 alpha[i] = ComplexFloat(wr[i],wi[i]);
05596         }
05597 
05598         if(info != 0)
05599         {
05600                 if(info < 0)
05601                 {
05602                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
05603                         char s[500];
05604                         sprintf(s, "%i'th argument of SGGEV is an illegal value!", -info);
05605                         Utility::RunTimeError(s);
05606                 }
05607                 else if(info > 0)  
05608                 {
05609                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
05610                         char s[500];
05611                         sprintf(s, "Lapack error!");
05612                         Utility::RunTimeError(s);
05613                 }
05614         }
05615 }
05616 
05617 void Lapack::Ggev(Matrix<float>& A, Matrix<float>& B, Vector<ComplexFloat>& alpha, Vector<float>& beta, Matrix<ComplexFloat>& Vl, Matrix<ComplexFloat>& Vr)
05618 {
05619         if(A.Rows() != A.Columns() || B.Rows() != B.Columns())
05620         {
05621                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
05622                 Utility::RunTimeError("Matrix is not square!");
05623         }
05624         if(A.Rows() != B.Rows())
05625         {
05626                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
05627                 Utility::RunTimeError("Matrix dimensions does not match!");
05628         }
05629         
05630         int info;
05631         Matrix<float> A2 = A.Clone();
05632         Matrix<float> B2 = B.Clone();
05633         int n = A2.Rows();
05634         Vector<float> wr(n);
05635         Vector<float> wi(n);
05636         beta = Vector<float>(n);
05637 
05638         Vl = Matrix<ComplexFloat>(n,n);
05639         Vr = Matrix<ComplexFloat>(n,n);
05640         Matrix<float> vl(n,n);
05641         Matrix<float> vr(n,n);
05642 
05643         int bs = 8*n+16;
05644         Vector<float> work(bs);
05645 
05646         SGGEV("V", "V", &n, A2.Data(), &n, B2.Data(), &n, wr.Data(), wi.Data(), beta.Data(), vl.Data(), &n, vr.Data(), &n, work.Data(), &bs, &info);
05647 
05648         alpha = Vector<ComplexFloat>(n);
05649         bool conj = false