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 "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 = internal::llt_inplace<Scalar, Upper>::blocked(A); 30 else ret = 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 } 73
