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 }
gtsam.examples.DogLegOptimizerExample.int
int
Definition: DogLegOptimizerExample.py:111
UPLO
#define UPLO(X)
Definition: gtsam/3rdparty/Eigen/blas/common.h:56
e
Array< double, 1, 3 > e(1./3., 0.5, 2.)
ret
int ret
Definition: 3rdparty/Eigen/lapack/cholesky.cpp:28
MatrixType
MatrixXf MatrixType
Definition: benchmark-blocking-sizes.cpp:52
B
MatrixType B(b, *n, *nrhs, *ldb)
EIGEN_LAPACK_FUNC
#define EIGEN_LAPACK_FUNC(FUNC, ARGLIST)
Definition: lapack_common.h:16
Eigen::Upper
@ Upper
Definition: Constants.h:211
a
Scalar * a
Definition: 3rdparty/Eigen/lapack/cholesky.cpp:26
b
Scalar * b
Definition: 3rdparty/Eigen/lapack/cholesky.cpp:56
n
int n
Definition: BiCGSTAB_simple.cpp:1
A
MatrixType A(a, *n, *n, *lda)
info
else if n * info
Definition: 3rdparty/Eigen/lapack/cholesky.cpp:18
SCALAR_SUFFIX_UP
#define SCALAR_SUFFIX_UP
Definition: blas/complex_double.cpp:12
lapack_common.h
Eigen::Lower
@ Lower
Definition: Constants.h:209
RealScalar
NumTraits< Scalar >::Real RealScalar
Definition: bench_gemm.cpp:47
xerbla_
EIGEN_WEAK_LINKING int xerbla_(const char *msg, int *info, int)
Definition: xerbla.cpp:15
lda
else if * lda(1, *n)) *info=-4;if(*info!=0
Definition: 3rdparty/Eigen/lapack/cholesky.cpp:19
max
#define max(a, b)
Definition: datatypes.h:20
INVALID
#define INVALID
Definition: gtsam/3rdparty/Eigen/blas/common.h:45
adjoint
void adjoint(const MatrixType &m)
Definition: adjoint.cpp:67
UP
#define UP
Definition: gtsam/3rdparty/Eigen/blas/common.h:39
Scalar
SCALAR Scalar
Definition: bench_gemm.cpp:46


gtsam
Author(s):
autogenerated on Fri Nov 1 2024 03:32:05