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