lu.cpp
Go to the documentation of this file.
00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2010-2011 Gael Guennebaud <gael.guennebaud@inria.fr>
00005 //
00006 // Eigen is free software; you can redistribute it and/or
00007 // modify it under the terms of the GNU Lesser General Public
00008 // License as published by the Free Software Foundation; either
00009 // version 3 of the License, or (at your option) any later version.
00010 //
00011 // Alternatively, you can redistribute it and/or
00012 // modify it under the terms of the GNU General Public License as
00013 // published by the Free Software Foundation; either version 2 of
00014 // the License, or (at your option) any later version.
00015 //
00016 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
00017 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00018 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
00019 // GNU General Public License for more details.
00020 //
00021 // You should have received a copy of the GNU Lesser General Public
00022 // License and a copy of the GNU General Public License along with
00023 // Eigen. If not, see <http://www.gnu.org/licenses/>.
00024 
00025 #include "common.h"
00026 #include <Eigen/LU>
00027 
00028 // computes an LU factorization of a general M-by-N matrix A using partial pivoting with row interchanges
00029 EIGEN_LAPACK_FUNC(getrf,(int *m, int *n, RealScalar *pa, int *lda, int *ipiv, int *info))
00030 {
00031   *info = 0;
00032         if(*m<0)                  *info = -1;
00033   else  if(*n<0)                  *info = -2;
00034   else  if(*lda<std::max(1,*m))   *info = -4;
00035   if(*info!=0)
00036   {
00037     int e = -*info;
00038     return xerbla_(SCALAR_SUFFIX_UP"GETRF", &e, 6);
00039   }
00040 
00041   if(*m==0 || *n==0)
00042     return 0;
00043 
00044   Scalar* a = reinterpret_cast<Scalar*>(pa);
00045   int nb_transpositions;
00046   int ret = Eigen::internal::partial_lu_impl<Scalar,ColMajor,int>
00047                  ::blocked_lu(*m, *n, a, *lda, ipiv, nb_transpositions);
00048 
00049   for(int i=0; i<std::min(*m,*n); ++i)
00050     ipiv[i]++;
00051 
00052   if(ret>=0)
00053     *info = ret+1;
00054 
00055   return 0;
00056 }
00057 
00058 //GETRS solves a system of linear equations
00059 //    A * X = B  or  A' * X = B
00060 //  with a general N-by-N matrix A using the LU factorization computed  by GETRF
00061 EIGEN_LAPACK_FUNC(getrs,(char *trans, int *n, int *nrhs, RealScalar *pa, int *lda, int *ipiv, RealScalar *pb, int *ldb, int *info))
00062 {
00063   *info = 0;
00064         if(OP(*trans)==INVALID)  *info = -1;
00065   else  if(*n<0)                 *info = -2;
00066   else  if(*nrhs<0)              *info = -3;
00067   else  if(*lda<std::max(1,*n))  *info = -5;
00068   else  if(*ldb<std::max(1,*n))  *info = -8;
00069   if(*info!=0)
00070   {
00071     int e = -*info;
00072     return xerbla_(SCALAR_SUFFIX_UP"GETRS", &e, 6);
00073   }
00074 
00075   Scalar* a = reinterpret_cast<Scalar*>(pa);
00076   Scalar* b = reinterpret_cast<Scalar*>(pb);
00077   MatrixType lu(a,*n,*n,*lda);
00078   MatrixType B(b,*n,*nrhs,*ldb);
00079 
00080   for(int i=0; i<*n; ++i)
00081     ipiv[i]--;
00082   if(OP(*trans)==NOTR)
00083   {
00084     B = PivotsType(ipiv,*n) * B;
00085     lu.triangularView<UnitLower>().solveInPlace(B);
00086     lu.triangularView<Upper>().solveInPlace(B);
00087   }
00088   else if(OP(*trans)==TR)
00089   {
00090     lu.triangularView<Upper>().transpose().solveInPlace(B);
00091     lu.triangularView<UnitLower>().transpose().solveInPlace(B);
00092     B = PivotsType(ipiv,*n).transpose() * B;
00093   }
00094   else if(OP(*trans)==ADJ)
00095   {
00096     lu.triangularView<Upper>().adjoint().solveInPlace(B);
00097     lu.triangularView<UnitLower>().adjoint().solveInPlace(B);
00098     B = PivotsType(ipiv,*n).transpose() * B;
00099   }
00100   for(int i=0; i<*n; ++i)
00101     ipiv[i]++;
00102 
00103   return 0;
00104 }


libicr
Author(s): Robert Krug
autogenerated on Mon Jan 6 2014 11:32:58