cholesky.cpp
Go to the documentation of this file.
00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2010-2011 Gael Guennebaud <gael.guennebaud@inria.fr>
00005 //
00006 // Eigen is free software; you can redistribute it and/or
00007 // modify it under the terms of the GNU Lesser General Public
00008 // License as published by the Free Software Foundation; either
00009 // version 3 of the License, or (at your option) any later version.
00010 //
00011 // Alternatively, you can redistribute it and/or
00012 // modify it under the terms of the GNU General Public License as
00013 // published by the Free Software Foundation; either version 2 of
00014 // the License, or (at your option) any later version.
00015 //
00016 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
00017 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00018 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
00019 // GNU General Public License for more details.
00020 //
00021 // You should have received a copy of the GNU Lesser General Public
00022 // License and a copy of the GNU General Public License along with
00023 // Eigen. If not, see <http://www.gnu.org/licenses/>.
00024 
00025 #include "lapack_common.h"
00026 #include <Eigen/Cholesky>
00027 
00028 // POTRF computes the Cholesky factorization of a real symmetric positive definite matrix A.
00029 EIGEN_LAPACK_FUNC(potrf,(char* uplo, int *n, RealScalar *pa, int *lda, int *info))
00030 {
00031   *info = 0;
00032         if(UPLO(*uplo)==INVALID) *info = -1;
00033   else  if(*n<0)                 *info = -2;
00034   else  if(*lda<std::max(1,*n))  *info = -4;
00035   if(*info!=0)
00036   {
00037     int e = -*info;
00038     return xerbla_(SCALAR_SUFFIX_UP"POTRF", &e, 6);
00039   }
00040 
00041   Scalar* a = reinterpret_cast<Scalar*>(pa);
00042   MatrixType A(a,*n,*n,*lda);
00043   int ret;
00044   if(UPLO(*uplo)==UP) ret = internal::llt_inplace<Upper>::blocked(A);
00045   else                ret = internal::llt_inplace<Lower>::blocked(A);
00046 
00047   if(ret>=0)
00048     *info = ret+1;
00049   
00050   return 0;
00051 }
00052 
00053 // POTRS solves a system of linear equations A*X = B with a symmetric
00054 // positive definite matrix A using the Cholesky factorization
00055 // A = U**T*U or A = L*L**T computed by DPOTRF.
00056 EIGEN_LAPACK_FUNC(potrs,(char* uplo, int *n, int *nrhs, RealScalar *pa, int *lda, RealScalar *pb, int *ldb, int *info))
00057 {
00058   *info = 0;
00059         if(UPLO(*uplo)==INVALID) *info = -1;
00060   else  if(*n<0)                 *info = -2;
00061   else  if(*nrhs<0)              *info = -3;
00062   else  if(*lda<std::max(1,*n))  *info = -5;
00063   else  if(*ldb<std::max(1,*n))  *info = -7;
00064   if(*info!=0)
00065   {
00066     int e = -*info;
00067     return xerbla_(SCALAR_SUFFIX_UP"POTRS", &e, 6);
00068   }
00069 
00070   Scalar* a = reinterpret_cast<Scalar*>(pa);
00071   Scalar* b = reinterpret_cast<Scalar*>(pb);
00072   MatrixType A(a,*n,*n,*lda);
00073   MatrixType B(b,*n,*nrhs,*ldb);
00074 
00075   if(UPLO(*uplo)==UP)
00076   {
00077     A.triangularView<Upper>().adjoint().solveInPlace(B);
00078     A.triangularView<Upper>().solveInPlace(B);
00079   }
00080   else
00081   {
00082     A.triangularView<Lower>().solveInPlace(B);
00083     A.triangularView<Lower>().adjoint().solveInPlace(B);
00084   }
00085 
00086   return 0;
00087 }


libicr
Author(s): Robert Krug
autogenerated on Mon Jan 6 2014 11:32:31