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