1*bf2c3715SXin Li // This file is part of Eigen, a lightweight C++ template library 2*bf2c3715SXin Li // for linear algebra. 3*bf2c3715SXin Li // 4*bf2c3715SXin Li // Copyright (C) 2010-2011 Gael Guennebaud <[email protected]> 5*bf2c3715SXin Li // 6*bf2c3715SXin Li // This Source Code Form is subject to the terms of the Mozilla 7*bf2c3715SXin Li // Public License v. 2.0. If a copy of the MPL was not distributed 8*bf2c3715SXin Li // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 9*bf2c3715SXin Li 10*bf2c3715SXin Li #include "common.h" 11*bf2c3715SXin Li #include <Eigen/LU> 12*bf2c3715SXin Li 13*bf2c3715SXin Li // computes an LU factorization of a general M-by-N matrix A using partial pivoting with row interchanges 14*bf2c3715SXin Li EIGEN_LAPACK_FUNC(getrf,(int *m, int *n, RealScalar *pa, int *lda, int *ipiv, int *info)) 15*bf2c3715SXin Li { 16*bf2c3715SXin Li *info = 0; 17*bf2c3715SXin Li if(*m<0) *info = -1; 18*bf2c3715SXin Li else if(*n<0) *info = -2; 19*bf2c3715SXin Li else if(*lda<std::max(1,*m)) *info = -4; 20*bf2c3715SXin Li if(*info!=0) 21*bf2c3715SXin Li { 22*bf2c3715SXin Li int e = -*info; 23*bf2c3715SXin Li return xerbla_(SCALAR_SUFFIX_UP"GETRF", &e, 6); 24*bf2c3715SXin Li } 25*bf2c3715SXin Li 26*bf2c3715SXin Li if(*m==0 || *n==0) 27*bf2c3715SXin Li return 0; 28*bf2c3715SXin Li 29*bf2c3715SXin Li Scalar* a = reinterpret_cast<Scalar*>(pa); 30*bf2c3715SXin Li int nb_transpositions; 31*bf2c3715SXin Li int ret = int(Eigen::internal::partial_lu_impl<Scalar,ColMajor,int> 32*bf2c3715SXin Li ::blocked_lu(*m, *n, a, *lda, ipiv, nb_transpositions)); 33*bf2c3715SXin Li 34*bf2c3715SXin Li for(int i=0; i<std::min(*m,*n); ++i) 35*bf2c3715SXin Li ipiv[i]++; 36*bf2c3715SXin Li 37*bf2c3715SXin Li if(ret>=0) 38*bf2c3715SXin Li *info = ret+1; 39*bf2c3715SXin Li 40*bf2c3715SXin Li return 0; 41*bf2c3715SXin Li } 42*bf2c3715SXin Li 43*bf2c3715SXin Li //GETRS solves a system of linear equations 44*bf2c3715SXin Li // A * X = B or A' * X = B 45*bf2c3715SXin Li // with a general N-by-N matrix A using the LU factorization computed by GETRF 46*bf2c3715SXin Li EIGEN_LAPACK_FUNC(getrs,(char *trans, int *n, int *nrhs, RealScalar *pa, int *lda, int *ipiv, RealScalar *pb, int *ldb, int *info)) 47*bf2c3715SXin Li { 48*bf2c3715SXin Li *info = 0; 49*bf2c3715SXin Li if(OP(*trans)==INVALID) *info = -1; 50*bf2c3715SXin Li else if(*n<0) *info = -2; 51*bf2c3715SXin Li else if(*nrhs<0) *info = -3; 52*bf2c3715SXin Li else if(*lda<std::max(1,*n)) *info = -5; 53*bf2c3715SXin Li else if(*ldb<std::max(1,*n)) *info = -8; 54*bf2c3715SXin Li if(*info!=0) 55*bf2c3715SXin Li { 56*bf2c3715SXin Li int e = -*info; 57*bf2c3715SXin Li return xerbla_(SCALAR_SUFFIX_UP"GETRS", &e, 6); 58*bf2c3715SXin Li } 59*bf2c3715SXin Li 60*bf2c3715SXin Li Scalar* a = reinterpret_cast<Scalar*>(pa); 61*bf2c3715SXin Li Scalar* b = reinterpret_cast<Scalar*>(pb); 62*bf2c3715SXin Li MatrixType lu(a,*n,*n,*lda); 63*bf2c3715SXin Li MatrixType B(b,*n,*nrhs,*ldb); 64*bf2c3715SXin Li 65*bf2c3715SXin Li for(int i=0; i<*n; ++i) 66*bf2c3715SXin Li ipiv[i]--; 67*bf2c3715SXin Li if(OP(*trans)==NOTR) 68*bf2c3715SXin Li { 69*bf2c3715SXin Li B = PivotsType(ipiv,*n) * B; 70*bf2c3715SXin Li lu.triangularView<UnitLower>().solveInPlace(B); 71*bf2c3715SXin Li lu.triangularView<Upper>().solveInPlace(B); 72*bf2c3715SXin Li } 73*bf2c3715SXin Li else if(OP(*trans)==TR) 74*bf2c3715SXin Li { 75*bf2c3715SXin Li lu.triangularView<Upper>().transpose().solveInPlace(B); 76*bf2c3715SXin Li lu.triangularView<UnitLower>().transpose().solveInPlace(B); 77*bf2c3715SXin Li B = PivotsType(ipiv,*n).transpose() * B; 78*bf2c3715SXin Li } 79*bf2c3715SXin Li else if(OP(*trans)==ADJ) 80*bf2c3715SXin Li { 81*bf2c3715SXin Li lu.triangularView<Upper>().adjoint().solveInPlace(B); 82*bf2c3715SXin Li lu.triangularView<UnitLower>().adjoint().solveInPlace(B); 83*bf2c3715SXin Li B = PivotsType(ipiv,*n).transpose() * B; 84*bf2c3715SXin Li } 85*bf2c3715SXin Li for(int i=0; i<*n; ++i) 86*bf2c3715SXin Li ipiv[i]++; 87*bf2c3715SXin Li 88*bf2c3715SXin Li return 0; 89*bf2c3715SXin Li } 90