Seldon user's guide

00001 // Copyright (C) 2003-2009 Marc Duruflé
00002 // Copyright (C) 2001-2009 Vivien Mallet
00003 //
00004 // This file is part of the linear-algebra library Seldon,
00005 // http://seldon.sourceforge.net/.
00006 //
00007 // Seldon is free software; you can redistribute it and/or modify it under the
00008 // terms of the GNU Lesser General Public License as published by the Free
00009 // Software Foundation; either version 2.1 of the License, or (at your option)
00010 // any later version.
00011 //
00012 // Seldon is distributed in the hope that it will be useful, but WITHOUT ANY
00013 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00014 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
00015 // more details.
00016 //
00017 // You should have received a copy of the GNU Lesser General Public License
00018 // along with Seldon. If not, see http://www.gnu.org/licenses/.
00019 
00020 
00021 #ifndef SELDON_FILE_LAPACK_LEAST_SQUARES_CXX
00022 
00023 /*
00024   Functions included in this file:
00025 
00026   xGEQRF   (GetQR, GetLQ)
00027   xGELQF   (GetQR, GetLQ)
00028   xGEQP3   (GetQR_Pivot)
00029   xORGQR   (GetQ_FromQR)
00030   xUNGQR   (GetQ_FromQR)
00031   xUNMQR   (MltQ_FromQR)
00032   xORMQR   (MltQ_FromQR)
00033   xORMQR + xTRSM   (SolveQR)
00034   ZUNMQR + ZTRSM   (SolveQR)
00035   xORMLQ + xTRSM   (SolveQR)
00036   ZUNMLQ + ZTRSM   (SolveQR)
00037   xTRSM + xORMLQ   (SolveLQ)
00038   ZTRSM + ZUNMLQ   (SolveLQ)
00039   xTRSM + xORMQR   (SolveLQ)
00040   ZTRSM + ZUNMQR   (SolveLQ)
00041 */
00042 
00043 namespace Seldon
00044 {
00045 
00046 
00048   // GETQR //
00049 
00050 
00051   /*** ColMajor ***/
00052 
00053 
00054   template<class Prop0, class Allocator0,
00055            class Allocator1>
00056   void GetQR(Matrix<float, Prop0, ColMajor, Allocator0>& A,
00057              Vector<float, VectFull, Allocator1>& tau,
00058              LapackInfo& info = lapack_info)
00059   {
00060     int m = A.GetM();
00061     int n = A.GetN();
00062     int lwork = max(m,n);
00063     Vector<float, VectFull, Allocator1> work(lwork);
00064     tau.Reallocate(min(m, n));
00065     sgeqrf_(&m, &n, A.GetData(), &m, tau.GetData(),
00066             work.GetData(), &lwork, &info.GetInfoRef());
00067   }
00068 
00069 
00070   template<class Prop0, class Allocator0,
00071            class Allocator1>
00072   void GetQR(Matrix<double, Prop0, ColMajor, Allocator0>& A,
00073              Vector<double, VectFull, Allocator1>& tau,
00074              LapackInfo& info = lapack_info)
00075   {
00076     int m = A.GetM();
00077     int n = A.GetN();
00078     int lwork = max(m,n);
00079     Vector<double, VectFull, Allocator1> work(lwork);
00080     tau.Reallocate(min(m, n));
00081     dgeqrf_(&m, &n, A.GetData(), &m, tau.GetData(),
00082             work.GetData(), &lwork, &info.GetInfoRef());
00083   }
00084 
00085 
00086   template<class Prop0, class Allocator0,
00087            class Allocator1>
00088   void GetQR(Matrix<complex<double>, Prop0, ColMajor, Allocator0>& A,
00089              Vector<complex<double>, VectFull, Allocator1>& tau,
00090              LapackInfo& info = lapack_info)
00091   {
00092     int m = A.GetM();
00093     int n = A.GetN();
00094     int lwork = max(m,n);
00095     Vector<complex<double>, VectFull, Allocator1> work(lwork);
00096     tau.Reallocate(min(m, n));
00097     zgeqrf_(&m, &n, A.GetData(), &m, tau.GetData(),
00098             work.GetData(), &lwork, &info.GetInfoRef());
00099   }
00100 
00101 
00102   /*** RowMajor ***/
00103 
00104 
00105   template<class Prop0, class Allocator0,
00106            class Allocator1>
00107   void GetQR(Matrix<float, Prop0, RowMajor, Allocator0>& A,
00108              Vector<float, VectFull, Allocator1>& tau,
00109              LapackInfo& info = lapack_info)
00110   {
00111     int m = A.GetM();
00112     int n = A.GetN();
00113     int lwork = max(m,n);
00114     Vector<float, VectFull, Allocator1> work(lwork);
00115     tau.Reallocate(min(m, n));
00116     // Factorization LQ of A^t.
00117     sgelqf_(&n, &m, A.GetData(), &n, tau.GetData(),
00118             work.GetData(), &lwork, &info.GetInfoRef());
00119   }
00120 
00121 
00122   template<class Prop0, class Allocator0,
00123            class Allocator1>
00124   void GetQR(Matrix<double, Prop0, RowMajor, Allocator0>& A,
00125              Vector<double, VectFull, Allocator1>& tau,
00126              LapackInfo& info = lapack_info)
00127   {
00128     int m = A.GetM();
00129     int n = A.GetN();
00130     int lwork = max(m,n);
00131     Vector<double, VectFull, Allocator1> work(lwork);
00132     tau.Reallocate(min(m, n));
00133     // Factorization LQ of A^t.
00134     dgelqf_(&n, &m, A.GetData(), &n, tau.GetData(),
00135             work.GetData(), &lwork, &info.GetInfoRef());
00136   }
00137 
00138 
00139   template<class Prop0, class Allocator0,
00140            class Allocator1>
00141   void GetQR(Matrix<complex<double>, Prop0, RowMajor, Allocator0>& A,
00142              Vector<complex<double>, VectFull, Allocator1>& tau,
00143              LapackInfo& info = lapack_info)
00144   {
00145     int m = A.GetM();
00146     int n = A.GetN();
00147     int lwork = max(m,n);
00148     Vector<complex<double>, VectFull, Allocator1> work(lwork);
00149     tau.Reallocate(min(m, n));
00150     // Factorization LQ of A^t.
00151     zgelqf_(&n, &m, A.GetData(), &n, tau.GetData(),
00152             work.GetData(), &lwork, &info.GetInfoRef());
00153   }
00154 
00155 
00156   // GETQR //
00158 
00159 
00161   // GETQR_PIVOT //
00162 
00163 
00164   /*** ColMajor ***/
00165 
00166 
00167   template<class Prop0, class Allocator0,
00168            class Allocator1>
00169   void GetQR_Pivot(Matrix<double, Prop0, ColMajor, Allocator0>& A,
00170                    Vector<double, VectFull, Allocator1>& tau,
00171                    Vector<int>& ipivot, LapackInfo& info = lapack_info)
00172   {
00173     int m = A.GetM();
00174     int n = A.GetN();
00175     int lwork = 4 * max(m, n);
00176     ipivot.Fill(0);
00177     Vector<double, VectFull, Allocator1> work(lwork);
00178     tau.Reallocate(min(m, n));
00179     dgeqp3_(&m, &n, A.GetData(), &m, ipivot.GetData(), tau.GetData(),
00180             work.GetData(), &lwork, &info.GetInfoRef());
00181   }
00182 
00183 
00184   // GETQR_PIVOT //
00186 
00187 
00189   // GETQ_FROMQR //
00190 
00191 
00192   /*** ColMajor ***/
00193 
00194   template<class Prop0, class Allocator0,
00195            class Allocator1>
00196   void GetQ_FromQR(Matrix<double, Prop0, ColMajor, Allocator0>& A,
00197                    Vector<double, VectFull, Allocator1>& tau,
00198                    LapackInfo& info = lapack_info)
00199   {
00200     int m = A.GetM();
00201     int n = A.GetN();
00202     int lwork = 2 * max(m, n);
00203     Vector<double, VectFull, Allocator1> work(lwork);
00204     dorgqr_(&m, &m, &n, A.GetData(), &m, tau.GetData(),
00205             work.GetData(), &lwork, &info.GetInfoRef());
00206   }
00207 
00208 
00209   template<class Prop0, class Allocator0,
00210            class Allocator1>
00211   void GetQ_FromQR(Matrix<complex<double>, Prop0, ColMajor, Allocator0>& A,
00212                    Vector<complex<double>, VectFull, Allocator1>& tau,
00213                    LapackInfo& info = lapack_info)
00214   {
00215     int m = A.GetM();
00216     int n = A.GetN();
00217     int lwork = 2 * max(m, n);
00218     Vector<double, VectFull, Allocator1> work(lwork);
00219     zungqr_(&m, &m, &n, A.GetDataVoid(), &m, tau.GetData(),
00220             work.GetData(), &lwork, &info.GetInfoRef());
00221   }
00222 
00223 
00224   template<class Prop0, class Allocator0,
00225            class Allocator1, class Allocator2, class Side, class Trans>
00226   void MltQ_FromQR(const Side& side, const Trans& trans,
00227                    Matrix<complex<double>, Prop0, ColMajor, Allocator0>& A,
00228                    Vector<complex<double>, VectFull, Allocator1>& tau,
00229                    Matrix<complex<double>, Prop0, ColMajor, Allocator2>& C,
00230                    LapackInfo& info = lapack_info)
00231   {
00232     int m = A.GetM();
00233     int n = A.GetN();
00234     int lwork = max(m, n);
00235     Vector<double, VectFull, Allocator1> work(lwork);
00236     char side_ = side.Char();
00237     char trans_ = trans.Char();
00238     int k = m;
00239     if (side_ == 'R')
00240       k = n;
00241 
00242     zunmqr_(&side, &trans, &m, &n, &k, A.GetDataVoid(), &m, tau.GetDataVoid(),
00243             C.GetDataVoid(), &m, work.GetData(), &lwork,
00244             &info.GetInfoRef());
00245   }
00246 
00247 
00248   // GETQ_FROMQR //
00250 
00251 
00253   // GETLQ //
00254 
00255 
00256   /*** ColMajor ***/
00257 
00258 
00259   template<class Prop0, class Allocator0,
00260            class Allocator1>
00261   void GetLQ(Matrix<float, Prop0, ColMajor, Allocator0>& A,
00262              Vector<float, VectFull, Allocator1>& tau,
00263              LapackInfo& info = lapack_info)
00264   {
00265     int m = A.GetM();
00266     int n = A.GetN();
00267     int lwork = max(m,n);
00268     Vector<float, VectFull, Allocator1> work(lwork);
00269     tau.Reallocate(min(m, n));
00270     sgelqf_(&m, &n, A.GetData(), &m, tau.GetData(),
00271             work.GetData(), &lwork, &info.GetInfoRef());
00272   }
00273 
00274 
00275   template<class Prop0, class Allocator0,
00276            class Allocator1>
00277   void GetLQ(Matrix<double, Prop0, ColMajor, Allocator0>& A,
00278              Vector<double, VectFull, Allocator1>& tau,
00279              LapackInfo& info = lapack_info)
00280   {
00281     int m = A.GetM();
00282     int n = A.GetN();
00283     int lwork = max(m,n);
00284     Vector<double, VectFull, Allocator1> work(lwork);
00285     tau.Reallocate(min(m, n));
00286     dgelqf_(&m, &n, A.GetData(), &m, tau.GetData(),
00287             work.GetData(), &lwork, &info.GetInfoRef());
00288   }
00289 
00290 
00291   template<class Prop0, class Allocator0,
00292            class Allocator1>
00293   void GetLQ(Matrix<complex<double>, Prop0, ColMajor, Allocator0>& A,
00294              Vector<complex<double>, VectFull, Allocator1>& tau,
00295              LapackInfo& info = lapack_info)
00296   {
00297     int m = A.GetM();
00298     int n = A.GetN();
00299     int lwork = max(m,n);
00300     Vector<complex<double>, VectFull, Allocator1> work(lwork);
00301     tau.Reallocate(min(m, n));
00302     zgelqf_(&m, &n, A.GetData(), &m, tau.GetData(),
00303             work.GetData(), &lwork, &info.GetInfoRef());
00304   }
00305 
00306 
00307   /*** RowMajor ***/
00308 
00309 
00310   template<class Prop0, class Allocator0,
00311            class Allocator1>
00312   void GetLQ(Matrix<float, Prop0, RowMajor, Allocator0>& A,
00313              Vector<float, VectFull, Allocator1>& tau,
00314              LapackInfo& info = lapack_info)
00315   {
00316     int m = A.GetM();
00317     int n = A.GetN();
00318     int lwork = max(m,n);
00319     Vector<float, VectFull, Allocator1> work(lwork);
00320     tau.Reallocate(min(m, n));
00321     // Factorization QR of A^t.
00322     sgeqrf_(&n, &m, A.GetData(), &n, tau.GetData(),
00323             work.GetData(), &lwork, &info.GetInfoRef());
00324   }
00325 
00326 
00327   template<class Prop0, class Allocator0,
00328            class Allocator1>
00329   void GetLQ(Matrix<double, Prop0, RowMajor, Allocator0>& A,
00330              Vector<double, VectFull, Allocator1>& tau,
00331              LapackInfo& info = lapack_info)
00332   {
00333     int m = A.GetM();
00334     int n = A.GetN();
00335     int lwork = max(m,n);
00336     Vector<double, VectFull, Allocator1> work(lwork);
00337     tau.Reallocate(min(m, n));
00338     // Factorization LQ of A^t.
00339     dgeqrf_(&n, &m, A.GetData(), &n, tau.GetData(),
00340             work.GetData(), &lwork, &info.GetInfoRef());
00341   }
00342 
00343 
00344   template<class Prop0, class Allocator0,
00345            class Allocator1>
00346   void GetLQ(Matrix<complex<double>, Prop0, RowMajor, Allocator0>& A,
00347              Vector<complex<double>, VectFull, Allocator1>& tau,
00348              LapackInfo& info = lapack_info)
00349   {
00350     int m = A.GetM();
00351     int n = A.GetN();
00352     int lwork = max(m,n);
00353     Vector<complex<double>, VectFull, Allocator1> work(lwork);
00354     tau.Reallocate(min(m, n));
00355     // Factorization LQ of A^t.
00356     zgeqrf_(&n, &m, A.GetData(), &n, tau.GetData(),
00357             work.GetData(), &lwork, &info.GetInfoRef());
00358   }
00359 
00360 
00361   // GETLQ //
00363 
00364 
00366   // MLTQ_FROMQR //
00367 
00368 
00369   /*** ColMajor ***/
00370 
00371 
00372   template<class Prop0, class Allocator0,
00373            class Allocator1, class Allocator2, class IsTranspose>
00374   void MltQ_FromQR(Matrix<double, Prop0, ColMajor, Allocator0>& A,
00375                    const IsTranspose& trans,
00376                    Vector<double, VectFull, Allocator1>& tau,
00377                    Vector<double, VectFull, Allocator2>& b,
00378                    LapackInfo& info = lapack_info)
00379   {
00380     int m = b.GetM();
00381     int n = 1;
00382     int k = tau.GetM();
00383     int lwork = max(m,n);
00384     Vector<double, VectFull, Allocator1> work(lwork);
00385     char side('L');
00386     char trans_(trans);
00387     dormqr_(&side, &trans_, &m, &n, &k, A.GetData(), &m, tau.GetData(),
00388             b.GetData(), &m, work.GetData(), &lwork,
00389             &info.GetInfoRef());
00390   }
00391 
00392 
00393   // MLTQ_FROMQR //
00395 
00396 
00398   // SOLVEQR //
00399 
00400 
00401   /*** ColMajor ***/
00402 
00403 
00404   template <class Prop0, class Allocator0,
00405             class Allocator1,class Allocator2>
00406   void SolveQR(const Matrix<float, Prop0, ColMajor, Allocator0>& A,
00407                const Vector<float, VectFull, Allocator1>& tau,
00408                Vector<float, VectFull, Allocator2>& b,
00409                LapackInfo& info = lapack_info)
00410   {
00411     int m = A.GetM();
00412     int n = A.GetN();
00413     int k = tau.GetM();
00414     int nrhs = 1, nb = b.GetM();
00415     char side('L');
00416     char trans('T');
00417     int lwork = max(m, n);
00418     Vector<float, VectFull, Allocator1> work(lwork);
00419     // Computes Q^t b.
00420     sormqr_(&side, &trans, &m, &nrhs, &k, A.GetData(),
00421             &m, tau.GetData(), b.GetData(),
00422             &m, work.GetData(), &lwork, &info.GetInfoRef());
00423 
00424     b.Resize(n);
00425     for (int i = nb; i < n; i++)
00426       b(i) = 0;
00427 
00428     // Then solves R x = Q^t b.
00429     float alpha(1);
00430     cblas_strsm(CblasColMajor, CblasLeft, CblasUpper, CblasNoTrans,
00431                 CblasNonUnit, b.GetM(), nrhs,
00432                 alpha, A.GetData(), A.GetM(), b.GetData(), b.GetM());
00433   }
00434 
00435 
00436   template <class Prop0, class Allocator0,
00437             class Allocator1,class Allocator2>
00438   void SolveQR(const Matrix<double, Prop0, ColMajor, Allocator0>& A,
00439                const Vector<double, VectFull, Allocator1>& tau,
00440                Vector<double, VectFull, Allocator2>& b,
00441                LapackInfo& info = lapack_info)
00442   {
00443     int m = A.GetM();
00444     int n = A.GetN();
00445     int k = tau.GetM();
00446     int nrhs = 1, nb = b.GetM();
00447     char side('L');
00448     char trans('T');
00449     int lwork = max(m, n);
00450     Vector<double, VectFull, Allocator1> work(lwork);
00451     // Computes Q^t b.
00452     dormqr_(&side, &trans, &lwork, &nrhs, &k, A.GetData(),
00453             &m, tau.GetData(), b.GetData(),
00454             &lwork, work.GetData(), &lwork, &info.GetInfoRef());
00455 
00456     b.Resize(n);
00457     for (int i = nb; i < n; i++)
00458       b(i) = 0;
00459 
00460     // Then solves R x = Q^t b.
00461     double alpha(1);
00462     cblas_dtrsm(CblasColMajor, CblasLeft, CblasUpper, CblasNoTrans,
00463                 CblasNonUnit, b.GetM(), nrhs,
00464                 alpha, A.GetData(), A.GetM(), b.GetData(), b.GetM());
00465   }
00466 
00467 
00468   template <class Prop0, class Allocator0,
00469             class Allocator1,class Allocator2>
00470   void SolveQR(const Matrix<complex<double>, Prop0, ColMajor, Allocator0>& A,
00471                const Vector<complex<double>, VectFull, Allocator1>& tau,
00472                Vector<complex<double>, VectFull, Allocator2>& b,
00473                LapackInfo& info = lapack_info)
00474   {
00475     int m = A.GetM();
00476     int n = A.GetN();
00477     int k = tau.GetM();
00478     int nrhs = 1, nb = b.GetM();
00479     char side('L');
00480     char trans('C');
00481     int lwork = max(m, n);
00482     Vector<complex<double>, VectFull, Allocator1> work(lwork);
00483     // Computes Q^t b.
00484     zunmqr_(&side, &trans, &m, &nrhs, &k, A.GetData(),
00485             &m, tau.GetData(), b.GetData(),
00486             &m, work.GetData(), &lwork, &info.GetInfoRef());
00487 
00488     b.Resize(n);
00489     for (int i = nb; i < n; i++)
00490       b(i) = 0;
00491 
00492     // Then solves R x = Q^t b.
00493     complex<double> alpha(1);
00494     cblas_ztrsm(CblasColMajor, CblasLeft, CblasUpper, CblasNoTrans,
00495                 CblasNonUnit, b.GetM(), nrhs,
00496                 &alpha, A.GetData(), A.GetM(), b.GetData(), b.GetM());
00497   }
00498 
00499 
00500   /*** RowMajor ***/
00501 
00502 
00503   template <class Prop0, class Allocator0,
00504             class Allocator1,class Allocator2>
00505   void SolveQR(const Matrix<float, Prop0, RowMajor, Allocator0>& A,
00506                const Vector<float, VectFull, Allocator1>& tau,
00507                Vector<float, VectFull, Allocator2>& b,
00508                LapackInfo& info = lapack_info)
00509   {
00510     int m = A.GetM();
00511     int n = A.GetN();
00512     int k = tau.GetM();
00513     int nrhs = 1, nb = b.GetM();
00514     char side('L');
00515     char trans('N');
00516     int lwork = max(m, n);
00517     Vector<float, VectFull, Allocator1> work(lwork);
00518     // Computes Q b.
00519     sormlq_(&side, &trans, &m, &nrhs, &k, A.GetData(),
00520             &n, tau.GetData(), b.GetData(),
00521             &m, work.GetData(), &lwork, &info.GetInfoRef());
00522 
00523     b.Resize(n);
00524     for (int i = nb; i < n; i++)
00525       b(i) = 0;
00526 
00527     // Solves L^t y = b.
00528     float alpha(1);
00529     cblas_strsm(CblasColMajor, CblasLeft, CblasLower, CblasTrans,
00530                 CblasNonUnit, n, nrhs,
00531                 alpha, A.GetData(), n, b.GetData(), b.GetM());
00532   }
00533 
00534 
00535   template <class Prop0, class Allocator0,
00536             class Allocator1,class Allocator2>
00537   void SolveQR(const Matrix<double, Prop0, RowMajor, Allocator0>& A,
00538                const Vector<double, VectFull, Allocator1>& tau,
00539                Vector<double, VectFull, Allocator2>& b,
00540                LapackInfo& info = lapack_info)
00541   {
00542     int m = A.GetM();
00543     int n = A.GetN();
00544     int k = tau.GetM();
00545     int nrhs = 1, nb = b.GetM();
00546     char side('L');
00547     char trans('N');
00548     int lwork = max(m, n);
00549     Vector<double, VectFull, Allocator1> work(lwork);
00550     // Computes Q b.
00551     dormlq_(&side, &trans, &m, &nrhs, &k, A.GetData(),
00552             &n, tau.GetData(), b.GetData(),
00553             &m, work.GetData(), &lwork, &info.GetInfoRef());
00554 
00555     b.Resize(n);
00556     for (int i = nb; i < n; i++)
00557       b(i) = 0;
00558 
00559     // Solves L^t y = b.
00560     double alpha(1);
00561     cblas_dtrsm(CblasColMajor, CblasLeft, CblasLower, CblasTrans,
00562                 CblasNonUnit, n, nrhs,
00563                 alpha, A.GetData(), n, b.GetData(), b.GetM());
00564   }
00565 
00566 
00567   template <class Prop0, class Allocator0,
00568             class Allocator1,class Allocator2>
00569   void SolveQR(const Matrix<complex<double>, Prop0, RowMajor, Allocator0>& A,
00570                const Vector<complex<double>, VectFull, Allocator1>& tau,
00571                Vector<complex<double>, VectFull, Allocator2>& b,
00572                LapackInfo& info = lapack_info)
00573   {
00574     int m = A.GetM();
00575     int n = A.GetN();
00576     int k = tau.GetM();
00577     int nrhs = 1, nb = b.GetM();
00578     char side('L');
00579     char trans('N');
00580     int lwork = max(m, n);
00581     Vector<complex<double>, VectFull, Allocator1> work(lwork);
00582     // Computes Q b.
00583     zunmlq_(&side, &trans, &m, &nrhs, &k, A.GetData(),
00584             &n, tau.GetData(), b.GetData(),
00585             &m, work.GetData(), &lwork, &info.GetInfoRef());
00586 
00587     b.Resize(n);
00588     for (int i = nb; i < n; i++)
00589       b(i) = 0;
00590 
00591     // Solves L^t y = b.
00592     complex<double> alpha(1);
00593     cblas_ztrsm(CblasColMajor, CblasLeft, CblasLower, CblasTrans,
00594                 CblasNonUnit, n, nrhs,
00595                 &alpha, A.GetData(), n, b.GetData(), b.GetM());
00596   }
00597 
00598 
00599   // SOLVEQR //
00601 
00602 
00604   // SOLVELQ //
00605 
00606 
00607   /*** ColMajor ***/
00608 
00609 
00610   template <class Prop0, class Allocator0,
00611             class Allocator1,class Allocator2>
00612   void SolveLQ(const Matrix<float, Prop0, ColMajor, Allocator0>& A,
00613                const Vector<float, VectFull, Allocator1>& tau,
00614                Vector<float, VectFull, Allocator2>& b,
00615                LapackInfo& info = lapack_info)
00616   {
00617     int m = A.GetM();
00618     int n = A.GetN();
00619     int k = tau.GetM();
00620     int nrhs = 1, nb = b.GetM();
00621     char side('L');
00622     char trans('T');
00623     int lwork = max(m, n);
00624     Vector<float, VectFull, Allocator1> work(lwork);
00625     // Solves L y = b.
00626     float alpha(1);
00627     cblas_strsm(CblasColMajor, CblasLeft, CblasLower, CblasNoTrans,
00628                 CblasNonUnit, m, nrhs,
00629                 alpha, A.GetData(), m, b.GetData(), b.GetM());
00630 
00631     b.Resize(n);
00632     for (int i = nb; i < n; i++)
00633       b(i) = 0;
00634 
00635     // Computes Q^t b.
00636     sormlq_(&side, &trans, &n, &nrhs, &k, A.GetData(),
00637             &m, tau.GetData(), b.GetData(),
00638             &n, work.GetData(), &lwork, &info.GetInfoRef());
00639   }
00640 
00641 
00642   template <class Prop0, class Allocator0,
00643             class Allocator1,class Allocator2>
00644   void SolveLQ(const Matrix<double, Prop0, ColMajor, Allocator0>& A,
00645                const Vector<double, VectFull, Allocator1>& tau,
00646                Vector<double, VectFull, Allocator2>& b,
00647                LapackInfo& info = lapack_info)
00648   {
00649     int m = A.GetM();
00650     int n = A.GetN();
00651     int k = tau.GetM();
00652     int nrhs = 1, nb = b.GetM();
00653     char side('L');
00654     char trans('T');
00655     int lwork = max(m, n);
00656     Vector<double, VectFull, Allocator1> work(lwork);
00657     // Solves L y = b.
00658     double alpha(1);
00659     cblas_dtrsm(CblasColMajor, CblasLeft, CblasLower, CblasNoTrans,
00660                 CblasNonUnit, m, nrhs,
00661                 alpha, A.GetData(), m, b.GetData(), b.GetM());
00662 
00663     b.Resize(n);
00664     for (int i = nb; i < n; i++)
00665       b(i) = 0;
00666 
00667     // Computes Q^t b.
00668     dormlq_(&side, &trans, &n, &nrhs, &k, A.GetData(),
00669             &m, tau.GetData(), b.GetData(),
00670             &n, work.GetData(), &lwork, &info.GetInfoRef());
00671   }
00672 
00673 
00674   template <class Prop0, class Allocator0,
00675             class Allocator1,class Allocator2>
00676   void SolveLQ(const Matrix<complex<double>, Prop0, ColMajor, Allocator0>& A,
00677                const Vector<complex<double>, VectFull, Allocator1>& tau,
00678                Vector<complex<double>, VectFull, Allocator2>& b,
00679                LapackInfo& info = lapack_info)
00680   {
00681     int m = A.GetM();
00682     int n = A.GetN();
00683     int k = tau.GetM();
00684     int nrhs = 1, nb = b.GetM();
00685     char side('L');
00686     char trans('C');
00687     int lwork = max(m, n);
00688     Vector<complex<double>, VectFull, Allocator1> work(lwork);
00689     // Solve L y = b.
00690     complex<double> alpha(1);
00691     cblas_ztrsm(CblasColMajor, CblasLeft, CblasLower, CblasNoTrans,
00692                 CblasNonUnit, m, nrhs,
00693                 &alpha, A.GetData(), m, b.GetData(), b.GetM());
00694 
00695     b.Resize(n);
00696     for (int i = nb; i < n; i++)
00697       b(i) = 0;
00698 
00699     // Computes Q^t.
00700     zunmlq_(&side, &trans, &n, &nrhs, &k, A.GetData(),
00701             &m, tau.GetData(), b.GetData(),
00702             &n, work.GetData(), &lwork, &info.GetInfoRef());
00703   }
00704 
00705 
00706   /*** RowMajor ***/
00707 
00708 
00709   template <class Prop0, class Allocator0,
00710             class Allocator1,class Allocator2>
00711   void SolveLQ(const Matrix<float, Prop0, RowMajor, Allocator0>& A,
00712                const Vector<float, VectFull, Allocator1>& tau,
00713                Vector<float, VectFull, Allocator2>& b,
00714                LapackInfo& info = lapack_info)
00715   {
00716     int m = A.GetM();
00717     int n = A.GetN();
00718     int k = tau.GetM();
00719     int nrhs = 1, nb = b.GetM();
00720     char side('L');
00721     char trans('N');
00722     int lwork = max(m, n);
00723     Vector<float, VectFull, Allocator1> work(lwork);
00724     // Solves R^t x = b.
00725     float alpha(1);
00726     cblas_strsm(CblasColMajor, CblasLeft, CblasUpper, CblasTrans,
00727                 CblasNonUnit, b.GetM(), nrhs,
00728                 alpha, A.GetData(), A.GetN(), b.GetData(), b.GetM());
00729 
00730     b.Resize(n);
00731     for (int i = nb; i < n; i++)
00732       b(i) = 0;
00733 
00734     // Multiplies by Q.
00735     sormqr_(&side, &trans, &n, &nrhs, &k, A.GetData(),
00736             &n, tau.GetData(), b.GetData(),
00737             &n, work.GetData(), &lwork, &info.GetInfoRef());
00738   }
00739 
00740 
00741   template <class Prop0, class Allocator0,
00742             class Allocator1,class Allocator2>
00743   void SolveLQ(const Matrix<double, Prop0, RowMajor, Allocator0>& A,
00744                const Vector<double, VectFull, Allocator1>& tau,
00745                Vector<double, VectFull, Allocator2>& b,
00746                LapackInfo& info = lapack_info)
00747   {
00748     int m = A.GetM();
00749     int n = A.GetN();
00750     int k = tau.GetM();
00751     int nrhs = 1, nb = b.GetM();
00752     char side('L');
00753     char trans('N');
00754     int lwork = max(m, n);
00755     Vector<double, VectFull, Allocator1> work(lwork);
00756     // Solves R^t x = b.
00757     double alpha(1);
00758     cblas_dtrsm(CblasColMajor, CblasLeft, CblasUpper, CblasTrans,
00759                 CblasNonUnit, b.GetM(), nrhs,
00760                 alpha, A.GetData(), A.GetN(), b.GetData(), b.GetM());
00761 
00762     b.Resize(n);
00763     for (int i = nb; i < n; i++)
00764       b(i) = 0;
00765 
00766     // Multiplies by Q.
00767     dormqr_(&side, &trans, &n, &nrhs, &k, A.GetData(),
00768             &n, tau.GetData(), b.GetData(),
00769             &n, work.GetData(), &lwork, &info.GetInfoRef());
00770   }
00771 
00772 
00773   template <class Prop0, class Allocator0,
00774             class Allocator1,class Allocator2>
00775   void SolveLQ(const Matrix<complex<double>, Prop0, RowMajor, Allocator0>& A,
00776                const Vector<complex<double>, VectFull, Allocator1>& tau,
00777                Vector<complex<double>, VectFull, Allocator2>& b,
00778                LapackInfo& info = lapack_info)
00779   {
00780     int m = A.GetM();
00781     int n = A.GetN();
00782     int k = tau.GetM();
00783     int nrhs = 1, nb = b.GetM();
00784     char side('L');
00785     char trans('C');
00786     int lwork = max(m, n);
00787     Vector<complex<double>, VectFull, Allocator1> work(lwork);
00788     // Solves R^t x = b.
00789     complex<double> alpha(1);
00790     cblas_ztrsm(CblasColMajor, CblasLeft, CblasUpper, CblasTrans,
00791                 CblasNonUnit, b.GetM(), nrhs,
00792                 &alpha, A.GetData(), A.GetN(), b.GetData(), b.GetM());
00793 
00794     b.Resize(n);
00795     for (int i = nb; i < n; i++)
00796       b(i) = 0;
00797 
00798     // Computes Q b.
00799     zunmqr_(&side, &trans, &n, &nrhs, &k, A.GetData(),
00800             &n, tau.GetData(), b.GetData(),
00801             &n, work.GetData(), &lwork, &info.GetInfoRef());
00802   }
00803 
00804 
00805   // SOLVELQ //
00807 
00808 
00809 } // end namespace
00810 
00811 #define SELDON_FILE_LAPACK_LEAST_SQUARES_CXX
00812 #endif
00813
All Classes Namespaces Files Functions Variables Typedefs
computation/interfaces/Lapack_LeastSquares.cxx
