Home | History | Annotate | Download | only in lapack 
      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 (at) 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