1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 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/Eigenvalues> 12 13 // computes eigen values and vectors of a general N-by-N matrix A 14 EIGEN_LAPACK_FUNC(syev,(char *jobz, char *uplo, int* n, Scalar* a, int *lda, Scalar* w, Scalar* /*work*/, int* lwork, int *info)) 15 { 16 // TODO exploit the work buffer 17 bool query_size = *lwork==-1; 18 19 *info = 0; 20 if(*jobz!='N' && *jobz!='V') *info = -1; 21 else if(UPLO(*uplo)==INVALID) *info = -2; 22 else if(*n<0) *info = -3; 23 else if(*lda<std::max(1,*n)) *info = -5; 24 else if((!query_size) && *lwork<std::max(1,3**n-1)) *info = -8; 25 26 if(*info!=0) 27 { 28 int e = -*info; 29 return xerbla_(SCALAR_SUFFIX_UP"SYEV ", &e, 6); 30 } 31 32 if(query_size) 33 { 34 *lwork = 0; 35 return 0; 36 } 37 38 if(*n==0) 39 return 0; 40 41 PlainMatrixType mat(*n,*n); 42 if(UPLO(*uplo)==UP) mat = matrix(a,*n,*n,*lda).adjoint(); 43 else mat = matrix(a,*n,*n,*lda); 44 45 bool computeVectors = *jobz=='V' || *jobz=='v'; 46 SelfAdjointEigenSolver<PlainMatrixType> eig(mat,computeVectors?ComputeEigenvectors:EigenvaluesOnly); 47 48 if(eig.info()==NoConvergence) 49 { 50 make_vector(w,*n).setZero(); 51 if(computeVectors) 52 matrix(a,*n,*n,*lda).setIdentity(); 53 //*info = 1; 54 return 0; 55 } 56 57 make_vector(w,*n) = eig.eigenvalues(); 58 if(computeVectors) 59 matrix(a,*n,*n,*lda) = eig.eigenvectors(); 60 61 return 0; 62 } 63
