Main Page | Namespace List | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members | Search for 
D: / cimpl1_beta / include

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.