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. 14EIGEN_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. 41EIGEN_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