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 "common.h" 11 #include <Eigen/Eigenvalues> 12 13 // computes an LU factorization of a general M-by-N matrix A using partial pivoting with row interchanges 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 nb = ILAENV( 1, 'SSYTRD', UPLO, N, -1, -1, -1 ) 29 // LWKOPT = MAX( 1, ( NB+2 )*N ) 30 // WORK( 1 ) = LWKOPT 31 // * 32 // IF( LWORK.LT.MAX( 1, 3*N-1 ) .AND. .NOT.LQUERY ) 33 // $ INFO = -8 34 // END IF 35 // * 36 // IF( INFO.NE.0 ) THEN 37 // CALL XERBLA( 'SSYEV ', -INFO ) 38 // RETURN 39 // ELSE IF( LQUERY ) THEN 40 // RETURN 41 // END IF 42 43 if(*info!=0) 44 { 45 int e = -*info; 46 return xerbla_(SCALAR_SUFFIX_UP"SYEV ", &e, 6); 47 } 48 49 if(query_size) 50 { 51 *lwork = 0; 52 return 0; 53 } 54 55 if(*n==0) 56 return 0; 57 58 PlainMatrixType mat(*n,*n); 59 if(UPLO(*uplo)==UP) mat = matrix(a,*n,*n,*lda).adjoint(); 60 else mat = matrix(a,*n,*n,*lda); 61 62 bool computeVectors = *jobz=='V' || *jobz=='v'; 63 SelfAdjointEigenSolver<PlainMatrixType> eig(mat,computeVectors?ComputeEigenvectors:EigenvaluesOnly); 64 65 if(eig.info()==NoConvergence) 66 { 67 vector(w,*n).setZero(); 68 if(computeVectors) 69 matrix(a,*n,*n,*lda).setIdentity(); 70 //*info = 1; 71 return 0; 72 } 73 74 vector(w,*n) = eig.eigenvalues(); 75 if(computeVectors) 76 matrix(a,*n,*n,*lda) = eig.eigenvectors(); 77 78 return 0; 79 } 80
