fortran用MKL库求对称矩阵本征问题

最新推荐文章于 2024-05-10 23:12:15 发布

置顶 xhh22900

最新推荐文章于 2024-05-10 23:12:15 发布

阅读量809

点赞数 1

分类专栏： Fortran C/C++ 笔记文章标签：算法 c++ 线性代数

本文链接：https://blog.csdn.net/xhh22900/article/details/117304235

版权

笔记同时被 3 个专栏收录

35 篇文章 3 订阅

订阅专栏

C/C++

31 篇文章 3 订阅

订阅专栏

Fortran

3 篇文章 3 订阅

订阅专栏

安装oneAPI（Intel Fortran编译器）

Intel oneAPI Base Toolkit
Intel oneAPI HPC Toolkit
可以选择下载安装或者在线安装，下载的话需要用教育邮箱登陆。

尝试MKL库

参考官方文档Symmetric Eigenproblems。
比如：

Fortran

program main
    implicit none
    integer*4,parameter::n=5,lda=5,lwmax=1000
!     .. Local Scalars ..
    INTEGER::INFO, LWORK
!
!     .. Local Arrays ..
    REAL*8::A( LDA, N ),  W( N ), WORK( LWMAX )
    DATA             A/&
       &1.96, 0.00, 0.00, 0.00, 0.00,&
      &-6.49, 3.80, 0.00, 0.00, 0.00,&
      &-0.47,-6.39, 4.17, 0.00, 0.00,&
      &-7.20, 1.50,-1.51, 5.70, 0.00,&
      &-0.65,-6.34, 2.67, 1.80,-7.10/
!
!     .. External Subroutines ..
    EXTERNAL         DSYEV
    EXTERNAL         PRINT_MATRIX
!
!     .. Intrinsic Functions ..
    INTRINSIC        INT, MIN
!
!     .. Executable Statements ..
    WRITE(*,*)'SSYEV Example Program Results'
!     
!     Query the optimal workspace.
!       
    LWORK = -1
    CALL DSYEV( 'Vectors', 'Upper', N, A, LDA, W, WORK, LWORK, INFO )
    LWORK = MIN( LWMAX, INT( WORK( 1 ) ) )
!
!     Solve eigenproblem.
!
    CALL DSYEV( 'Vectors', 'Upper', N, A, LDA, W, WORK, LWORK, INFO )
!
!     Check for convergence.
!
    IF( INFO.GT.0 ) THEN
        WRITE(*,*)'The algorithm failed to compute eigenvalues.'
        STOP
    END IF
!
!     Print eigenvalues.

    !call insert_sort(w,N)
    CALL PRINT_MATRIX( 'Eigenvalues', 1, N, W, 1 )
!      
!     Print eigenvectors.
!
    CALL PRINT_MATRIX( 'Eigenvectors (stored columnwise)', N, N, A,LDA )
    STOP
end program
!
!     End of SSYEV Example.
!
!  =============================================================================
!
!     Auxiliary routine: printing a matrix.
!
subroutine PRINT_MATRIX( DESC, M, N, A, LDA )
    implicit none
    CHARACTER*(*)    DESC
    INTEGER          M, N, LDA
    REAL*8             A( LDA, * )
    !
    INTEGER          I, J
    !
    WRITE(*,*)
    WRITE(*,*) DESC
    DO I = 1, M
        WRITE(*,9998) ( A( I, J ), J = 1, N )
    END DO
    !
    9998 FORMAT( 11(:,1X,e30.20) )
    RETURN
end subroutine

注意编译时后面要加上-mkl，如

ifort test.f90 -o test -mkl

最后的结果为

