xref: /aosp_15_r20/external/eigen/lapack/lu.cpp (revision bf2c37156dfe67e5dfebd6d394bad8b2ab5804d4)
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