Hi, thanks for your interest in the new Extended Eigensolver Package! Let me do some research and I will get back to you in a day with some comprehensive information to help you out.

Dear Nikita,

Let us assume we want to solve a real generalized eigenvalue problem.

In the case of ijob=10 we need to compute  a double complex precision matrix  Ze*B-A and the representation  of the output matrix Ze*B-A  depends on the solver you are going to use.  For example, if we going to use the MKL PARDISO, we have to compute Ze*B-A in the compressed sparse row format. If you are using a sparse direct solver,  we need to factorize  Ze*B-A with the help of this direct solver.  In the case of a preconditioned iterative solver, you might compute a preconditioner for the linear system (Ze*A-B) X= Y.    Of course, you can compute a preconditioner or performs the matrix factorization later when FEAST ask you to solve the system.

Here is a pseudocode for a real symmetric problem:

ijob=-1 ! initialization

do while (ijob/=0) 

    call dfeast_srci(ijob, N, Ze, work, workc, Aq, Sq, fpm,  epsout,loop, Emin, Emax, M0,  E, X, mode,res,info)    

    select case(ijob)

    case(10) !! Form (ZeB-A)  and Factorize the complex matrix (ZeB-A)

                   !!  Let (ssa,sisa,sjsa) be the CSR representation  of the input matrix A, 

                   !!   (ssb,sisb,sjsb) be the CSR representation  of the input matrix B

                   !!   and let (saz,isaz,jsaz) be the  CSR representation of (Ze*A-B) 

                   call <compute the resolvent matrix > (N, ssa,sisa,sjsa,Ze,ssb,sisb,sjsb,saz,isaz,jsaz) !! get saz

                  !!  factorize the matrix Ze*B-A or

                  !!   compute a preconditioner for (Ze*A-B) if a preconditioned iterative solver is used

                  PHASE=22

                  call  pardiso(….,saz,isaz,jsaz,….)         

    case(11) !! Solve the complex linear system (ZeB-A)workc(1:N,1:M0) = workc(1:N,1:M0) result in workc

                  PHASE=33

                  call  pardiso(…, saz,isaz,jsaz, … ,workc(:,1:M0), …)         

     case(30) !! Perform multiplication A*X(1:N,i:j) result in work(1:N,i:j)

                     !! where i=fpm(24) and j=fpm(24)+fpm(25)−1

                    call  <sparse matrix-matrix multiply>   (…, ssa, sjsa, sisa, …  , X(1,fpm(24)), ..,  work(1,fpm(24)))

     case(40) !! Perform multiplication B*X(1:N,i:j) result in work(1:N,i:j)

                     !! where i=fpm(24) and j=fpm(24)+fpm(25)−1

                    call  <sparse matrix-matrix multiply>   (…, ssb, sjsb, sisb, …  , X(1,fpm(24)), ..,  work(1,fpm(24))) 

      end select

end do

Hope it helps

All the best

Sergey

Hi Nikita,

I'd assume that  it is not enough memory for one or more input arrays.  I'd recommend to check the type and size of each input array. If it doesn't help,  please provide a reproducer.

All the best

Sergey

Dear Nikita,

I investigated the issue and the problem appears because of a bug in  RCI  Extended Eigensolver  interfaces from the MKL LP64 binaries. Everything should work well if you link with MKL ILP64 binaries like this

icc -I$MKL_ROOT/__release_lnx/mkl/include -DMKL_INT="long long" code.cpp $MKL_ROOT/__release_lnx/mkl/lib/intel64/libmkl_intel_ilp64.a -Wl,--start-group $MKL_ROOT/__release_lnx/mkl/lib/intel64/lib/libmkl_intel_thread.a $MKL_ROOT/__release_lnx/mkl/lib/intel64/libmkl_core.a -Wl,--end-group   $MKL_ROOT/__release_lnx/compiler/lib/intel64/libiomp5.so -lpthread

Besides  ILP64 RCI interfaces work a bit faster since the computational kernels are ILP64 and  additional wrappers are not needed.

Unfortunately the Intel MKL doesn’t contain any iterative solver for complex linear systems. Feel free to submit a feature request. But you can try the MKL Out-of-Core PARDISO  with MKL Extended Eigensolver in the case of large sparse matrices.

Thanks for catching this issue.

All the best

Sergey

the issue is escalated - here is the issue number for you reference : DPD200329829

1. we will let you know here into this thread when the fix of the issue would be available.

2 - you can alos will find this info into Release Note.