SSYEV Example Program Results
 
 Eigenvalues
    -0.11065575232626278179E+02    -0.62287466937218827212E+01     0.86402803023585872388E+00     0.88654570265779408800E+01     0.16094836840924127586E+02

 Eigenvectors (stored columnwise)
    -0.29806697142941551704E+00    -0.60751344955327046815E+00     0.40261995251418697395E+00    -0.37448098436165483394E+00     0.48963726915317773436E+00
    -0.50779843411371539119E+00    -0.28796757291700009196E+00    -0.40658568098689418235E+00    -0.35716882040996411618E+00    -0.60525527259440614625E+00
    -0.81606186905410313392E-01    -0.38432041912929776339E+00    -0.65996550612182380835E+00     0.50076383196156670774E+00     0.39914829453170985740E+00
    -0.35892968587850466143E-02    -0.44672977485209214299E+00     0.45528986722532166498E+00     0.62036521460984062060E+00    -0.45637458574365619146E+00
    -0.80412957790481354170E+00     0.44803170564264482856E+00     0.17245847437922906531E+00     0.31076842049946928892E+00     0.16224757660015381999E+00

由于输出的本征值并不总是按从小到大的顺序的，所以有时我们需要一个子程序对本征值数组排序，排序的算法很多，比较简单的如插入排序：

subroutine insert_sort(A,num)
   implicit none
   real*8::A(*),key
   integer*4::num,i,j,k,temp
   !num=size(A)
   do j=1,num
       key=A(j)
       i=j-1
       temp=j
       do k=i,1,-1
           if(A(k)<key)then
               A(k+1)=A(k)
               temp=k
           endif
       enddo
       A(temp)=key
   enddo
end subroutine

C++

类似地

#include <stdio.h>
#include <stdlib.h>
#include <mkl.h>
//extern void dsyev( char* jobz, char* uplo, int* n, double* a, int* lda,
           //     double* w, double* work, int* lwork, int* info );
/* Auxiliary routines prototypes */
#define N 5 
#define LDA N
void print_matrix( char* desc, int m, int n, double* a, int lda );

int main(){
	int n = N, lda = LDA, info, lwork;
	double wkopt;
	double* work;
	/* Local arrays */
	double w[N];
	double a[LDA*N] = {
		1.96,  0.00,  0.00,  0.00,  0.00,
		-6.49,  3.80,  0.00,  0.00,  0.00,
		-0.47, -6.39,  4.17,  0.00,  0.00,
		-7.20,  1.50, -1.51,  5.70,  0.00,
		-0.65, -6.34,  2.67,  1.80, -7.10
	};
	printf( " DSYEV Example Program Results\n" );
    /* Query and allocate the optimal workspace */
    lwork = -1;
    dsyev( "Vectors", "Upper", &n, a, &lda, w, &wkopt, &lwork, &info );
    lwork = (int)wkopt;
    work = (double*)malloc( lwork*sizeof(double) );
    /* Solve eigenproblem */
    dsyev( "Vectors", "Upper", &n, a, &lda, w, work, &lwork, &info );
    /* Check for convergence */
    if( info > 0 ) {
        printf( "The algorithm failed to compute eigenvalues.\n" );
        exit( 1 );
    }
    /* Print eigenvalues */
    print_matrix( "Eigenvalues", 1, n, w, 1 );
    /* Print eigenvectors */
    print_matrix( "Eigenvectors (stored columnwise)", n, n, a, lda );
    /* Free workspace */
    free( (void*)work );
    exit( 0 );
}
void print_matrix( char* desc, int m, int n, double* a, int lda ) {
        int i, j;
        printf( "\n %s\n", desc );
        for( i = 0; i < m; i++ ) {
                for( j = 0; j < n; j++ ) printf( " %6.2f", a[i+j*lda] );
                printf( "\n" );
        }
}

编译时也需要加上-mkl

icc demo1.cpp -o demo1 -mkl

结果为

 DSYEV Example Program Results

 Eigenvalues
 -11.07  -6.23   0.86   8.87  16.09

 Eigenvectors (stored columnwise)
  -0.30  -0.61   0.40  -0.37   0.49
  -0.51  -0.29  -0.41  -0.36  -0.61
  -0.08  -0.38  -0.66   0.50   0.40
  -0.00  -0.45   0.46   0.62  -0.46
  -0.80   0.45   0.17   0.31   0.16