3rdparty/Eigen/lapack/cholesky.cpp
Go to the documentation of this file.
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 "lapack_common.h"
11 #include <Eigen/Cholesky>
12 
13 // POTRF computes the Cholesky factorization of a real symmetric positive definite matrix A.
14 EIGEN_LAPACK_FUNC(potrf,(char* uplo, int *n, RealScalar *pa, int *lda, int *info))
15 {
16  *info = 0;
17  if(UPLO(*uplo)==INVALID) *info = -1;
18  else if(*n<0) *info = -2;
19  else if(*lda<std::max(1,*n)) *info = -4;
20  if(*info!=0)
21  {
22  int e = -*info;
23  return xerbla_(SCALAR_SUFFIX_UP"POTRF", &e, 6);
24  }
25 
26  Scalar* a = reinterpret_cast<Scalar*>(pa);
27  MatrixType A(a,*n,*n,*lda);
28  int ret;
29  if(UPLO(*uplo)==UP) ret = int(internal::llt_inplace<Scalar, Upper>::blocked(A));
30  else ret = int(internal::llt_inplace<Scalar, Lower>::blocked(A));
31 
32  if(ret>=0)
33  *info = ret+1;
34 
35  return 0;
36 }
37 
38 // POTRS solves a system of linear equations A*X = B with a symmetric
39 // positive definite matrix A using the Cholesky factorization
40 // A = U**T*U or A = L*L**T computed by DPOTRF.
41 EIGEN_LAPACK_FUNC(potrs,(char* uplo, int *n, int *nrhs, RealScalar *pa, int *lda, RealScalar *pb, int *ldb, int *info))
42 {
43  *info = 0;
44  if(UPLO(*uplo)==INVALID) *info = -1;
45  else if(*n<0) *info = -2;
46  else if(*nrhs<0) *info = -3;
47  else if(*lda<std::max(1,*n)) *info = -5;
48  else if(*ldb<std::max(1,*n)) *info = -7;
49  if(*info!=0)
50  {
51  int e = -*info;
52  return xerbla_(SCALAR_SUFFIX_UP"POTRS", &e, 6);
53  }
54 
55  Scalar* a = reinterpret_cast<Scalar*>(pa);
56  Scalar* b = reinterpret_cast<Scalar*>(pb);
57  MatrixType A(a,*n,*n,*lda);
58  MatrixType B(b,*n,*nrhs,*ldb);
59 
60  if(UPLO(*uplo)==UP)
61  {
62  A.triangularView<Upper>().adjoint().solveInPlace(B);
63  A.triangularView<Upper>().solveInPlace(B);
64  }
65  else
66  {
67  A.triangularView<Lower>().solveInPlace(B);
68  A.triangularView<Lower>().adjoint().solveInPlace(B);
69  }
70 
71  return 0;
72 }
#define SCALAR_SUFFIX_UP
SCALAR Scalar
Definition: bench_gemm.cpp:46
#define max(a, b)
Definition: datatypes.h:20
void adjoint(const MatrixType &m)
Definition: adjoint.cpp:67
int n
MatrixXf MatrixType
else if n * info
#define EIGEN_LAPACK_FUNC(FUNC, ARGLIST)
Definition: lapack_common.h:16
EIGEN_WEAK_LINKING int xerbla_(const char *msg, int *info, int)
Definition: xerbla.cpp:15
Array< double, 1, 3 > e(1./3., 0.5, 2.)
NumTraits< Scalar >::Real RealScalar
Definition: bench_gemm.cpp:47
MatrixType A(a, *n, *n, *lda)
else if * lda(1, *n)) *info=-4;if(*info!=0
MatrixType B(b, *n, *nrhs, *ldb)


gtsam
Author(s):
autogenerated on Tue Jul 4 2023 02:34:01