2011-07-06 274 views
1

我想計算不同大小矩陣的特徵值和特徵向量。我正在使用一段非常簡單的Fortran90代碼,並且正在編譯鏈接到我的機器中提供的適用於英特爾MKL軟件包的Lapack庫,該機器在Ubuntu中運行。代碼「matrix_diag_01.f90」附在消息的末尾。 「隨機」模塊只包含Numerical Recipes中的「跑」隨機數發生器。代碼編譯得很好使用在Fortran中使用mkl lapack庫的特徵值和特徵向量

ifort -I $(MKLPATH) -o matrix_diag_01 matrix_diag_01.f90 
     random.f90 $(MKLPATH)libmkl_lapack95.a -Wl,--start-group 
     $(MKLPATH)libmkl_intel_lp64.a $(MKLPATH)libmkl_lapack.a 
     $(MKLPATH)libmkl_intel_thread.a $(MKLPATH)libmkl_core.a 
     -Wl,--end-group -lguide -lpthread 

當給出小的矩陣時,可執行文件很好地工作。但是,對於大小爲3000x3000的矩陣,它會產生一些奇怪的行爲。首先它出現這個錯誤

MKL ERROR : Parameter 8 was incorrect on entry to SSYEVD 

但是,在調用SSYEVD時只有3個參數。其次,它返回特徵向量,但不返回特徵值。我通過在另一臺具有更大內存的機器上編譯進行檢查,但結果是一樣的。

任何人都可以請幫忙嗎?

謝謝!

PROGRAM matrix_diag_01 

USE random 

IMPLICIT NONE 

INTERFACE 
    SUBROUTINE diag(mat,n) 
     INTEGER n 
     REAL,DIMENSION(n,n) :: mat 
    END SUBROUTINE 
END INTERFACE 

INTEGER n,i,j,iseed 
REAL, DIMENSION(:), ALLOCATABLE :: w 
REAL, DIMENSION(:,:), ALLOCATABLE :: mat 

write (*,*) ' Please enter size of matrix' 
read (*,*) n 
write (*,*) ' Please type seed' 
read (*,*) iseed 

allocate (mat(n,n)) 

do i = 1,n 
    do j = 1,n 
     mat(i,j) = ran(iseed) 
    end do 
end do 

call diag(mat,n) 

stop 
END 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
SUBROUTINE diag(mat,n) 

    USE mkl95_lapack 
    USE mkl95_precision 

    IMPLICIT NONE 

    CHARACTER(len=1) :: jobz = 'V' 

    INTEGER n,i 

    REAL,DIMENSION(n,n) :: mat 
    REAL,DIMENSION(:,:),ALLOCATABLE :: matt,a 
    REAL,DIMENSION(:),ALLOCATABLE :: w 

    allocate (matt(n,n),a(n,n),w(n)) 

    matt = mat*transpose(mat) 
    a = sqrt(matt) 
    open (unit=7,file="matrix.dat",status="unknown") 
    do i = 1,n 
    write (7,100) a(i,:) 
    end do 
    close (unit=7) 

    call syevd(a,w,jobz) 

    open (unit=8,file="eig_val.dat",status="unknown") 
    do i = 1,n 
    write (8,100) w(i) 
    end do 
    close (unit=8) 

    open (unit=9,file="eig_vec.dat",status="unknown") 
    do i = 1,n 
    write (9,100) a(i,:) 
    end do 
    close (unit=9) 

    return 
100 format(5000f16.5) 
end 
+0

我沒有解決MKL錯誤,雖然沒有足夠的可用內存聽起來不可能。關於不返回的特徵值,你有沒有試過用'info'參數調用'ssyevd'? – eriktous

回答

2

如果你看一下您所呼叫的函數的源代碼,它本身發出以下電話: http://software.intel.com/sites/products/documentation/hpc/mkl/lapack/mkl_lapack_examples/ssyevd_ex.f.htm

LWORK = MIN(LWMAX, INT(WORK(1))) 

CALL SSYEVD('Vectors', 'Upper', N, A, LDA, W, WORK, LWORK, 
     $    IWORK, LIWORK, INFO) 

想必參數LWORK是錯誤的。

+0

非常感謝您的回答,whoplisp,但是我調用函數和鏈接庫(即Fortran 95)的方式不能直接控制LWORK參數。請參閱http://software.intel.com/sites/products/documentation/hpc/mkl/mklman/lse/functn_syevd.htm – ricoamor

+0

是的,因爲您正在調用包裝。也許嘗試調用Fortran77版本。另請參閱有關限制的功能文檔。完整的3000x3000矩陣似乎相當大。即使使用double,也許會出現穩定或舍入錯誤的問題。 – whoplisp

+0

非常感謝。這似乎是未來的方向。只需使用Fortran 77代碼的這一部分,並忘記包裝。它適用於3000x3000矩陣。謝謝! – ricoamor

相關問題