1c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// This file is part of Eigen, a lightweight C++ template library 2c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// for linear algebra. 3c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// 4c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Copyright (C) 2011 Gael Guennebaud <gael.guennebaud@inria.fr> 5c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// 6c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// This Source Code Form is subject to the terms of the Mozilla 7c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Public License v. 2.0. If a copy of the MPL was not distributed 8c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 9c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 10c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#include "common.h" 11c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#include <Eigen/Eigenvalues> 12c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 13c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// computes an LU factorization of a general M-by-N matrix A using partial pivoting with row interchanges 14c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan KamathEIGEN_LAPACK_FUNC(syev,(char *jobz, char *uplo, int* n, Scalar* a, int *lda, Scalar* w, Scalar* /*work*/, int* lwork, int *info)) 15c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 16c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // TODO exploit the work buffer 17c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bool query_size = *lwork==-1; 18c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 19c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath *info = 0; 20c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(*jobz!='N' && *jobz!='V') *info = -1; 21c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else if(UPLO(*uplo)==INVALID) *info = -2; 22c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else if(*n<0) *info = -3; 23c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else if(*lda<std::max(1,*n)) *info = -5; 24c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else if((!query_size) && *lwork<std::max(1,3**n-1)) *info = -8; 25c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 26c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// if(*info==0) 27c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// { 28c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// int nb = ILAENV( 1, 'SSYTRD', UPLO, N, -1, -1, -1 ) 29c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// LWKOPT = MAX( 1, ( NB+2 )*N ) 30c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// WORK( 1 ) = LWKOPT 31c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// * 32c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// IF( LWORK.LT.MAX( 1, 3*N-1 ) .AND. .NOT.LQUERY ) 33c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// $ INFO = -8 34c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// END IF 35c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// * 36c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// IF( INFO.NE.0 ) THEN 37c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// CALL XERBLA( 'SSYEV ', -INFO ) 38c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// RETURN 39c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// ELSE IF( LQUERY ) THEN 40c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// RETURN 41c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// END IF 42c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 43c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(*info!=0) 44c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 45c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int e = -*info; 46c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return xerbla_(SCALAR_SUFFIX_UP"SYEV ", &e, 6); 47c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 48c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 49c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(query_size) 50c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 51c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath *lwork = 0; 52c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 0; 53c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 54c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 55c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(*n==0) 56c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 0; 57c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 58c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath PlainMatrixType mat(*n,*n); 59c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(UPLO(*uplo)==UP) mat = matrix(a,*n,*n,*lda).adjoint(); 60c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else mat = matrix(a,*n,*n,*lda); 61c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 62c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bool computeVectors = *jobz=='V' || *jobz=='v'; 63c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath SelfAdjointEigenSolver<PlainMatrixType> eig(mat,computeVectors?ComputeEigenvectors:EigenvaluesOnly); 64c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 65c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(eig.info()==NoConvergence) 66c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 67c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath vector(w,*n).setZero(); 68c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(computeVectors) 69c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath matrix(a,*n,*n,*lda).setIdentity(); 70c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath //*info = 1; 71c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 0; 72c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 73c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 74c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath vector(w,*n) = eig.eigenvalues(); 75c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(computeVectors) 76c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath matrix(a,*n,*n,*lda) = eig.eigenvectors(); 77c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 78c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 0; 79c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 80