TABLE OF CONTENTS
- ABINIT/jacobi
- ABINIT/m_hide_lapack
- m_cgtools/xhesv_cplex
- m_hide_lapack/cginv
- m_hide_lapack/dzegdi
- m_hide_lapack/dzgefa
- m_hide_lapack/lubksb
- m_hide_lapack/ludcmp
- m_hide_lapack/matr3eigval
- m_hide_lapack/matrginv
- m_hide_lapack/test_xginv
- m_hide_lapack/wrap_CGEEV
- m_hide_lapack/wrap_CHEEV
- m_hide_lapack/wrap_CHPEV
- m_hide_lapack/wrap_ZGEEV
- m_hide_lapack/wrap_ZHEEV
- m_hide_lapack/wrap_ZHEEVX
- m_hide_lapack/wrap_ZHEGV
- m_hide_lapack/wrap_ZHEGVX
- m_hide_lapack/wrap_ZHPEV
- m_hide_lapack/xheev_cplex
- m_hide_lapack/xheevx_cplex
- m_hide_lapack/xhegv_cplex
- m_hide_lapack/xhegvx_cplex
- m_hide_lapack/zginv
- m_hide_lapack/zhpd_invert
ABINIT/jacobi [ Functions ]
NAME
jacobi
FUNCTION
Computes all eigenvalues and eigenvectors of a real symmetric matrix a, which is of size n by n, stored in a physical np by np array. On output, elements of a above the diagonal are destroyed. d returns the eigenvalues of a in its first n elements. v is a matrix with the same logical and physical dimensions as a, whose columns contain, on output, the normalized eigenvectors of a. nrot returns the number of Jacobi rotations that were required.
INPUTS
OUTPUT
NOTES
This routine is deprecated, use Lapack API
SOURCE
3131 subroutine jacobi(a,n,np,d,v,nrot) 3132 3133 !Arguments 3134 integer :: n,np,nrot 3135 real*8 :: a(np,np),d(np),v(np,np) 3136 !Local variables 3137 integer, parameter :: NMAX=500 3138 integer i,ip,iq,j 3139 real*8 c,g,h,s,sm,t,tau,theta,tresh,b(NMAX),z(NMAX) 3140 do ip=1,n 3141 do iq=1,n 3142 v(ip,iq)=0. 3143 enddo 3144 v(ip,ip)=1. 3145 enddo 3146 do ip=1,n 3147 b(ip)=a(ip,ip) 3148 d(ip)=b(ip) 3149 z(ip)=0. 3150 enddo 3151 nrot=0 3152 do i=1,50 3153 sm=0. 3154 do ip=1,n-1 3155 do iq=ip+1,n 3156 sm=sm+abs(a(ip,iq)) 3157 enddo 3158 enddo 3159 if(sm.eq.0.)return 3160 if(i.lt.4)then 3161 tresh=0.2*sm/n**2 3162 else 3163 tresh=0. 3164 endif 3165 do ip=1,n-1 3166 do iq=ip+1,n 3167 g=100.*abs(a(ip,iq)) 3168 if((i.gt.4).and.(abs(d(ip))+g.eq.abs(d(ip))) & 3169 & .and.(abs(d(iq))+g.eq.abs(d(iq))))then 3170 a(ip,iq)=0. 3171 else if(abs(a(ip,iq)).gt.tresh)then 3172 h=d(iq)-d(ip) 3173 if(abs(h)+g.eq.abs(h))then 3174 t=a(ip,iq)/h 3175 else 3176 theta=0.5*h/a(ip,iq) 3177 t=1./(abs(theta)+sqrt(1.+theta**2)) 3178 if(theta.lt.0.)t=-t 3179 endif 3180 c=1./sqrt(1+t**2) 3181 s=t*c 3182 tau=s/(1.+c) 3183 h=t*a(ip,iq) 3184 z(ip)=z(ip)-h 3185 z(iq)=z(iq)+h 3186 d(ip)=d(ip)-h 3187 d(iq)=d(iq)+h 3188 a(ip,iq)=0. 3189 do j=1,ip-1 3190 g=a(j,ip) 3191 h=a(j,iq) 3192 a(j,ip)=g-s*(h+g*tau) 3193 a(j,iq)=h+s*(g-h*tau) 3194 enddo 3195 do j=ip+1,iq-1 3196 g=a(ip,j) 3197 h=a(j,iq) 3198 a(ip,j)=g-s*(h+g*tau) 3199 a(j,iq)=h+s*(g-h*tau) 3200 enddo 3201 do j=iq+1,n 3202 g=a(ip,j) 3203 h=a(iq,j) 3204 a(ip,j)=g-s*(h+g*tau) 3205 a(iq,j)=h+s*(g-h*tau) 3206 enddo 3207 do j=1,n 3208 g=v(j,ip) 3209 h=v(j,iq) 3210 v(j,ip)=g-s*(h+g*tau) 3211 v(j,iq)=h+s*(g-h*tau) 3212 enddo 3213 nrot=nrot+1 3214 endif 3215 enddo 3216 enddo 3217 do ip=1,n 3218 b(ip)=b(ip)+z(ip) 3219 d(ip)=b(ip) 3220 z(ip)=0. 3221 enddo 3222 enddo 3223 write(std_out,*) 'too many iterations in jacobi' 3224 3225 end subroutine jacobi
ABINIT/m_hide_lapack [ Modules ]
NAME
m_hide_lapack
FUNCTION
ABINIT Linear Algebra Subroutine Interfaces. This modules provides interfaces performing the overloading of commonly used Lapack routines. The main purpose of this module is to create a layer between abinit routines and Lapack procedures. This layer can be used to hide the parallel Scalapack version. In this case, only the MPI commutator has to be provided in input as the wrapper will take care of the initialization of the Scalapack grid as well as of the distribution of the matrix. Note that this allows one to reduce the CPU time per processor but not the memory allocated since the entire matrix has to be provided in input. The interfaces are very similar to the Lapack F77 version (neither F90 constructs nor F90 assumed size arrays are used). The main simplification with respect to the F77 version of Lapack is that the work arrays are allocated inside the wrapper with optimal size thus reducing the number of input argcomm_scalapackuments that has to be passed. Leading dimensions have been removed from the interface whenever possible. In F90 one can pass the array descriptor if the routines should operate on a slice of the local array (seldom done in abinit). Using array descriptor is OK but will it likely slow-down the calculation as some compilers perform a copy of the input-output data. If efficiency is a concern, then the F77 call should be used
COPYRIGHT
Copyright (C) 1992-2024 ABINIT group (MG, GMR, XG) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
TODO
1) Use a function to define the size of the Scalapack block according to some heuristic method. 2) Define a threshold below which Scalapack is not used although the MPI communicator is passed. 3) Split MPI communicator for Scalapack (.i.e. use a max size for the Scalapack comm; see abi_linalg_init). 4) On certain networks, xmpi_sum might crash due to the size of the MPI packet. This problem should be solved in hide_mpi (Module containing a private global variable defining a threshold above which the input array is split into smaller chunks.
SOURCE
42 #if defined HAVE_CONFIG_H 43 #include "config.h" 44 #endif 45 46 #include "abi_common.h" 47 48 MODULE m_hide_lapack 49 50 use defs_basis 51 use m_abicore 52 use m_xmpi 53 use m_errors 54 use m_slk 55 56 use m_linalg_interfaces 57 58 use m_time, only : cwtime 59 use m_fstrings, only : firstchar 60 !use m_slk, only : matrix_scalapack, processor_scalapack 61 62 implicit none 63 64 private 65 66 ! ========================== 67 ! Complex Hermitian matrices 68 ! ========================== 69 ! The _cplex version receives matrices declared as arr(cplex,N,M) 70 ! and calls the complex/real version depending on cplex (used e.g. for istwf_k = 2) 71 72 public :: xheev ! Computes all the eigenvalues and, optionally, eigenvectors of a complex Hermitian matrix. 73 public :: xheev_cplex 74 75 public :: xhpev ! Computes all the eigenvalues and, optionally, eigenvectors of a complex Hermitian matrix 76 ! in packed storage (Scalapack version not available) 77 78 public :: xhegv ! Compute all the eigenvalues, and optionally, the eigenvectors of a complex generalized 79 ! Hermitian-definite eigenproblem, of the form: 80 ! A*x=(lambda)*B*x, A*Bx=(lambda)*x, or B*A*x=(lambda)*x 81 82 public :: xheevx ! Computes selected eigenvalues and, optionally, eigenvectors of a complex Hermitian matrix A. 83 ! Eigenvalues and eigenvectors can be selected by specifying either a range of values or a range of 84 ! indices for the desired eigenvalues. 85 86 public :: xheevx_cplex 87 88 public :: xhegvx ! Computes selected eigenvalues, and optionally, eigenvectors 89 ! of a complex generalized Hermitian-definite eigenproblem, 90 ! of the form A*x=(lambda)*B*x, A*Bx=(lambda)*x, or B*A*x=(lambda)*x. 91 ! Eigenvalues and eigenvectors can be selected by specifying either a range of values or a range of 92 ! indices for the desired eigenvalues. 93 public :: xhegvx_cplex 94 95 96 public :: xhesv_cplex ! Solve A * X = B, where A is an N-by-N Hermitian matrix and X and B are N-by-NRHS matrices. 97 98 99 ! ============================== 100 ! Complex non-symmetric matrices 101 ! ============================== 102 103 public :: xgeev ! Computes for a complex nonsymmetric matrix A, the eigenvalues and, optionally, 104 ! the left and/or right eigenvectors. 105 106 public :: xginv ! Invert a general matrix of complex elements by means of LU factorization. 107 108 109 public :: xhdp_invert ! Invert a Hermitian positive definite matrix. 110 111 interface xheev 112 module procedure wrap_CHEEV 113 module procedure wrap_ZHEEV 114 end interface xheev 115 116 interface xhpev 117 module procedure wrap_CHPEV 118 module procedure wrap_ZHPEV 119 end interface xhpev 120 121 interface xhegv 122 module procedure wrap_ZHEGV 123 end interface xhegv 124 125 public :: xhegv_cplex 126 127 interface xheevx 128 module procedure wrap_ZHEEVX 129 end interface xheevx 130 131 interface xhegvx 132 module procedure wrap_ZHEGVX 133 end interface xhegvx 134 135 interface xgeev 136 module procedure wrap_CGEEV 137 module procedure wrap_ZGEEV 138 end interface xgeev 139 140 interface xginv 141 module procedure cginv 142 module procedure zginv 143 end interface xginv 144 145 interface xhdp_invert 146 module procedure zhpd_invert 147 end interface xhdp_invert 148 149 public :: matrginv ! Invert a general matrix of real*8 elements. 150 public :: matr3eigval ! Find the eigenvalues of a real symmetric 3x3 matrix, entered in full storage mode. 151 152 !FIXME This procedures are deprecated, use lapack API 153 public :: jacobi ! Computes all eigenvalues and eigenvectors of a real symmetric matrix a, 154 public :: ludcmp 155 public :: lubksb 156 public :: dzgedi 157 public :: dzgefa 158 159 !---------------------------------------------------------------------- 160 ! support for unitary tests and profiling. 161 162 type,public :: latime_t 163 character(len=500) :: testname 164 integer :: msize 165 real(dp) :: ctime 166 real(dp) :: wtime 167 real(dp) :: max_abserr=-one 168 real(dp) :: gflops 169 end type latime_t 170 171 public :: test_xginv 172 173 !---------------------------------------------------------------------- 174 ! private variables 175 176 integer,private,parameter :: SLK_BLOCK_SIZE = 24 177 ! Default block size for Scalapack distribution. 178 ! As recommended by Intel MKL, a more sensible default than the previous value of 40 179 180 CONTAINS !=========================================================================================================================
m_cgtools/xhesv_cplex [ Functions ]
[ Top ] [ m_cgtools ] [ Functions ]
NAME
xhesv_cplex
FUNCTION
ZHESV computes the solution to a complex (real) system of linear equations
A * X = B,
where A is an N-by-N Hermitian matrix and X and B are N-by-NRHS matrices.
The value of cplex (1 or 2) defines whether we have a complex Hermitian or real symmetric matrix
The diagonal pivoting method is used to factor A as
A = U * D * U**H, if UPLO = 'U', or
A = L * D * L**H, if UPLO = 'L',
where U (or L) is a product of permutation and unit upper (lower)
triangular matrices, and D is Hermitian and block diagonal with
1-by-1 and 2-by-2 diagonal blocks. The factored form of A is then
used to solve the system of equations A * X = B.
INPUTS
UPLO is CHARACTER*1
= 'U': Upper triangle of A is stored;
= 'L': Lower triangle of A is stored.
N is INTEGER
The number of linear equations, i.e., the order of the
matrix A. N >= 0.
NRHS is INTEGER
The number of right hand sides, i.e., the number of columns
of the matrix B. NRHS >= 0.
A is COMPLEX*16 array, dimension (LDA,N)
On entry, the Hermitian matrix A. If UPLO = 'U', the leading
N-by-N upper triangular part of A contains the upper
triangular part of the matrix A, and the strictly lower
triangular part of A is not referenced. If UPLO = 'L', the
leading N-by-N lower triangular part of A contains the lower
triangular part of the matrix A, and the strictly upper
triangular part of A is not referenced.
On exit, if INFO = 0, the block diagonal matrix D and the
multipliers used to obtain the factor U or L from the
factorization A = U*D*U**H or A = L*D*L**H as computed by
ZHETRF.
B is COMPLEX*16 array, dimension (LDB,NRHS)
On entry, the N-by-NRHS right hand side matrix B.
On exit, if INFO = 0, the N-by-NRHS solution matrix X.
INFO is INTEGER
= 0: successful exit
< 0: if INFO = -i, the i-th argument had an illegal value
> 0: if INFO = i, D(i,i) is exactly zero. The factorization
has been completed, but the block diagonal matrix D is
exactly singular, so the solution could not be computed.
SOURCE
3847 subroutine xhesv_cplex(UPLO, cplex, N, NRHS, A, B, msg, info) 3848 3849 !Arguments ------------------------------------ 3850 character(len=1),intent(in) :: UPLO 3851 integer,intent(in) :: cplex, N, NRHS 3852 real(dp),intent(inout) :: A(cplex, N, N) 3853 real(dp),intent(inout) :: B(cplex, N, NRHS) 3854 character(len=*),intent(out) :: msg 3855 integer,intent(out) :: info 3856 3857 !Local variables ------------------------------ 3858 !scalars 3859 integer :: lwork, lda, ldb 3860 integer,allocatable :: ipiv(:) 3861 real(dp),allocatable :: work(:,:) 3862 3863 !************************************************************************ 3864 3865 if (all(cplex /= [1, 2])) then 3866 write(msg,'(a,i0)')" Wrong value for cplex: ",cplex 3867 info = 1; return 3868 end if 3869 3870 lda = N; ldb = N 3871 3872 ABI_MALLOC(ipiv, (N)) 3873 ABI_MALLOC(work, (cplex, 1)) 3874 lwork = -1 3875 3876 if (cplex == 2) then 3877 ! Complex version 3878 call zhesv(uplo, N, NRHS, A, lda, ipiv, B, ldb, work, lwork, info) 3879 lwork = int(work(1, 1)) 3880 ABI_REMALLOC(work, (cplex, lwork)) 3881 3882 call zhesv(uplo, N, NRHS, A, lda, ipiv, B, ldb, work, lwork, info) 3883 else 3884 ! Read version 3885 call dsysv(uplo, N, NRHS, A, lda, ipiv, B, ldb, work, lwork, info) 3886 lwork = int(work(1, 1)) 3887 ABI_REMALLOC(work, (cplex, lwork)) 3888 3889 call dsysv(uplo, N, NRHS, A, lda, ipiv, B, ldb, work, lwork, info) 3890 end if 3891 3892 ABI_FREE(ipiv) 3893 ABI_FREE(work) 3894 3895 if (info < 0) then 3896 write(msg,'(a,i0,a)')" The ",-info,"-th argument of ZPOTRI had an illegal value." 3897 else if (info > 0) then 3898 write(msg, "(a,i0,4a)") & 3899 " D(i,i) is exactly zero for i= ", info, ch10, & 3900 "The factorization has been completed, but the block diagonal matrix D is ", ch10, & 3901 "exactly singular, so the solution could not be computed." 3902 end if 3903 3904 end subroutine xhesv_cplex
m_hide_lapack/cginv [ Functions ]
[ Top ] [ m_hide_lapack ] [ Functions ]
NAME
cginv
FUNCTION
Invert a general matrix of complex elements in single precision. CGETRF computes an LU factorization of a general N-by-N matrix A using partial pivoting with row interchanges. CGETRI computes the inverse of a matrix using the LU factorization computed by CGETRF.
INPUTS
n=size of complex matrix a a=matrix of complex elements [comm]=MPI communicator for ScaLAPACK inversion. Only available if the code has been compiled with Scalapack support. To avoid wasting CPU time the scalapack initialization is avoided if the number of processors in 1, in this case the sequential LAPACK routine is called.
SIDE EFFECTS
a(n,n)= array of complex elements, input, inverted at output
TODO
Add Scalapack version, matrix_scalapack has to be modified by adding a single precision complex buffer.
SOURCE
2527 subroutine cginv(a, n, comm) 2528 2529 !Arguments ------------------------------------ 2530 !scalars 2531 integer,intent(in) :: n 2532 integer,optional,intent(in) :: comm 2533 !arrays 2534 complex(spc),intent(inout) :: a(n,n) 2535 2536 !Local variables------------------------------- 2537 !scalars 2538 integer :: lwork,info,nprocs 2539 logical :: use_scalapack 2540 character(len=500) :: msg 2541 !arrays 2542 integer,allocatable :: ipiv(:) 2543 complex(spc),allocatable :: work(:) 2544 #ifdef HAVE_LINALG_SCALAPACK 2545 !integer :: ierr,istwf_k,ipiv_size,liwork 2546 !integer,allocatable :: iwork(:) 2547 !type(matrix_scalapack) :: Slk_mat 2548 !type(processor_scalapack) :: Slk_processor 2549 #endif 2550 2551 ! ************************************************************************* 2552 2553 use_scalapack=.FALSE. 2554 if (PRESENT(comm)) then 2555 nprocs = xmpi_comm_size(comm) 2556 ! TODO 2557 !#ifdef HAVE_LINALG_SCALAPACK 2558 ! use_scalapack = (nprocs>1) 2559 !#endif 2560 end if 2561 2562 SELECT CASE(use_scalapack) 2563 2564 CASE (.FALSE.) 2565 ABI_MALLOC(ipiv, (n)) 2566 2567 call CGETRF(n,n,a,n,ipiv,info) ! P* L* U Factorization. 2568 2569 if (info < 0) then 2570 write(msg,'(a,i0,a)')" The ",-info,"-th argument of CGETRF had an illegal value." 2571 ABI_ERROR(msg) 2572 end if 2573 2574 if (info > 0) then 2575 write(msg,'(3a,i0,4a)')& 2576 "The matrix that has been passed in argument is probably either singular or nearly singular.",ch10,& 2577 "U(i,i) in the P*L*U factorization is exactly zero for i = ",info,ch10,& 2578 "The factorization has been completed but the factor U is exactly singular.",ch10,& 2579 "Division by zero will occur if it is used to solve a system of equations." 2580 ABI_ERROR(msg) 2581 end if 2582 2583 lwork=MAX(1,n) 2584 ABI_MALLOC(work,(lwork)) 2585 2586 call CGETRI(n,a,n,ipiv,work,lwork,info) ! Inverts U and the computes inv(A) 2587 2588 if (info < 0) then 2589 write(msg,'(a,i0,a)')" The ",-info,"-th argument of CGETRI had an illegal value." 2590 ABI_ERROR(msg) 2591 end if 2592 2593 if (info > 0) then 2594 write(msg,'(3a,i0,a)')& 2595 "The matrix that has been passed to this subroutine is probably either singular or nearly singular.",ch10,& 2596 "U(i,i) for i= ",info," is exactly zero; the matrix is singular and its inverse could not be computed." 2597 ABI_ERROR(msg) 2598 end if 2599 2600 ABI_FREE(ipiv) 2601 ABI_FREE(work) 2602 RETURN 2603 2604 CASE (.TRUE.) 2605 2606 #if 0 2607 ! FIXME matrix_scalapack does not have a single precision complex buffer 2608 2609 #ifdef HAVE_LINALG_SCALAPACK 2610 call Slk_processor%init(comm) 2611 istwf_k=1 2612 2613 ! Initialize and fill Scalapack matrix from the global one. 2614 call Slk_mat%init(n,n,Slk_processor,istwf_k) 2615 2616 ! IMPORTANT NOTE: PZGETRF requires square block decomposition i.e., MB_A = NB_A. 2617 if ( Slk_mat%descript%tab(MB_)/=Slk_mat%descript%tab(NB_) ) then 2618 msg ="PZGETRF requires square block decomposition i.e., MB_A = NB_A." 2619 ABI_ERROR(msg) 2620 end if 2621 2622 !!call slk_matrix_from_global_dpc_2D(Slk_mat,"All",a) 2623 2624 ipiv_size = my_locr(Slk_mat) + Slk_mat%descript%tab(MB_) 2625 ABI_MALLOC(ipiv,(ipiv_size)) 2626 2627 call PCGETRF(Slk_mat%sizeb_global(1),Slk_mat%sizeb_global(2),Slk_mat%buffer_cplx_sp,& 2628 & 1,1,Slk_mat%descript%tab,ipiv,info) ! P * L * U Factorization. 2629 2630 if (info/=0) then 2631 write(msg,'(a,i0)')"PCGETRF returned info= ",info 2632 ABI_ERROR(msg) 2633 end if 2634 2635 ! Get optimal size of workspace for PCGETRI. 2636 lwork=-1; liwork=-1 2637 ABI_MALLOC(work,(1)) 2638 ABI_MALLOC(iwork,(1)) 2639 2640 call PCGETRI(Slk_mat%sizeb_global(1),Slk_mat%buffer_cplx_sp,1,1,Slk_mat%descript%tab,ipiv,& 2641 & work,lwork,iwork,liwork,info) 2642 2643 ABI_CHECK(info==0,"PZGETRI: Error during compuation of workspace size") 2644 2645 lwork = NINT(DBLE(work(1))); liwork=iwork(1) 2646 ABI_FREE(work) 2647 ABI_FREE(iwork) 2648 2649 ! Solve the problem. 2650 ABI_MALLOC(work,(lwork)) 2651 ABI_MALLOC(iwork,(liwork)) 2652 2653 call PCGETRI(Slk_mat%sizeb_global(1),Slk_mat%buffer_cplx_sp,1,1,Slk_mat%descript%tab,ipiv,& 2654 & work,lwork,iwork,liwork,info) 2655 2656 if (info/=0) then 2657 write(msg,'(a,i0)')"PZGETRI returned info= ",info 2658 ABI_ERROR(msg) 2659 end if 2660 2661 ABI_FREE(work) 2662 ABI_FREE(iwork) 2663 ABI_FREE(ipiv) 2664 2665 ! Reconstruct the global matrix from the distributed one. 2666 a = czero 2667 !! call slk_matrix_to_global_dpc_2D(Slk_mat,"All",a) ! Fill the entries calculated by this node. 2668 call Slk_mat%free() 2669 2670 call xmpi_sum(a,comm,ierr) ! Fill the remaing entries of the global matrix 2671 call Slk_processor%free() 2672 2673 RETURN 2674 #endif 2675 2676 #endif 2677 2678 ABI_BUG("You should not be here!") 2679 2680 END SELECT 2681 2682 end subroutine cginv
m_hide_lapack/dzegdi [ Functions ]
[ Top ] [ m_hide_lapack ] [ Functions ]
NAME
dzgedi
FUNCTION
This routine is the clone of zgefa.F90 using real*8 a(2) instead of complex*16 for the purpose of ABINIT
INPUTS
OUTPUT
SIDE EFFECTS
NOTES
SOURCE
3408 subroutine dzgedi(a,lda,n,ipvt,det,work,job) 3409 3410 integer :: lda,n,ipvt(n),job 3411 real*8 :: a(2,lda,n),det(2,2),work(2,n) 3412 ! 3413 ! zgedi computes the determinant and inverse of a matrix 3414 ! using the factors computed by zgeco or zgefa. 3415 ! 3416 ! on entry 3417 ! 3418 ! a complex*16(lda, n) 3419 ! the output from zgeco or zgefa. 3420 ! 3421 ! lda integer 3422 ! the leading dimension of the array a . 3423 ! 3424 ! n integer 3425 ! the order of the matrix a . 3426 ! 3427 ! ipvt integer(n) 3428 ! the pivot vector from zgeco or zgefa. 3429 ! 3430 ! work complex*16(n) 3431 ! work vector. contents destroyed. 3432 ! 3433 ! job integer 3434 ! = 11 both determinant and inverse. 3435 ! = 01 inverse only. 3436 ! = 10 determinant only. 3437 ! 3438 ! on return 3439 ! 3440 ! a inverse of original matrix if requested. 3441 ! otherwise unchanged. 3442 ! 3443 ! det complex*16(2) 3444 ! determinant of original matrix if requested. 3445 ! otherwise not referenced. 3446 ! determinant = det(1) * 10.0**det(2) 3447 ! with 1.0 .le. cabs1(det(1)) .lt. 10.0 3448 ! or det(1) .eq. 0.0 . 3449 ! 3450 ! error condition 3451 ! 3452 ! a division by zero will occur if the input factor contains 3453 ! a zero on the diagonal and the inverse is requested. 3454 ! it will not occur if the subroutines are called correctly 3455 ! and if zgeco has set rcond .gt. 0.0 or zgefa has set 3456 ! info .eq. 0 . 3457 ! 3458 ! linpack. this version dated 08/14/78 . 3459 ! cleve moler, university of new mexico, argonne national lab. 3460 ! 3461 ! subroutines and functions 3462 ! 3463 ! internal variables 3464 ! 3465 double precision :: r(2),rk(2),rkj(2) 3466 double precision :: ten,rinv2,rabs 3467 integer :: i,j,k,kb,kp1,l,nm1 3468 ! 3469 ! compute determinant 3470 ! 3471 if (job/10 .eq. 0) go to 70 3472 det(1,1) = 1.0d0; det(2,1) = 0.0d0 3473 det(1,2) = 0.0d0; det(2,2) = 0.0d0 3474 ten = 10.0d0 3475 do i = 1, n 3476 if (ipvt(i) .ne. i) then 3477 det(1,1) = -det(1,1) 3478 det(2,1) = -det(2,1) 3479 end if 3480 r(1)=det(1,1); r(2)=det(2,1) 3481 det(1,1) = r(1)*a(1,i,i)-r(2)*a(2,i,i) 3482 det(2,1) = r(2)*a(1,i,i)+r(1)*a(2,i,i) 3483 ! ...exit 3484 rabs = abs(det(1,1))+abs(det(2,1)) 3485 if (rabs .eq. 0.0d0) go to 60 3486 10 continue 3487 rabs = abs(det(1,1))+abs(det(2,1)) 3488 if (rabs .ge. 1.0d0) go to 20 3489 det(1,1) = ten*det(1,1); det(2,1) = ten*det(2,1) 3490 det(1,2) = det(1,2) - 1.0d0 3491 go to 10 3492 20 continue 3493 30 continue 3494 rabs = abs(det(1,1))+abs(det(2,1)) 3495 if (rabs .lt. ten) go to 40 3496 det(1,1) = det(1,1)/ten; det(2,1) = det(2,1)/ten 3497 det(1,2) = det(1,2) + 1.0d0 3498 go to 30 3499 40 continue 3500 end do 3501 60 continue 3502 70 continue 3503 ! 3504 ! compute inverse(u) 3505 ! 3506 if (mod(job,10) .eq. 0) go to 150 3507 do 100 k = 1, n 3508 !a(k,k) = (1.0d0,0.0d0)/a(k,k) 3509 !t = -a(k,k) 3510 !call zscal(k-1,t,a(1,k),1) 3511 rinv2 = 1.d0/(a(1,k,k)**2+a(2,k,k)**2) 3512 a(1,k,k) = rinv2*a(1,k,k) 3513 a(2,k,k) = -rinv2*a(2,k,k) 3514 rk(1) = -a(1,k,k); rk(2) = -a(2,k,k) 3515 do i=1,k-1 3516 r(1)=a(1,i,k) 3517 r(2)=a(2,i,k) 3518 a(1,i,k)=rk(1)*r(1)-rk(2)*r(2) 3519 a(2,i,k)=rk(1)*r(2)+rk(2)*r(1) 3520 end do 3521 kp1 = k + 1 3522 if (n .lt. kp1) go to 90 3523 do 80 j = kp1, n 3524 !t = a(k,j) 3525 !a(k,j) = (0.0d0,0.0d0) 3526 !call zaxpy(k,t,a(1,k),1,a(1,j),1) 3527 rkj(1) = a(1,k,j); rkj(2) = a(2,k,j) 3528 a(1,k,j) = 0.d0; a(2,k,j) = 0.d0 3529 do i=1,k 3530 a(1,i,j)=rkj(1)*a(1,i,k)-rkj(2)*a(2,i,k)+a(1,i,j) 3531 a(2,i,j)=rkj(2)*a(1,i,k)+rkj(1)*a(2,i,k)+a(2,i,j) 3532 end do 3533 80 continue 3534 90 continue 3535 100 continue 3536 do i=1,n 3537 end do 3538 ! 3539 ! form inverse(u)*inverse(l) 3540 ! 3541 nm1 = n - 1 3542 if (nm1 .lt. 1) go to 140 3543 do 130 kb = 1, nm1 3544 k = n - kb 3545 kp1 = k + 1 3546 do 110 i = kp1, n 3547 work(1,i) = a(1,i,k); work(2,i) = a(2,i,k) 3548 a(1,i,k) = 0.0d0; a(2,i,k) = 0.d0 3549 110 continue 3550 do 120 j = kp1, n 3551 r(1) = work(1,j); r(2) = work(2,j) 3552 !call zaxpy(n,t,a(1,j),1,a(1,k),1) 3553 do i=1,n 3554 a(1,i,k)=r(1)*a(1,i,j)-r(2)*a(2,i,j)+a(1,i,k) 3555 a(2,i,k)=r(2)*a(1,i,j)+r(1)*a(2,i,j)+a(2,i,k) 3556 end do 3557 120 continue 3558 l = ipvt(k) 3559 if (l .ne. k) then 3560 !call zswap(n,a(1,k),1,a(1,l),1) 3561 do i=1,n 3562 r(1) = a(1,i,k); r(2) = a(2,i,k) 3563 a(1,i,k) = a(1,i,l); a(2,i,k) = a(2,i,l) 3564 a(1,i,l) = r(1); a(2,i,l) = r(2) 3565 end do 3566 end if 3567 130 continue 3568 140 continue 3569 150 continue 3570 3571 end subroutine dzgedi
m_hide_lapack/dzgefa [ Functions ]
[ Top ] [ m_hide_lapack ] [ Functions ]
NAME
dzgefa
FUNCTION
This routine is the clone of zgefa.F90 using real*8 a(2) instead of complex*16 for the purpose of ABINIT (2008,TD)
INPUTS
OUTPUT
SIDE EFFECTS
NOTES
SOURCE
3592 subroutine dzgefa(a,lda,n,ipvt,info) 3593 3594 use m_linalg_interfaces 3595 3596 !Arguments 3597 integer :: lda,n,ipvt(n),info 3598 real*8 :: a(2,lda,n) 3599 ! 3600 ! zgefa factors a complex*16 matrix by gaussian elimination. 3601 ! 3602 ! dzgefa is usually called by zgeco, but it can be called 3603 ! directly with a saving in time if rcond is not needed. 3604 ! (time for zgeco) = (1 + 9/n)*(time for zgefa) . 3605 ! 3606 ! on entry 3607 ! 3608 ! a complex*16(lda, n) 3609 ! the matrix to be factored. 3610 ! 3611 ! lda integer 3612 ! the leading dimension of the array a . 3613 ! 3614 ! n integer 3615 ! the order of the matrix a . 3616 ! 3617 ! on return 3618 ! 3619 ! a an upper triangular matrix and the multipliers 3620 ! which were used to obtain it. 3621 ! the factorization can be written a = l*u where 3622 ! l is a product of permutation and unit lower 3623 ! triangular matrices and u is upper triangular. 3624 ! 3625 ! ipvt integer(n) 3626 ! an integer vector of pivot indices. 3627 ! 3628 ! info integer 3629 ! = 0 normal value. 3630 ! = k if u(k,k) .eq. 0.0 . this is not an error 3631 ! condition for this subroutine, but it does 3632 ! indicate that zgesl or zgedi will divide by zero 3633 ! if called. use rcond in zgeco for a reliable 3634 ! indication of singularity. 3635 ! 3636 ! linpack. this version dated 08/14/78 . 3637 ! cleve moler, university of new mexico, argonne national lab. 3638 ! 3639 ! subroutines and functions 3640 ! 3641 ! internal variables 3642 ! 3643 !Local variables 3644 real*8 :: r(2),rk(2),rlj(2) 3645 real*8 :: rinv2,rmax,rabs 3646 integer :: i,j,k,kp1,l,nm1 3647 3648 ! 3649 ! gaussian elimination with partial pivoting 3650 ! 3651 info = 0 3652 nm1 = n - 1 3653 if (nm1 .lt. 1) go to 70 3654 do 60 k = 1, nm1 3655 kp1 = k + 1 3656 ! 3657 ! find l = pivot index 3658 ! 3659 !l = izamax(n-k+1,a(k,k),1) + k - 1 3660 rmax=0.d0 3661 l=0 3662 do i=k,n 3663 rabs=abs(a(1,i,k))+abs(a(2,i,k)) 3664 if(rmax<=rabs) then 3665 rmax=rabs 3666 l=i 3667 end if 3668 end do 3669 ipvt(k) = l 3670 ! 3671 ! zero pivot implies this column already triangularized 3672 ! 3673 if (abs(a(1,l,k))+abs(a(2,l,k)) .eq. 0.0d0) go to 40 3674 ! 3675 ! interchange if necessary 3676 ! 3677 if (l .eq. k) go to 10 3678 r(1) = a(1,l,k); r(2) = a(2,l,k) 3679 a(1,l,k) = a(1,k,k); a(2,l,k) = a(2,k,k) 3680 a(1,k,k) = r(1); a(2,k,k) = r(2) 3681 10 continue 3682 ! 3683 ! compute multipliers 3684 ! 3685 rinv2 = 1.d0/(a(1,k,k)**2+a(2,k,k)**2) 3686 rk(1) = -rinv2*a(1,k,k) 3687 rk(2) = rinv2*a(2,k,k) 3688 !call zscal(n-k,t,a(k+1,k),1) 3689 do i=k+1,n 3690 r(1)=a(1,i,k) 3691 r(2)=a(2,i,k) 3692 a(1,i,k)=rk(1)*r(1)-rk(2)*r(2) 3693 a(2,i,k)=rk(1)*r(2)+rk(2)*r(1) 3694 end do 3695 ! 3696 ! row elimination with column indexing 3697 ! 3698 do j = kp1, n 3699 rlj(1) = a(1,l,j); rlj(2) = a(2,l,j) 3700 if (l .eq. k) go to 20 3701 a(1,l,j) = a(1,k,j); a(2,l,j) = a(2,k,j) 3702 a(1,k,j) = rlj(1); a(2,k,j) = rlj(2) 3703 20 continue 3704 !call zaxpy(n-k,t,a(1,k+1,k),1,a(1,k+1,j),1) 3705 do i=k+1,n 3706 a(1,i,j)=rlj(1)*a(1,i,k)-rlj(2)*a(2,i,k)+a(1,i,j) 3707 a(2,i,j)=rlj(2)*a(1,i,k)+rlj(1)*a(2,i,k)+a(2,i,j) 3708 end do 3709 end do 3710 go to 50 3711 40 continue 3712 info = k 3713 50 continue 3714 60 continue 3715 70 continue 3716 ipvt(n) = n 3717 if (abs(a(1,n,n))+abs(a(2,n,n)) .eq. 0.0d0) info = n 3718 3719 end subroutine dzgefa
m_hide_lapack/lubksb [ Functions ]
[ Top ] [ m_hide_lapack ] [ Functions ]
NAME
lubksb
FUNCTION
Solves the set of n linear equations A . X = B. Here a is input, not as the matrix A but rather as its LU decomposition, determined by the routine ludcmp. indx is input as the permutation vector returned by ludcmp. b(1:n) is input as the right-hand side vector B, and returns with the solution vector X. a, n, np, and indx are not modified by this routine and can be left in place for successive calls with different right-hand sides b. This routine takes into account the possibility that b will begin with many zero elements, so it is efficient for use in matrix inversion.
INPUTS
OUTPUT
SIDE EFFECTS
NOTES
This routine is deprecated, use lapack API
SOURCE
3352 SUBROUTINE lubksb(a,n,np,indx,b) 3353 3354 INTEGER n,np,indx(n) 3355 REAL*8 a(np,np),b(n) 3356 3357 INTEGER i,ii,j,ll 3358 REAL*8 sum 3359 ! write(std_out,*) 'ENTERING LUBKSB...' 3360 ! write(std_out,201) ((a(i,j),j=1,n),i=1,n) 3361 ! 201 FORMAT('A in lubksb ',/,3F16.8,/,3F16.8,/,3F16.8) 3362 3363 ii=0 3364 do i=1,n 3365 ll=indx(i) 3366 sum=b(ll) 3367 b(ll)=b(i) 3368 if (ii.ne.0)then 3369 do j=ii,i-1 3370 sum=sum-a(i,j)*b(j) 3371 enddo 3372 else if (sum.ne.0.) then 3373 ii=i 3374 endif 3375 b(i)=sum 3376 enddo 3377 do i=n,1,-1 3378 sum=b(i) 3379 do j=i+1,n 3380 sum=sum-a(i,j)*b(j) 3381 enddo 3382 b(i)=sum/a(i,i) 3383 enddo 3384 ! write(std_out,*) 'LEAVING LUBKSB...' 3385 return 3386 3387 END SUBROUTINE LUBKSB
m_hide_lapack/ludcmp [ Functions ]
[ Top ] [ m_hide_lapack ] [ Functions ]
NAME
ludcmp
FUNCTION
Given a matrix a(1:n,1:n), with physical dimension np by np, this routine replaces it by the LU decomposition of a rowwise permutation of itself. a and n are input. a is output, arranged as in equation (2.3.14) above; indx(1:n) is an output vector that records the row permutation effected by the partial pivoting; id is output as +- 1 depending on whether the number of row interchanges was even or odd, respectively. This routine is used in combination with lubksb to solve linear equations or invert a matrix.
INPUTS
OUTPUT
NOTES
This routine is depreacted, use lapack API
SOURCE
3251 SUBROUTINE ludcmp(a,n,np,indx,id,info) 3252 3253 INTEGER n,np,indx(n),NMAX,id,info 3254 REAL*8 a(np,np),TINY 3255 PARAMETER (NMAX=500,TINY=1.0e-20) 3256 3257 INTEGER i,imax,j,k 3258 REAL*8 aamax,dum,sum,vv(NMAX) 3259 3260 ! write(std_out,*) 'ENTERING LUDCMP...' 3261 ! write(std_out,*) 'in ludcmp n=',n,' np=',np 3262 ! write(std_out,201) ((a(i,j),j=1,n),i=1,n) 3263 ! 201 FORMAT('A in ludcmp ',/,3F16.8,/,3F16.8,/,3F16.8) 3264 id=1 3265 info=0 3266 do i=1,n 3267 aamax=0. 3268 do j=1,n 3269 if (abs(a(i,j)).gt.aamax) aamax=abs(a(i,j)) 3270 enddo 3271 if (aamax.eq.0.) then 3272 write(std_out,*) 'LUDCMP: singular matrix !!!' 3273 do j=1,3 3274 write(std_out,*) (a(j,k),k=1,3) 3275 enddo 3276 info=1 3277 return 3278 ! stop 'singular matrix in ludcmp' 3279 endif 3280 vv(i)=1./aamax 3281 enddo 3282 do j=1,n 3283 do i=1,j-1 3284 sum=a(i,j) 3285 do k=1,i-1 3286 sum=sum-a(i,k)*a(k,j) 3287 enddo 3288 a(i,j)=sum 3289 enddo 3290 aamax=0. 3291 do i=j,n 3292 sum=a(i,j) 3293 do k=1,j-1 3294 sum=sum-a(i,k)*a(k,j) 3295 enddo 3296 a(i,j)=sum 3297 dum=vv(i)*abs(sum) 3298 if (dum.ge.aamax) then 3299 imax=i 3300 aamax=dum 3301 endif 3302 enddo 3303 if (j.ne.imax)then 3304 do k=1,n 3305 dum=a(imax,k) 3306 a(imax,k)=a(j,k) 3307 a(j,k)=dum 3308 enddo 3309 id=-id 3310 vv(imax)=vv(j) 3311 endif 3312 indx(j)=imax 3313 if(a(j,j).eq.0.)a(j,j)=TINY 3314 if(j.ne.n)then 3315 dum=1./a(j,j) 3316 do i=j+1,n 3317 a(i,j)=a(i,j)*dum 3318 enddo 3319 endif 3320 enddo 3321 ! write(std_out,*) 'LEAVING LUDCMP...' 3322 return 3323 END SUBROUTINE ludcmp
m_hide_lapack/matr3eigval [ Functions ]
[ Top ] [ m_hide_lapack ] [ Functions ]
NAME
matr3eigval
FUNCTION
Find the eigenvalues of a real symmetric 3x3 matrix, entered in full storage mode.
INPUTS
matr(3,3)=real symmetric 3x3 matrix
OUTPUT
eigval(3)=three eigenvalues
SOURCE
3080 subroutine matr3eigval(eigval,matr) 3081 3082 !Arguments ------------------------------------ 3083 !arrays 3084 real(dp),intent(in) :: matr(3,3) 3085 real(dp),intent(out) :: eigval(3) 3086 3087 !Local variables------------------------------- 3088 !scalars 3089 integer :: ier 3090 !arrays 3091 real(dp) :: eigvec(2,3,3),matrx(2,6),zhpev1(2,2*3-1),zhpev2(3*3-2) 3092 3093 ! ************************************************************************* 3094 3095 matrx(1,1)=matr(1,1) 3096 matrx(1,2)=matr(1,2) 3097 matrx(1,3)=matr(2,2) 3098 matrx(1,4)=matr(1,3) 3099 matrx(1,5)=matr(2,3) 3100 matrx(1,6)=matr(3,3) 3101 matrx(2,:)=zero 3102 3103 call ZHPEV ('V','U',3,matrx,eigval,eigvec,3,zhpev1,zhpev2,ier) 3104 !write(std_out,*)' eigval=',eigval 3105 3106 end subroutine matr3eigval
m_hide_lapack/matrginv [ Functions ]
[ Top ] [ m_hide_lapack ] [ Functions ]
NAME
matrginv
FUNCTION
Invert a general matrix of real*8 elements.
INPUTS
lda=leading dimension of complex matrix a n=size of complex matrix a a=matrix of real elements
OUTPUT
a=inverse of a input matrix
SIDE EFFECTS
a(lda,n)= array of real elements, input, inverted at output
SOURCE
2964 subroutine matrginv(a,lda,n) 2965 2966 !Arguments ------------------------------------ 2967 !scalars 2968 integer,intent(in) :: lda,n 2969 !arrays 2970 real(dp),intent(inout) :: a(lda,n) 2971 2972 !Local variables------------------------------- 2973 !scalars 2974 integer :: ierr,nwork 2975 #if defined HAVE_LINALG_ESSL 2976 real(dp) :: rcond 2977 #endif 2978 character(len=500) :: message 2979 !arrays 2980 integer,allocatable :: ipvt(:) 2981 #if defined HAVE_LINALG_ESSL 2982 real(dp) :: det(2) 2983 #elif defined HAVE_LINALG_ASL 2984 real(dp) :: det(2) 2985 #endif 2986 real(dp),allocatable :: work(:) 2987 2988 ! ************************************************************************* 2989 2990 #if defined HAVE_LINALG_ESSL 2991 nwork=200*n 2992 #else 2993 nwork=n 2994 #endif 2995 2996 ABI_MALLOC(work,(nwork)) 2997 ABI_MALLOC(ipvt,(n)) 2998 2999 #if defined HAVE_LINALG_ESSL 3000 3001 call dgeicd(a,lda,n,0,rcond,det,work,nwork) 3002 if(abs(rcond)==zero) then 3003 write(message, '(7a)' )& 3004 ' The matrix that has been passed in argument of this subroutine',ch10,& 3005 ' is probably either singular or nearly singular.',ch10,& 3006 ' The ESSL routine dgeicd failed.',ch10,& 3007 ' Action: Contact ABINIT group ' 3008 ABI_ERROR(message) 3009 end if 3010 3011 #elif defined HAVE_LINALG_ASL 3012 3013 call dbgmlu(a,lda,n,ipvt,ierr) 3014 if(ierr /= 0) then 3015 write(message, '(7a)' ) ch10,& 3016 ' The matrix that has been passed in argument of this subroutine',ch10,& 3017 ' is probably either singular or nearly singular.',ch10,& 3018 ' The ASL routine dbgmlu failed.',ch10,& 3019 ' Action: Contact ABINIT group ' 3020 ABI_ERROR(message) 3021 end if 3022 3023 call dbgmdi(a,lda,n,ipvt,det,-1,work,ierr) 3024 3025 if(ierr /= 0) then 3026 write(message, '(7a)' ) & 3027 ' The matrix that has been passed in argument of this subroutine',ch10,& 3028 ' is probably either singular or nearly singular.',ch10,& 3029 ' The ASL routine dbgmdi failed.',ch10,& 3030 ' Action: Contact ABINIT group ' 3031 ABI_ERROR(message) 3032 end if 3033 3034 #else 3035 3036 call dgetrf(n,n,a,lda,ipvt,ierr) 3037 if(ierr /= 0) then 3038 write(message, '(7a)' ) & 3039 ' The matrix that has been passed in argument of this subroutine',ch10,& 3040 ' is probably either singular or nearly singular.',ch10,& 3041 ' The LAPACK routine dgetrf failed.',ch10,& 3042 ' Action: Contact ABINIT group ' 3043 ABI_ERROR(message) 3044 end if 3045 3046 call dgetri(n,a,lda,ipvt,work,n,ierr) 3047 3048 if(ierr /= 0) then 3049 write(message, '(7a)' ) & 3050 ' The matrix that has been passed in argument of this subroutine',ch10,& 3051 ' is probably either singular or nearly singular.',ch10,& 3052 ' The LAPACK routine dgetri failed.',ch10,& 3053 ' Action: Contact ABINIT group ' 3054 ABI_ERROR(message) 3055 end if 3056 3057 #endif 3058 3059 ABI_FREE(work) 3060 ABI_FREE(ipvt) 3061 3062 end subroutine matrginv
m_hide_lapack/test_xginv [ Functions ]
[ Top ] [ m_hide_lapack ] [ Functions ]
NAME
test_xginv
FUNCTION
m_hide_lapack/wrap_CGEEV [ Functions ]
[ Top ] [ m_hide_lapack ] [ Functions ]
NAME
wrap_CGEEV
FUNCTION
wrap_CGEEV computes for an N-by-N complex nonsymmetric matrix A, the eigenvalues and, optionally, the left and/or right eigenvectors using single precision arithmetic. [PRIVATE] The right eigenvector v(j) of A satisfies: A * v(j) = lambda(j) * v(j) where lambda(j) is its eigenvalue. The left eigenvector u(j) of A satisfies u(j)**H * A = lambda(j) * u(j)**H where u(j)**H denotes the conjugate transpose of u(j). The computed eigenvectors are normalized to have Euclidean norm equal to 1 and largest component real.
INPUTS
JOBVL (input) CHARACTER*1
= 'N': left eigenvectors of A are not computed;
= 'V': left eigenvectors of are computed.
JOBVR (input) CHARACTER*1
= 'N': right eigenvectors of A are not computed;
= 'V': right eigenvectors of A are computed.
N (input) INTEGER
The order of the matrix A. N >= 0.
LDA (input) INTEGER
The leading dimension of the array A. LDA >= max(1,N).
LDVL (input) INTEGER
The leading dimension of the array VL. LDVL >= 1; if
JOBVL = 'V', LDVL >= N.
LDVR (input) INTEGER
The leading dimension of the array VR. LDVR >= 1; if
JOBVR = 'V', LDVR >= N.
OUTPUT
W (output) COMPLEX(SPC) array, dimension (N)
W contains the computed eigenvalues.
VL (output) COMPLEX(SCP) array, dimension (LDVL,N)
If JOBVL = 'V', the left eigenvectors u(j) are stored one
after another in the columns of VL, in the same order
as their eigenvalues.
If JOBVL = 'N', VL is not referenced.
u(j) = VL(:,j), the j-th column of VL.
VR (output) COMPLEX(SPC) array, dimension (LDVR,N)
If JOBVR = 'V', the right eigenvectors v(j) are stored one
after another in the columns of VR, in the same order
as their eigenvalues.
If JOBVR = 'N', VR is not referenced.
v(j) = VR(:,j), the j-th column of VR.
See also SIDE EFFECTS
SIDE EFFECTS
A (input/output) COMPLEX(SPC) array, dimension (LDA,N)
On entry, the N-by-N matrix A.
On exit, A has been overwritten.
SOURCE
2330 subroutine wrap_CGEEV(jobvl, jobvr, n, a, lda, w, vl, ldvl, vr, ldvr) 2331 2332 !Arguments ------------------------------------ 2333 !scalars 2334 integer,intent(in) :: n,lda,ldvl,ldvr 2335 character(len=*),intent(in) :: jobvl,jobvr 2336 !arrays 2337 complex(spc),intent(inout) :: a(lda,n) 2338 complex(spc),intent(out) :: w(n) 2339 complex(spc),intent(out) :: vl(ldvl,n) 2340 complex(spc),intent(out) :: vr(ldvr,n) 2341 2342 !Local variables ------------------------------ 2343 !scalars 2344 integer :: info,lwork 2345 character(len=500) :: msg 2346 !arrays 2347 real(sp),allocatable :: rwork(:) 2348 complex(spc),allocatable :: work(:) 2349 2350 !************************************************************************ 2351 2352 lwork = MAX(1,2*n) 2353 2354 ABI_MALLOC(work,(lwork)) 2355 ABI_MALLOC(rwork,(2*n)) 2356 2357 call CGEEV(jobvl,jobvr,n,a,lda,w,vl,ldvl,vr,ldvr,work,lwork,rwork,info) 2358 2359 if (info < 0) then 2360 write(msg,'(a,i0,a)')" The ",-info,"-th argument of CGEEV had an illegal value." 2361 ABI_ERROR(msg) 2362 end if 2363 2364 if (info > 0) then 2365 write(msg,'(3a,i0,a,i0,a)')& 2366 "CGEEV: The QR algorithm failed to compute all the eigenvalues, and no eigenvectors have been computed;",ch10,& 2367 "Elements ",info+1,":",n," of W contain eigenvalues which have converged. " 2368 ABI_ERROR(msg) 2369 end if 2370 2371 ABI_FREE(work) 2372 ABI_FREE(rwork) 2373 2374 end subroutine wrap_CGEEV
m_hide_lapack/wrap_CHEEV [ Functions ]
[ Top ] [ m_hide_lapack ] [ Functions ]
NAME
wrap_CHEEV
FUNCTION
wrap_CHEEV computes the eigenvalues and, optionally, the eigenvectors of a complex Hermitian matrix in single precision. [PRIVATE]
INPUTS
JOBZ (input) CHARACTER*1
= 'N': Compute eigenvalues only;
= 'V': Compute eigenvalues and eigenvectors.
UPLO (input) CHARACTER*1
= 'U': Upper triangle of A is stored;
= 'L': Lower triangle of A is stored.
N (input) INTEGER
The order of the matrix A. N >= 0.
OUTPUT
W (output) REAL(SP) array, dimension (N)
If INFO = 0, the eigenvalues in ascending order.
See also SIDE EFFECTS
SIDE EFFECTS
A (input/output) COMPLEX(SPC) array, dimension (N, N)
On entry, the Hermitian matrix A. If UPLO = 'U', the
leading N-by-N upper triangular part of A contains the
upper triangular part of the matrix A. If UPLO = 'L',
the leading N-by-N lower triangular part of A contains
the lower triangular part of the matrix A.
On exit, if JOBZ = 'V', then if INFO = 0, A contains the
orthonormal eigenvectors of the matrix A.
If JOBZ = 'N', then on exit the lower triangle (if UPLO='L')
or the upper triangle (if UPLO='U') of A, including the
diagonal, is destroyed.
SOURCE
224 subroutine wrap_CHEEV(jobz, uplo, n, a, w) 225 226 !Arguments ------------------------------------ 227 !scalars 228 integer,intent(in) :: n 229 character(len=*),intent(in) :: jobz,uplo 230 !scalars 231 real(sp),intent(out) :: w(n) 232 complex(spc),intent(inout) :: a(n,n) 233 234 !Local variables ------------------------------ 235 !scalars 236 integer :: lwork,info 237 character(len=500) :: msg 238 !arrays 239 real(sp),allocatable :: rwork(:) 240 complex(spc),allocatable :: work(:) 241 242 !************************************************************************ 243 244 lwork = MAX(1,2*n-1) 245 246 ABI_MALLOC(work, (lwork)) 247 ABI_MALLOC(rwork, (MAX(1,3*n-2))) 248 249 call CHEEV(jobz,uplo,n,a,n,w,work,lwork,rwork,info) 250 251 if (info < 0) then 252 write(msg,'(a,i0,a)')"The ",-info,"-th argument of CHEEV had an illegal value." 253 ABI_ERROR(msg) 254 end if 255 256 if (info > 0) then 257 write(msg,'(2a,i0,a)')& 258 "CHEEV: the algorithm failed to converge; ",ch10,& 259 info," off-diagonal elements of an intermediate tridiagonal form did not converge to zero. " 260 ABI_ERROR(msg) 261 end if 262 263 ABI_FREE(rwork) 264 ABI_FREE(work) 265 266 !TODO scaLAPACK version (complex single precision buffer is needed in matrix_scalapack) 267 268 end subroutine wrap_CHEEV
m_hide_lapack/wrap_CHPEV [ Functions ]
[ Top ] [ m_hide_lapack ] [ Functions ]
NAME
wrap_CHPEV
FUNCTION
wrap_CHPEV computes all the eigenvalues and, optionally, eigenvectors of a complex Hermitian matrix in packed storage. Scalapack version is not available. [PRIVATE].
INPUTS
JOBZ (input) CHARACTER*1
= 'N': Compute eigenvalues only;
= 'V': Compute eigenvalues and eigenvectors.
UPLO (input) CHARACTER*1
= 'U': Upper triangle of A is stored;
= 'L': Lower triangle of A is stored.
N (input) INTEGER
The order of the matrix A. N >= 0.
LDZ (input) INTEGER
The leading dimension of the array Z. LDZ >= 1, and if
JOBZ = 'V', LDZ >= max(1,N).
OUTPUT
W (output) REAL(SP) array, dimension (N)
If INFO = 0, the eigenvalues in ascending order.
Z (output) COMPLEX(SPC) array, dimension (LDZ, N)
If JOBZ = 'V', then if INFO = 0, Z contains the orthonormal
eigenvectors of the matrix A, with the i-th column of Z
holding the eigenvector associated with W(i).
If JOBZ = 'N', then Z is not referenced.
See also SIDE EFFECTS
SIDE EFFECTS
AP (input/output) COMPLEX(SPC) array, dimension (N*(N+1)/2)
On entry, the upper or lower triangle of the Hermitian matrix
A, packed columnwise in a linear array. The j-th column of A
is stored in the array AP as follows:
if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n.
On exit, AP is overwritten by values generated during the
reduction to tridiagonal form. If UPLO = 'U', the diagonal
and first superdiagonal of the tridiagonal matrix T overwrite
the corresponding elements of A, and if UPLO = 'L', the
diagonal and first subdiagonal of T overwrite the
corresponding elements of A.
SOURCE
648 subroutine wrap_CHPEV(jobz, uplo, n, ap, w, z, ldz) 649 650 !Arguments ------------------------------------ 651 !scalars 652 integer,intent(in) :: n,ldz 653 character(len=*),intent(in) :: jobz,uplo 654 !arrays 655 real(sp),intent(out) :: w(n) 656 complex(spc),intent(inout) :: ap(n*(n+1)/2) 657 complex(spc),intent(out) :: z(ldz,n) 658 659 !Local variables ------------------------------ 660 !scalars 661 integer :: info 662 character(len=500) :: msg 663 !arrays 664 real(sp),allocatable :: rwork(:) 665 complex(spc),allocatable :: work(:) 666 667 !************************************************************************ 668 669 ABI_MALLOC(work, (MAX(1,2*n-1))) 670 ABI_MALLOC(rwork, (MAX(1,3*n-2))) 671 672 call CHPEV( JOBZ, UPLO, N, AP, W, Z, LDZ, WORK, RWORK, INFO ) 673 674 if (info < 0) then 675 write(msg,'(a,i0,a)')" The ",-info,"-th argument of ZHEEV had an illegal value." 676 ABI_ERROR(msg) 677 end if 678 679 if (info > 0) then 680 write(msg,'(2a,i0,a)')& 681 "ZHPEV: the algorithm failed to converge; ",ch10,& 682 info," off-diagonal elements of an intermediate tridiagonal form did not converge to zero. " 683 ABI_ERROR(msg) 684 end if 685 686 ABI_FREE(rwork) 687 ABI_FREE(work) 688 689 end subroutine wrap_CHPEV
m_hide_lapack/wrap_ZGEEV [ Functions ]
[ Top ] [ m_hide_lapack ] [ Functions ]
NAME
wrap_ZGEEV
FUNCTION
wrap_ZGEEV computes for an N-by-N complex nonsymmetric matrix A, the eigenvalues and, optionally, the left and/or right eigenvectors using double precision arithmetic. [PRIVATE] The right eigenvector v(j) of A satisfies: A * v(j) = lambda(j) * v(j) where lambda(j) is its eigenvalue. The left eigenvector u(j) of A satisfies u(j)**H * A = lambda(j) * u(j)**H where u(j)**H denotes the conjugate transpose of u(j). The computed eigenvectors are normalized to have Euclidean norm equal to 1 and largest component real. No scalapack version is available (PZGEEV is not provided by the Scalapack team)
INPUTS
JOBVL (input) CHARACTER*1
= 'N': left eigenvectors of A are not computed;
= 'V': left eigenvectors of are computed.
JOBVR (input) CHARACTER*1
= 'N': right eigenvectors of A are not computed;
= 'V': right eigenvectors of A are computed.
N (input) INTEGER
The order of the matrix A. N >= 0.
LDA (input) INTEGER
The leading dimension of the array A. LDA >= max(1,N).
LDVL (input) INTEGER
The leading dimension of the array VL. LDVL >= 1; if
JOBVL = 'V', LDVL >= N.
LDVR (input) INTEGER
The leading dimension of the array VR. LDVR >= 1; if
JOBVR = 'V', LDVR >= N.
OUTPUT
W (output) COMPLEX(DPC) array, dimension (N)
W contains the computed eigenvalues.
VL (output) COMPLEX(DPC) array, dimension (LDVL,N)
If JOBVL = 'V', the left eigenvectors u(j) are stored one
after another in the columns of VL, in the same order
as their eigenvalues.
If JOBVL = 'N', VL is not referenced.
u(j) = VL(:,j), the j-th column of VL.
VR (output) COMPLEX(DPC) array, dimension (LDVR,N)
If JOBVR = 'V', the right eigenvectors v(j) are stored one
after another in the columns of VR, in the same order
as their eigenvalues.
If JOBVR = 'N', VR is not referenced.
v(j) = VR(:,j), the j-th column of VR.
See also SIDE EFFECTS
SIDE EFFECTS
A (input/output) COMPLEX(DPC) array, dimension (LDA,N)
On entry, the N-by-N matrix A.
On exit, A has been overwritten.
SOURCE
2444 subroutine wrap_ZGEEV(jobvl,jobvr,n,a,lda,w,vl,ldvl,vr,ldvr) 2445 2446 !Arguments ------------------------------------ 2447 !scalars 2448 integer,intent(in) :: n,lda,ldvl,ldvr 2449 character(len=*),intent(in) :: jobvl,jobvr 2450 !arrays 2451 complex(dpc),intent(inout) :: a(lda,n) 2452 complex(dpc),intent(out) :: w(n) 2453 complex(dpc),intent(out) :: vl(ldvl,n) 2454 complex(dpc),intent(out) :: vr(ldvr,n) 2455 2456 !Local variables ------------------------------ 2457 !scalars 2458 integer :: info,lwork 2459 logical :: use_scalapack 2460 character(len=500) :: msg 2461 !arrays 2462 real(dp),allocatable :: rwork(:) 2463 complex(dpc),allocatable :: work(:) 2464 2465 !************************************************************************ 2466 2467 use_scalapack=.FALSE. 2468 2469 SELECT CASE(use_scalapack) 2470 CASE (.FALSE.) 2471 2472 lwork = MAX(1,2*n) 2473 ABI_MALLOC(work,(lwork)) 2474 ABI_MALLOC(rwork,(2*n)) 2475 2476 call ZGEEV(jobvl,jobvr,n,a,lda,w,vl,ldvl,vr,ldvr,work,lwork,rwork,info) 2477 2478 if (info < 0) then 2479 write(msg,'(a,i0,a)')" The ",-info,"-th argument of ZGEEV had an illegal value." 2480 ABI_ERROR(msg) 2481 end if 2482 2483 if (info > 0) then 2484 write(msg,'(3a,i0,a,i0,a)')& 2485 "ZGEEV: The QR algorithm failed to compute all the eigenvalues, and no eigenvectors have been computed;",ch10,& 2486 "Elements ",info+1,":",n," of W contain eigenvalues which have converged. " 2487 ABI_ERROR(msg) 2488 end if 2489 2490 ABI_FREE(work) 2491 ABI_FREE(rwork) 2492 RETURN 2493 2494 CASE (.TRUE.) 2495 ABI_BUG("You should not be here!") 2496 END SELECT 2497 2498 end subroutine wrap_ZGEEV
m_hide_lapack/wrap_ZHEEV [ Functions ]
[ Top ] [ m_hide_lapack ] [ Functions ]
NAME
wrap_ZHEEV
FUNCTION
wrap_ZHEEV computes the eigenvalues and, optionally, the eigenvectors of a complex Hermitian matrix in double precision. [PRIVATE]
INPUTS
JOBZ (input) CHARACTER*1
= 'N': Compute eigenvalues only;
= 'V': Compute eigenvalues and eigenvectors.
UPLO (input) CHARACTER*1
= 'U': Upper triangle of A is stored;
= 'L': Lower triangle of A is stored.
N (input) INTEGER
The order of the matrix A. N >= 0.
[comm]=MPI communicator for ScaLAPACK inversion. Only available if the code has been compiled with Scalapack support.
To avoid wasting CPU time the scalapack initialization is avoided if the number of processors in 1,
in this case the sequential LAPACK routine is called.
OUTPUT
W (output) REAL(DP) array, dimension (N)
If INFO = 0, the eigenvalues in ascending order.
See also SIDE EFFECTS
SIDE EFFECTS
A (input/output) COMPLEX(DPC) array, dimension (N, N)
On entry, the Hermitian matrix A. If UPLO = 'U', the
leading N-by-N upper triangular part of A contains the
upper triangular part of the matrix A. If UPLO = 'L',
the leading N-by-N lower triangular part of A contains
the lower triangular part of the matrix A.
On exit, if JOBZ = 'V', then if INFO = 0, A contains the
orthonormal eigenvectors of the matrix A.
If JOBZ = 'N', then on exit the lower triangle (if UPLO='L')
or the upper triangle (if UPLO='U') of A, including the
diagonal, is destroyed.
SOURCE
317 subroutine wrap_ZHEEV(jobz, uplo, n, a, w, comm) 318 319 !Arguments ------------------------------------ 320 !scalars 321 integer,intent(in) :: n 322 integer,optional,intent(in) :: comm 323 character(len=*),intent(in) :: jobz,uplo 324 !arrays 325 complex(dpc),intent(inout) :: a(n,n) 326 real(dp),intent(out) :: w(n) 327 328 !Local variables ------------------------------ 329 !scalars 330 integer :: lwork,info,nprocs 331 logical :: use_scalapack 332 character(len=500) :: msg 333 !arrays 334 real(dp),allocatable :: rwork(:) 335 complex(dpc),allocatable :: work(:) 336 #ifdef HAVE_LINALG_SCALAPACK 337 integer :: ierr,istwf_k 338 logical :: want_eigenvectors 339 type(matrix_scalapack) :: Slk_mat,Slk_vec 340 type(processor_scalapack) :: Slk_processor 341 #endif 342 !************************************************************************ 343 344 use_scalapack=.FALSE. 345 if (PRESENT(comm)) then 346 nprocs = xmpi_comm_size(comm) 347 #ifdef HAVE_LINALG_SCALAPACK 348 use_scalapack = (nprocs>1) 349 #endif 350 end if 351 352 SELECT CASE(use_scalapack) 353 CASE (.FALSE.) 354 355 lwork = MAX(1,2*n-1) 356 ABI_MALLOC(work, (lwork)) 357 ABI_MALLOC(rwork, (MAX(1,3*n-2))) 358 359 call ZHEEV(jobz,uplo,n,a,n,w,work,lwork,rwork,info) 360 361 if (info < 0) then 362 write(msg,'(a,i0,a)')" The ",-info,"-th argument of ZHEEV had an illegal value." 363 ABI_ERROR(msg) 364 end if 365 366 if (info > 0) then 367 write(msg,'(2a,i0,a)')& 368 "ZHEEV: the algorithm failed to converge; ",ch10,& 369 info," off-diagonal elements of an intermediate tridiagonal form did not converge to zero. " 370 ABI_ERROR(msg) 371 end if 372 373 ABI_FREE(rwork) 374 ABI_FREE(work) 375 RETURN 376 377 CASE (.TRUE.) 378 #ifdef HAVE_LINALG_SCALAPACK 379 call Slk_processor%init(comm) 380 istwf_k=1 381 382 ! Initialize and fill Scalapack matrix from the global one. 383 call Slk_mat%init(n,n,Slk_processor,istwf_k) 384 call slk_matrix_from_global_dpc_2D(Slk_mat,uplo,a) 385 386 want_eigenvectors = firstchar(jobz,(/"V","v"/)) 387 if (want_eigenvectors) then 388 ! Initialize the distributed vectors. 389 call Slk_vec%init(n,n,Slk_processor,istwf_k) 390 end if 391 392 ! Solve the problem with scaLAPACK. 393 call slk_mat%heev(jobz, uplo, Slk_vec, w) 394 call Slk_mat%free() 395 396 if (want_eigenvectors) then ! A is overwritten with the eigenvectors 397 a = czero 398 call slk_matrix_to_global_dpc_2D(Slk_vec,"All",a) ! Fill the entries calculated by this node. 399 call Slk_vec%free() 400 call xmpi_sum(a,comm,ierr) ! Fill the remaing entries of the global matrix 401 end if 402 403 call Slk_processor%free() 404 RETURN 405 #endif 406 407 ABI_BUG("You should not be here!") 408 END SELECT 409 410 end subroutine wrap_ZHEEV
m_hide_lapack/wrap_ZHEEVX [ Functions ]
[ Top ] [ m_hide_lapack ] [ Functions ]
NAME
wrap_ZHEEVX
FUNCTION
wrap_ZHEEVX computes selected eigenvalues and, optionally, eigenvectors of a complex Hermitian matrix A. Eigenvalues and eigenvectors can be selected by specifying either a range of values or a range of indices for the desired eigenvalues.
INPUTS
JOBZ (input) CHARACTER*1
= 'N': Compute eigenvalues only;
= 'V': Compute eigenvalues and eigenvectors.
RANGE (input) CHARACTER*1
= 'A': all eigenvalues will be found.
= 'V': all eigenvalues in the half-open interval (VL,VU]
will be found.
= 'I': the IL-th through IU-th eigenvalues will be found.
UPLO (input) CHARACTER*1
= 'U': Upper triangle of A is stored;
= 'L': Lower triangle of A is stored.
N (input) INTEGER
The order of the matrix A. N >= 0.
LDA (input) INTEGER
The leading dimension of the array A. LDA >= max(1,N).
VL (input) REAL(DP)
VU (input) REAL(DP)
If RANGE='V', the lower and upper bounds of the interval to
be searched for eigenvalues. VL < VU.
Not referenced if RANGE = 'A' or 'I'.
IL (input) INTEGER
IU (input) INTEGER
If RANGE='I', the indices (in ascending order) of the
smallest and largest eigenvalues to be returned.
1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
Not referenced if RANGE = 'A' or 'V'.
ABSTOL (input) REAL(DP)
The absolute error tolerance for the eigenvalues.
An approximate eigenvalue is accepted as converged
when it is determined to lie in an interval [a,b]
of width less than or equal to
ABSTOL + EPS * max( |a|,|b| ) ,
where EPS is the machine precision. If ABSTOL is less than
or equal to zero, then EPS*|T| will be used in its place,
where |T| is the 1-norm of the tridiagonal matrix obtained
by reducing A to tridiagonal form.
Eigenvalues will be computed most accurately when ABSTOL is
set to twice the underflow threshold 2*DLAMCH('S'), not zero.
If this routine returns with INFO>0, indicating that some
eigenvectors did not converge, try setting ABSTOL to
2*DLAMCH('S').
See "Computing Small Singular Values of Bidiagonal Matrices
with Guaranteed High Relative Accuracy," by Demmel and
Kahan, LAPACK Working Note #3.
LDZ (input) INTEGER
The leading dimension of the array Z. LDZ >= 1, and if
JOBZ = 'V', LDZ >= max(1,N).
[comm]=MPI communicator for ScaLAPACK inversion. Only available if the code has been compiled with Scalapack support.
To avoid wasting CPU time the scalapack initialization is avoided if the number of processors in 1,
in this case the sequential LAPACK routine is called.
OUTPUT
M (output) INTEGER
The total number of eigenvalues found. 0 <= M <= N.
If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
W (output) REAL(DP) array, dimension (N)
On normal exit, the first M elements contain the selected
eigenvalues in ascending order.
Z (output) COMPLEX(DPC) array, dimension (LDZ, max(1,M))
If JOBZ = 'V', then if INFO = 0, the first M columns of Z
contain the orthonormal eigenvectors of the matrix A
corresponding to the selected eigenvalues, with the i-th
column of Z holding the eigenvector associated with W(i).
If an eigenvector fails to converge, then that column of Z
contains the latest approximation to the eigenvector, and the
index of the eigenvector is returned in IFAIL.
If JOBZ = 'N', then Z is not referenced.
Note: the user must ensure that at least max(1,M) columns are
supplied in the array Z; if RANGE = 'V', the exact value of M
is not known in advance and an upper bound must be used.
See also SIDE EFFECTS
SIDE EFFECTS
A (input/output) COMPLEX(DPC) array, dimension (N, N)
On entry, the Hermitian matrix A. If UPLO = 'U', the
leading N-by-N upper triangular part of A contains the
upper triangular part of the matrix A. If UPLO = 'L',
the leading N-by-N lower triangular part of A contains
the lower triangular part of the matrix A.
On exit, the lower triangle (if UPLO='L') or the upper
triangle (if UPLO='U') of A, including the diagonal, is
destroyed.
SOURCE
1357 subroutine wrap_ZHEEVX(jobz,range,uplo,n,a,vl,vu,il,iu,abstol,m,w,z,ldz,comm) 1358 1359 !Arguments ------------------------------------ 1360 !scalars 1361 integer,intent(in) :: il,iu,ldz,n 1362 integer,optional,intent(in) :: comm 1363 integer,intent(inout) :: m 1364 real(dp),intent(in) :: abstol,vl,vu 1365 character(len=*),intent(in) :: jobz,range,uplo 1366 !arrays 1367 real(dp),intent(out) :: w(n) 1368 complex(dpc),intent(out) :: z(ldz,m) 1369 complex(dpc),intent(inout) :: a(n,n) 1370 1371 !Local variables ------------------------------ 1372 !scalars 1373 integer :: lwork,info,nprocs 1374 logical :: use_scalapack 1375 character(len=500) :: msg 1376 !arrays 1377 integer,allocatable :: ifail(:),iwork(:) 1378 real(dp),allocatable :: rwork(:) 1379 complex(dpc),allocatable :: work(:) 1380 #ifdef HAVE_LINALG_SCALAPACK 1381 integer :: ierr,istwf_k 1382 logical :: want_eigenvectors 1383 type(matrix_scalapack) :: Slk_mat,Slk_vec 1384 type(processor_scalapack) :: Slk_processor 1385 #endif 1386 1387 !************************************************************************ 1388 1389 use_scalapack=.FALSE. 1390 if (PRESENT(comm)) then 1391 nprocs = xmpi_comm_size(comm) 1392 #ifdef HAVE_LINALG_SCALAPACK 1393 use_scalapack = (nprocs>1) 1394 #endif 1395 end if 1396 1397 SELECT CASE(use_scalapack) 1398 CASE (.FALSE.) ! Standard LAPACK call. 1399 1400 lwork = MAX(1,2*n) 1401 ABI_MALLOC(work,(lwork)) 1402 ABI_MALLOC(rwork,(7*n)) 1403 ABI_MALLOC(iwork,(5*n)) 1404 ABI_MALLOC(ifail,(n)) 1405 1406 call ZHEEVX(jobz,range,uplo,n,a,n,vl,vu,il,iu,abstol,m,w,z,ldz,work,lwork,rwork,iwork,ifail,info) 1407 1408 if (info < 0) then 1409 write(msg,'(a,i0,a)')" The ",-info,"-th argument of ZHEEVX had an illegal value." 1410 ABI_ERROR(msg) 1411 end if 1412 1413 if (info > 0) then 1414 write(msg,'(2a,i0,a)')"ZHEEVX: the algorithm failed to converge; ",ch10,& 1415 info,"eigenvectors failed to converge. " 1416 ABI_ERROR(msg) 1417 end if 1418 1419 ABI_FREE(iwork) 1420 ABI_FREE(ifail) 1421 ABI_FREE(rwork) 1422 ABI_FREE(work) 1423 RETURN 1424 1425 CASE (.TRUE.) 1426 1427 #ifdef HAVE_LINALG_SCALAPACK 1428 call Slk_processor%init(comm) 1429 istwf_k=1 1430 1431 ! Initialize and fill Scalapack matrix from the global one. 1432 call Slk_mat%init(n,n,Slk_processor,istwf_k) 1433 call slk_matrix_from_global_dpc_2D(Slk_mat,uplo,a) 1434 1435 want_eigenvectors = firstchar(jobz,(/"V","v"/)) 1436 if (want_eigenvectors) then 1437 ! Initialize the distributed vectors. 1438 call Slk_vec%init(n,n,Slk_processor,istwf_k) 1439 end if 1440 1441 ! Solve the problem. 1442 call slk_mat%pzheevx(jobz,range,uplo,vl,vu,il,iu,abstol,Slk_vec,m,w) 1443 call Slk_mat%free() 1444 1445 if (want_eigenvectors) then ! A is overwritten with the eigenvectors 1446 z = czero 1447 call slk_matrix_to_global_dpc_2D(Slk_vec,"All",z) ! Fill the entries calculated by this node. 1448 call Slk_vec%free() 1449 call xmpi_sum(z,comm,ierr) ! Fill the remaing entries of the global matrix 1450 end if 1451 1452 call Slk_processor%free() 1453 RETURN 1454 #endif 1455 1456 ABI_BUG("You should not be here!") 1457 END SELECT 1458 1459 end subroutine wrap_ZHEEVX
m_hide_lapack/wrap_ZHEGV [ Functions ]
[ Top ] [ m_hide_lapack ] [ Functions ]
NAME
wrap_ZHEGV
FUNCTION
wrap_ZHEGV computes all the eigenvalues, and optionally, the eigenvectors of a complex generalized
Hermitian-definite eigenproblem, of the form
A*x=(lambda)*B*x (1),
A*Bx=(lambda)*x, (2), or
B*A*x=(lambda)*x (3).
Here A and B are assumed to be Hermitian and B is also positive definite.
INPUTS
ITYPE (input) INTEGER Specifies the problem type to be solved:
= 1: A*x = (lambda)*B*x
= 2: A*B*x = (lambda)*x
= 3: B*A*x = (lambda)*x
JOBZ (input) CHARACTER*1
= "N": Compute eigenvalues only;
= "V": Compute eigenvalues and eigenvectors.
UPLO (input) CHARACTER*1
= "U": Upper triangle of A and B are stored;
= "L": Lower triangle of A and B are stored.
N (input) INTEGER
The order of the matrices A and B. N >= 0.
[comm]=MPI communicator for ScaLAPACK inversion. Only available if the code has been compiled with Scalapack support.
To avoid wasting CPU time the scalapack initialization is avoided if the number of processors in 1,
in this case the sequential LAPACK routine is called.
OUTPUT
W (output) REAL(DP) array, dimension (N)
If INFO = 0, the eigenvalues in ascending order.
See also SIDE EFFECTS
SIDE EFFECTS
A (input/output) COMPLEX(DPC) array, dimension (N, N)
On entry, the Hermitian matrix A. If UPLO = "U", the leading N-by-N upper triangular part of A
<S-F1>contains the upper triangular part of the matrix A.
If UPLO = "L", the leading N-by-N lower triangular part of A contains the lower triangular part of the matrix A.
On exit, if JOBZ = "V", then A contains the matrix Z of eigenvectors.
The eigenvectors are normalized as follows: if ITYPE = 1 or 2,
Z**H*B*Z = I; if ITYPE = 3, Z**H*inv(B)*Z = I.
If JOBZ = "N", then on exit the upper triangle (if UPLO="U") or the lower triangle
(if UPLO="L") of A, including the diagonal, is destroyed.
B (input/output) COMPLEX*16 array, dimension (LDB, N)
On entry, the Hermitian positive definite matrix B.
If UPLO = "U", the leading N-by-N upper triangular part of B contains the upper triangular part of the matrix B.
If UPLO = "L", the leading N-by-N lower triangular part of B contains the lower triangular part of the matrix B.
On exit, if INFO <= N, the part of B containing the matrix is overwritten by the triangular
factor U or L from the Cholesky factorization B = U**H*U or B = L*L**H.
SOURCE
917 subroutine wrap_ZHEGV(itype, jobz, uplo, n, a, b, w, comm) 918 919 !Arguments ------------------------------------ 920 !scalars 921 integer,intent(in) :: n,itype 922 integer,optional,intent(in) :: comm 923 character(len=*),intent(in) :: jobz,uplo 924 !arrays 925 complex(dpc),intent(inout) :: a(n,n),b(n,n) 926 real(dp),intent(out) :: w(n) 927 928 !Local variables ------------------------------ 929 !scalars 930 integer :: lwork,info,nprocs,ii 931 logical :: use_scalapack 932 character(len=500) :: msg 933 !arrays 934 real(dp),allocatable :: rwork(:) 935 complex(dpc),allocatable :: work(:) 936 #ifdef HAVE_LINALG_SCALAPACK 937 integer :: ierr,istwf_k 938 type(matrix_scalapack) :: Slk_matA,Slk_matB 939 type(processor_scalapack) :: Slk_processor 940 #endif 941 !************************************************************************ 942 943 use_scalapack=.FALSE. 944 if (PRESENT(comm)) then 945 nprocs = xmpi_comm_size(comm) 946 #ifdef HAVE_LINALG_SCALAPACK 947 use_scalapack = (nprocs>1) 948 #endif 949 end if 950 951 SELECT CASE(use_scalapack) 952 CASE (.FALSE.) 953 lwork = MAX(1,2*n-1) 954 955 ABI_MALLOC(work, (lwork)) 956 ABI_MALLOC(rwork, (MAX(1,3*n-2))) 957 958 call ZHEGV(itype,jobz,uplo,n,a,n,b,n,w,work,lwork,rwork,info) 959 960 if (info < 0) then 961 write(msg,'(a,i0,a)')" The ",-info,"-th argument of ZHEGV had an illegal value." 962 ABI_ERROR(msg) 963 end if 964 965 if (info > 0) then 966 if (info<= n) then 967 write(msg,'(2a,i0,a)')& 968 "ZHEGV failed to converge: ",ch10,& 969 info," off-diagonal elements of an intermediate tridiagonal form did not converge to zero. " 970 else 971 ii = info -n 972 write(msg,'(3a,i0,3a)')& 973 "ZHEGV failed to converge: ",ch10,& 974 "The leading minor of order ",ii," of B is not positive definite. ",ch10,& 975 "The factorization of B could not be completed and no eigenvalues or eigenvectors were computed." 976 end if 977 ABI_ERROR(msg) 978 end if 979 980 ABI_FREE(rwork) 981 ABI_FREE(work) 982 RETURN 983 984 CASE (.TRUE.) 985 986 #ifdef HAVE_LINALG_SCALAPACK 987 call Slk_processor%init(comm) 988 istwf_k=1 989 990 ! Initialize and fill Scalapack matrix from the global one. 991 call Slk_matA%init(n,n,Slk_processor,istwf_k) 992 call slk_matrix_from_global_dpc_2D(Slk_matA,uplo,a) 993 994 call Slk_matB%init(n,n,Slk_processor,istwf_k) 995 call slk_matrix_from_global_dpc_2D(Slk_matB,uplo,b) 996 997 ! Solve the problem with scaLAPACK. 998 ABI_ERROR("slk_pZHEGV not yet coded") 999 ! TODO 1000 !% call slk_pzhegv(itype,jobz,uplo,Slk_matA,Slk_matB,w) 1001 1002 call Slk_matB%free() 1003 1004 if (firstchar(jobz,(/"V","v"/))) then ! A is overwritten with the eigenvectors 1005 a = czero 1006 call slk_matrix_to_global_dpc_2D(Slk_matA,"All",a) ! Fill the entries calculated by this node. 1007 call xmpi_sum(a,comm,ierr) ! Fill the remaing entries of the global matrix 1008 end if 1009 1010 call Slk_matA%free() 1011 call Slk_processor%free() 1012 RETURN 1013 #endif 1014 1015 ABI_BUG("You should not be here!") 1016 END SELECT 1017 1018 end subroutine wrap_ZHEGV
m_hide_lapack/wrap_ZHEGVX [ Functions ]
[ Top ] [ m_hide_lapack ] [ Functions ]
NAME
wrap_ZHEGVX
FUNCTION
wrap_ZHEGVX - compute selected eigenvalues, and optionally, eigenvectors of a complex generalized Hermitian-definite eigenproblem, of the form A*x=(lambda)*B*x, A*Bx=(lambda)*x, or B*A*x=(lambda)*x. Eigenvalues and eigenvectors can be selected by specifying either a range of values or a range of indices for the desired eigenvalues.
INPUTS
ITYPE (input) INTEGER Specifies the problem type to be solved:
= 1: A*x = (lambda)*B*x
= 2: A*B*x = (lambda)*x
= 3: B*A*x = (lambda)*x
JOBZ (input) CHARACTER*1
= 'N': Compute eigenvalues only;
= 'V': Compute eigenvalues and eigenvectors.
RANGE (input) CHARACTER*1
= 'A': all eigenvalues will be found.
= 'V': all eigenvalues in the half-open interval (VL,VU]
will be found.
= 'I': the IL-th through IU-th eigenvalues will be found.
UPLO (input) CHARACTER*1
= 'U': Upper triangle of A is stored;
= 'L': Lower triangle of A is stored.
N (input) INTEGER
The order of the matrices A and B. N >= 0.
LDA (input) INTEGER
The leading dimension of the array A. LDA >= max(1,N).
VL (input) REAL(DP)
VU (input) REAL(DP)
If RANGE='V', the lower and upper bounds of the interval to
be searched for eigenvalues. VL < VU.
Not referenced if RANGE = 'A' or 'I'.
IL (input) INTEGER
IU (input) INTEGER
If RANGE='I', the indices (in ascending order) of the
smallest and largest eigenvalues to be returned.
1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
Not referenced if RANGE = 'A' or 'V'.
ABSTOL (input) REAL(DP)
The absolute error tolerance for the eigenvalues.
An approximate eigenvalue is accepted as converged
when it is determined to lie in an interval [a,b]
of width less than or equal to
ABSTOL + EPS * max( |a|,|b| ) ,
where EPS is the machine precision. If ABSTOL is less than
or equal to zero, then EPS*|T| will be used in its place,
where |T| is the 1-norm of the tridiagonal matrix obtained
by reducing A to tridiagonal form.
Eigenvalues will be computed most accurately when ABSTOL is
set to twice the underflow threshold 2*DLAMCH('S'), not zero.
If this routine returns with INFO>0, indicating that some
eigenvectors did not converge, try setting ABSTOL to
2*DLAMCH('S').
LDZ (input) INTEGER
The leading dimension of the array Z. LDZ >= 1, and if
JOBZ = 'V', LDZ >= max(1,N).
[comm]=MPI communicator for ScaLAPACK inversion. Only available if the code has been compiled with Scalapack support.
To avoid wasting CPU time the scalapack initialization is avoided if the number of processors in 1,
in this case the sequential LAPACK routine is called.
OUTPUT
M (output) INTEGER
The total number of eigenvalues found. 0 <= M <= N.
If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
W (output) REAL(DP) array, dimension (N)
On normal exit, the first M elements contain the selected
eigenvalues in ascending order.
Z (output) COMPLEX(DPC) array, dimension (LDZ, max(1,M))
If JOBZ = 'V', then if INFO = 0, the first M columns of Z
contain the orthonormal eigenvectors of the matrix A
corresponding to the selected eigenvalues, with the i-th
column of Z holding the eigenvector associated with W(i).
The eigenvectors are normalized as follows: if ITYPE = 1 or 2, Z**T*B*Z = I; if ITYPE = 3, Z**T*inv(B)*Z = I.
If an eigenvector fails to converge, then that column of Z
contains the latest approximation to the eigenvector, and the
index of the eigenvector is returned in IFAIL.
If JOBZ = 'N', then Z is not referenced.
Note: the user must ensure that at least max(1,M) columns are
supplied in the array Z; if RANGE = 'V', the exact value of M
is not known in advance and an upper bound must be used.
See also SIDE EFFECTS
SIDE EFFECTS
A (input/output) COMPLEX(DPC) array, dimension (N, N)
On entry, the Hermitian matrix A. If UPLO = 'U', the
leading N-by-N upper triangular part of A contains the
upper triangular part of the matrix A. If UPLO = "L",
the leading N-by-N lower triangular part of A contains
the lower triangular part of the matrix A.
On exit, the lower triangle (if UPLO="L") or the upper
triangle (if UPLO="U") of A, including the diagonal, is
destroyed.
B (input/output) COMPLEX(DPC) array, dimension (LDB, N)
On entry, the Hermitian matrix B. If UPLO = "U", the leading N-by-N upper triangular part
of B contains the upper triangular part of the matrix B.
If UPLO = "L", the leading N-by-N lower triangular part of B contains the lower triangular part of the matrix B.
On exit, if INFO <= N, the part of B containing the matrix is overwritten by the triangular factor
U or L from the Cholesky factorization B = U**H*U or B = L*L**H.
SOURCE
1846 subroutine wrap_ZHEGVX(itype,jobz,range,uplo,n,a,b,vl,vu,il,iu,abstol,m,w,z,ldz,comm) 1847 1848 !Arguments ------------------------------------ 1849 !scalars 1850 integer,intent(in) :: il,iu,ldz,n,itype 1851 integer,optional,intent(in) :: comm 1852 integer,intent(inout) :: m 1853 real(dp),intent(in) :: abstol,vl,vu 1854 character(len=*),intent(in) :: jobz,range,uplo 1855 !arrays 1856 real(dp),intent(out) :: w(n) 1857 !complex(dpc),intent(out) :: z(ldz,n) 1858 complex(dpc),intent(out) :: z(ldz,m) 1859 complex(dpc),intent(inout) :: a(n,n),b(n,n) 1860 1861 !Local variables ------------------------------ 1862 !scalars 1863 integer :: lwork,info,nprocs,ii 1864 logical :: use_scalapack 1865 character(len=500) :: msg 1866 !arrays 1867 integer,allocatable :: ifail(:),iwork(:) 1868 real(dp),allocatable :: rwork(:) 1869 complex(dpc),allocatable :: work(:) 1870 #ifdef HAVE_LINALG_SCALAPACK 1871 integer :: ierr,istwf_k 1872 logical :: want_eigenvectors 1873 type(matrix_scalapack) :: Slk_matA,Slk_matB,Slk_vec 1874 type(processor_scalapack) :: Slk_processor 1875 #endif 1876 1877 !************************************************************************ 1878 1879 use_scalapack=.FALSE. 1880 if (PRESENT(comm)) then 1881 nprocs = xmpi_comm_size(comm) 1882 #ifdef HAVE_LINALG_SCALAPACK 1883 use_scalapack = (nprocs>1) 1884 #endif 1885 end if 1886 1887 SELECT CASE(use_scalapack) 1888 CASE (.FALSE.) 1889 ! Standard LAPACK call. 1890 lwork = MAX(1,2*n) 1891 ABI_MALLOC(work,(lwork)) 1892 ABI_MALLOC(rwork,(7*n)) 1893 ABI_MALLOC(iwork,(5*n)) 1894 ABI_MALLOC(ifail,(n)) 1895 1896 call ZHEGVX(itype,jobz,range,uplo,n,a,n,b,n,vl,vu,il,iu,abstol,m,w,z,ldz,work,lwork,rwork,iwork,ifail,info) 1897 1898 if (info < 0) then 1899 write(msg,'(a,i0,a)')" The ",-info,"-th argument of ZHEGVX had an illegal value." 1900 ABI_ERROR(msg) 1901 end if 1902 1903 if (info > 0) then 1904 if (info<= n) then 1905 write(msg,'(a,i0,a)')"ZHEGVX failed to converge: ",info," eigenvectors failed to converge. " 1906 else 1907 ii = info -n 1908 write(msg,'(3a,i0,3a)')& 1909 "ZHEEVX failed to converge: ",ch10,& 1910 "The leading minor of order ",ii," of B is not positive definite. ",ch10,& 1911 "The factorization of B could not be completed and no eigenvalues or eigenvectors were computed." 1912 end if 1913 ABI_ERROR(msg) 1914 end if 1915 1916 ABI_FREE(iwork) 1917 ABI_FREE(ifail) 1918 ABI_FREE(rwork) 1919 ABI_FREE(work) 1920 RETURN 1921 1922 CASE (.TRUE.) 1923 1924 #ifdef HAVE_LINALG_SCALAPACK 1925 call Slk_processor%init(comm) 1926 istwf_k=1 1927 1928 ! Initialize and fill Scalapack matrix from the global one. 1929 call Slk_matA%init(n,n,Slk_processor,istwf_k) 1930 call slk_matrix_from_global_dpc_2D(Slk_matA,uplo,a) 1931 1932 call Slk_matB%init(n,n,Slk_processor,istwf_k) 1933 call slk_matrix_from_global_dpc_2D(Slk_matB,uplo,b) 1934 1935 want_eigenvectors = firstchar(jobz,(/"V","v"/)) 1936 if (want_eigenvectors) then ! Initialize the distributed vectors. 1937 call Slk_vec%init(n,n,Slk_processor,istwf_k) 1938 end if 1939 1940 ! Solve the problem. 1941 ABI_ERROR("slk_pZHEGVX not coded yet") 1942 ! TODO write the scaLAPACK wrapper. 1943 !call slk_pZHEGVX(itype,jobz,range,uplo,Slk_matA,Slk_matB,vl,vu,il,iu,abstol,Slk_vec,m,w) 1944 1945 call Slk_matA%free() 1946 call Slk_matB%free() 1947 1948 if (want_eigenvectors) then ! A is overwritten with the eigenvectors 1949 z = czero 1950 call slk_matrix_to_global_dpc_2D(Slk_vec,"All",z) ! Fill the entries calculated by this node. 1951 call Slk_vec%free() 1952 call xmpi_sum(z,comm,ierr) ! Fill the remaing entries of the global matrix 1953 end if 1954 1955 call Slk_processor%free() 1956 1957 RETURN 1958 #endif 1959 1960 ABI_BUG("You should not be here!") 1961 END SELECT 1962 1963 end subroutine wrap_ZHEGVX
m_hide_lapack/wrap_ZHPEV [ Functions ]
[ Top ] [ m_hide_lapack ] [ Functions ]
NAME
wrap_ZHPEV
FUNCTION
wrap_ZHPEV computes all the eigenvalues and, optionally, eigenvectors of a complex Hermitian matrix in packed storage. Scalapack version is not available. [PRIVATE].
INPUTS
JOBZ (input) CHARACTER*1
= 'N': Compute eigenvalues only;
= 'V': Compute eigenvalues and eigenvectors.
UPLO (input) CHARACTER*1
= 'U': Upper triangle of A is stored;
= 'L': Lower triangle of A is stored.
N (input) INTEGER
The order of the matrix A. N >= 0.
LDZ (input) INTEGER
The leading dimension of the array Z. LDZ >= 1, and if
JOBZ = 'V', LDZ >= max(1,N).
[comm]=MPI communicator for ScaLAPACK inversion. Only available if the code has been compiled with Scalapack support.
To avoid wasting CPU time the scalapack initialization is avoided if the number of processors in 1,
in this case the sequential LAPACK routine is called. Note that scalapack does not provide native
support for packed symmetric matrices. Threfore we have to distribute the full matrix among the nodes.
in order to perform the calculation in parallel.
OUTPUT
W (output) REAL(DP) array, dimension (N)
If INFO = 0, the eigenvalues in ascending order.
Z (output) COMPLEX(DPC) array, dimension (LDZ, N)
If JOBZ = 'V', then if INFO = 0, Z contains the orthonormal
eigenvectors of the matrix A, with the i-th column of Z
holding the eigenvector associated with W(i).
If JOBZ = 'N', then Z is not referenced.
See also SIDE EFFECTS
SIDE EFFECTS
AP (input/output) COMPLEX(DPC) array, dimension (N*(N+1)/2)
On entry, the upper or lower triangle of the Hermitian matrix
A, packed columnwise in a linear array. The j-th column of A
is stored in the array AP as follows:
if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n.
On exit, AP is overwritten by values generated during the
reduction to tridiagonal form. If UPLO = 'U', the diagonal
and first superdiagonal of the tridiagonal matrix T overwrite
the corresponding elements of A, and if UPLO = 'L', the
diagonal and first subdiagonal of T overwrite the
corresponding elements of A. Unchanged if ScaLAPACK is used.
SOURCE
754 subroutine wrap_ZHPEV(jobz, uplo, n, ap, w, z, ldz, comm) 755 756 !Arguments ------------------------------------ 757 !scalars 758 integer,intent(in) :: n,ldz 759 integer,optional,intent(in) :: comm 760 character(len=*),intent(in) :: jobz,uplo 761 !arrays 762 real(dp),intent(out) :: w(n) 763 complex(dpc),intent(inout) :: ap(n*(n+1)/2) 764 complex(dpc),intent(out) :: z(ldz,n) 765 766 !Local variables ------------------------------ 767 !scalars 768 integer :: info,nprocs 769 logical :: use_scalapack 770 character(len=500) :: msg 771 !arrays 772 real(dp),allocatable :: rwork(:) 773 complex(dpc),allocatable :: work(:) 774 #ifdef HAVE_LINALG_SCALAPACK 775 integer :: ierr,istwf_k 776 logical :: want_eigenvectors 777 type(matrix_scalapack) :: Slk_mat,Slk_vec 778 type(processor_scalapack) :: Slk_processor 779 #endif 780 781 !************************************************************************ 782 783 use_scalapack=.FALSE. 784 if (PRESENT(comm)) then 785 nprocs = xmpi_comm_size(comm) 786 #ifdef HAVE_LINALG_SCALAPACK 787 use_scalapack = (nprocs>1) 788 #endif 789 end if 790 791 SELECT CASE(use_scalapack) 792 CASE (.FALSE.) 793 ABI_MALLOC(work, (MAX(1,2*n-1))) 794 ABI_MALLOC(rwork, (MAX(1,3*n-2))) 795 796 call ZHPEV(jobz,uplo,n,ap,w,z,ldz,work,rwork,info) 797 798 if (info < 0) then 799 write(msg,'(a,i0,a)')" The ",-info,"-th argument of ZHPEV had an illegal value." 800 ABI_ERROR(msg) 801 end if 802 803 if (info > 0) then 804 write(msg,'(2a,i0,a)')& 805 "ZHPEV: the algorithm failed to converge; ",ch10,& 806 info," off-diagonal elements of an intermediate tridiagonal form did not converge to zero. " 807 ABI_ERROR(msg) 808 end if 809 810 ABI_FREE(rwork) 811 ABI_FREE(work) 812 RETURN 813 814 CASE (.TRUE.) 815 816 #ifdef HAVE_LINALG_SCALAPACK 817 call Slk_processor%init(comm) 818 istwf_k=1 819 820 ! Initialize and fill Scalapack matrix from the global one. 821 call Slk_mat%init(n,n,Slk_processor,istwf_k) 822 call slk_matrix_from_global_dpc_1Dp(Slk_mat,uplo,ap) 823 824 want_eigenvectors = firstchar(jobz,(/"V","v"/)) 825 if (want_eigenvectors) then 826 ! Initialize the distributed vectors. 827 call Slk_vec%init(n,n,Slk_processor,istwf_k) 828 end if 829 830 ! Solve the problem with scaLAPACK. 831 call slk_mat%heev(jobz,uplo,Slk_vec,w) 832 call Slk_mat%free() 833 834 if (want_eigenvectors) then ! Collect the eigenvectors. 835 z = zero 836 call slk_matrix_to_global_dpc_2D(Slk_vec,"All",z) ! Fill the entries calculated by this node. 837 call Slk_vec%free() 838 call xmpi_sum(z,comm,ierr) ! Fill the remaing entries of the global matrix 839 end if 840 841 call Slk_processor%free() 842 843 RETURN 844 #endif 845 846 ABI_BUG("You should not be here!") 847 END SELECT 848 849 end subroutine wrap_ZHPEV
m_hide_lapack/xheev_cplex [ Functions ]
[ Top ] [ m_hide_lapack ] [ Functions ]
NAME
xheev_cplex
FUNCTION
xheev_cplex computes the eigenvalues and, optionally, the eigenvectors of a (complex Hermitian| real symmetric) matrix in double precision.
INPUTS
JOBZ (input) CHARACTER*1
= 'N': Compute eigenvalues only;
= 'V': Compute eigenvalues and eigenvectors.
UPLO (input) CHARACTER*1
= 'U': Upper triangle of A is stored;
= 'L': Lower triangle of A is stored.
CPLEX= Size of the first dimension of the A matrix.
1 for a real symmetric matrix.
2 for complex Hermitian matrix stored in a real array with real and imaginary part.
N (input) INTEGER
The order of the matrix A. N >= 0.
[comm]=MPI communicator for ScaLAPACK inversion. Only available if the code has been compiled with Scalapack support.
To avoid wasting CPU time the scalapack initialization is avoided if the number of processors in 1,
in this case the sequential LAPACK routine is called.
OUTPUT
W (output) REAL(DP) array, dimension (N)
If INFO = 0, the eigenvalues in ascending order.
See also SIDE EFFECTS
SIDE EFFECTS
A (input/output) REAL(DP) array, dimension (CPLEX, N, N)
On entry, the (complex Hermitian|Real symmetric) matrix A. If UPLO = 'U', the
leading N-by-N upper triangular part of A contains the
upper triangular part of the matrix A. If UPLO = 'L',
the leading N-by-N lower triangular part of A contains
the lower triangular part of the matrix A.
On exit, if JOBZ = 'V', then if INFO = 0, A contains the
orthonormal eigenvectors of the matrix A.
If JOBZ = 'N', then on exit the lower triangle (if UPLO='L')
or the upper triangle (if UPLO='U') of A, including the
diagonal, is destroyed.
SOURCE
464 subroutine xheev_cplex(jobz, uplo, cplex, n, a, w, msg, ierr, comm) 465 466 !Arguments ------------------------------------ 467 !scalars 468 integer,intent(in) :: n,cplex 469 integer,optional,intent(in) :: comm 470 character(len=*),intent(in) :: jobz,uplo 471 integer,intent(out) :: ierr 472 character(len=*),intent(out) :: msg 473 !arrays 474 real(dp),intent(inout) :: a(cplex,n,n) 475 real(dp),intent(out) :: w(n) 476 477 !Local variables ------------------------------ 478 !scalars 479 integer :: lwork,nprocs 480 logical :: use_scalapack 481 !arrays 482 real(dp),allocatable :: rwork(:) 483 real(dp),allocatable :: work_real(:) 484 complex(dpc),allocatable :: work_cplx(:) 485 #ifdef HAVE_LINALG_SCALAPACK 486 !integer :: istwf_k 487 !logical :: want_eigenvectors 488 !type(matrix_scalapack) :: Slk_mat,Slk_vec 489 !type(processor_scalapack) :: Slk_processor 490 #endif 491 !************************************************************************ 492 493 use_scalapack=.FALSE. 494 if (PRESENT(comm)) then 495 nprocs = xmpi_comm_size(comm) 496 #ifdef HAVE_LINALG_SCALAPACK 497 use_scalapack = (nprocs>1) 498 #endif 499 end if 500 501 if (ALL(cplex/= [1, 2])) then 502 write(msg,'(a,i0)')" Wrong value for cplex: ",cplex 503 ierr = 1; return 504 end if 505 506 SELECT CASE(use_scalapack) 507 CASE (.FALSE.) 508 if (cplex==1) then 509 ! Real symmetric case. 510 lwork = MAX(1,3*n-1) 511 ABI_MALLOC(work_real,(lwork)) 512 513 call DSYEV(jobz,uplo,n,a,n,w,work_real,lwork,ierr) 514 515 if (ierr < 0) then 516 write(msg,'(a,i0,a)')" The ",-ierr,"-th argument of DSYEV had an illegal value." 517 end if 518 519 if (ierr > 0) then 520 write(msg,'(2a,i0,a)')& 521 "DSYEV: the algorithm failed to converge; ",ch10,& 522 ierr," off-diagonal elements of an intermediate tridiagonal form did not converge to zero. " 523 end if 524 525 ABI_FREE(work_real) 526 RETURN 527 528 else 529 ! Hermitian case. 530 lwork = MAX(1,2*n-1) 531 532 ABI_MALLOC(work_cplx, (lwork)) 533 ABI_MALLOC(rwork, (MAX(1,3*n-2))) 534 535 call ZHEEV(jobz,uplo,n,a,n,w,work_cplx,lwork,rwork,ierr) 536 537 if (ierr < 0) then 538 write(msg,'(a,i0,a)')" The ",-ierr,"-th argument of ZHEEV had an illegal value." 539 end if 540 541 if (ierr > 0) then 542 write(msg,'(2a,i0,a)')& 543 "ZHEEV: the algorithm failed to converge; ",ch10,& 544 ierr," off-diagonal elements of an intermediate tridiagonal form did not converge to zero. " 545 end if 546 547 ABI_FREE(rwork) 548 ABI_FREE(work_cplx) 549 RETURN 550 end if ! cplex 551 552 CASE (.TRUE.) 553 554 #ifdef HAVE_LINALG_SCALAPACK 555 ABI_ERROR("Not coded yet") 556 557 !call Slk_processor%init(comm) 558 !istwf_k=1 559 ! 560 !! Initialize and fill Scalapack matrix from the global one. 561 !call Slk_mat%init(n,n,Slk_processor,istwf_k) 562 ! 563 !call slk_matrix_from_global_dpc_2D(Slk_mat,uplo,a) 564 ! 565 !want_eigenvectors = firstchar(jobz,(/"V","v"/)) 566 !if (want_eigenvectors) then ! Initialize the distributed vectors. 567 ! call Slk_vec%init(n,n,Slk_processor,istwf_k) 568 !end if 569 ! 570 !! Solve the problem with scaLAPACK. 571 !call slk_mat%heev(jobz,uplo,Slk_vec,w) 572 !call Slk_mat%free() 573 ! 574 !if (want_eigenvectors) then ! A is overwritten with the eigenvectors 575 ! a = czero 576 ! call slk_matrix_to_global_dpc_2D(Slk_vec,"All",a) ! Fill the entries calculated by this node. 577 ! call Slk_vec%free() 578 ! call xmpi_sum(a,comm,ierr) ! Fill the remaing entries of the global matrix 579 !end if 580 ! 581 !call Slk_processor%free() 582 583 RETURN 584 #endif 585 586 ABI_BUG("You should not be here!") 587 END SELECT 588 589 end subroutine xheev_cplex
m_hide_lapack/xheevx_cplex [ Functions ]
[ Top ] [ m_hide_lapack ] [ Functions ]
NAME
xheevx_cplex
FUNCTION
xheevx_cplex computes selected eigenvalues and, optionally, eigenvectors of a (real symmetric|complex Hermitian) matrix A. Eigenvalues and eigenvectors can be selected by specifying either a range of values or a range of indices for the desired eigenvalues.
INPUTS
JOBZ (input) CHARACTER*1
= 'N': Compute eigenvalues only;
= 'V': Compute eigenvalues and eigenvectors.
RANGE (input) CHARACTER*1
= 'A': all eigenvalues will be found.
= 'V': all eigenvalues in the half-open interval (VL,VU]
will be found.
= 'I': the IL-th through IU-th eigenvalues will be found.
UPLO (input) CHARACTER*1
= 'U': Upper triangle of A is stored;
= 'L': Lower triangle of A is stored.
CPLEX Size of the first dimension of the matrix A.
1 for real symmetric matrix
2 for complex Hermitian matrix.
N (input) INTEGER
The order of the matrix A. N >= 0.
LDA (input) INTEGER
The leading dimension of the array A. LDA >= max(1,N).
VL (input) REAL(DP)
VU (input) REAL(DP)
If RANGE='V', the lower and upper bounds of the interval to
be searched for eigenvalues. VL < VU.
Not referenced if RANGE = 'A' or 'I'.
IL (input) INTEGER
IU (input) INTEGER
If RANGE='I', the indices (in ascending order) of the
smallest and largest eigenvalues to be returned.
1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
Not referenced if RANGE = 'A' or 'V'.
ABSTOL (input) REAL(DP)
The absolute error tolerance for the eigenvalues.
An approximate eigenvalue is accepted as converged
when it is determined to lie in an interval [a,b]
of width less than or equal to
ABSTOL + EPS * max( |a|,|b| ) ,
where EPS is the machine precision. If ABSTOL is less than
or equal to zero, then EPS*|T| will be used in its place,
where |T| is the 1-norm of the tridiagonal matrix obtained
by reducing A to tridiagonal form.
Eigenvalues will be computed most accurately when ABSTOL is
set to twice the underflow threshold 2*DLAMCH('S'), not zero.
If this routine returns with INFO>0, indicating that some
eigenvectors did not converge, try setting ABSTOL to
2*DLAMCH('S').
See "Computing Small Singular Values of Bidiagonal Matrices
with Guaranteed High Relative Accuracy," by Demmel and
Kahan, LAPACK Working Note #3.
LDZ (input) INTEGER
The leading dimension of the array Z. LDZ >= 1, and if
JOBZ = 'V', LDZ >= max(1,N).
[comm]=MPI communicator for ScaLAPACK inversion. Only available if the code has been compiled with Scalapack support.
To avoid wasting CPU time the scalapack initialization is avoided if the number of processors in 1,
in this case the sequential LAPACK routine is called.
OUTPUT
M (output) INTEGER
The total number of eigenvalues found. 0 <= M <= N.
If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
W (output) REAL(DP) array, dimension (N)
On normal exit, the first M elements contain the selected
eigenvalues in ascending order.
Z (output) REAL(DP) array, dimension (CPLEX, LDZ, max(1,M))
If JOBZ = 'V', then if INFO = 0, the first M columns of Z
contain the orthonormal eigenvectors of the matrix A
corresponding to the selected eigenvalues, with the i-th
column of Z holding the eigenvector associated with W(i).
If an eigenvector fails to converge, then that column of Z
contains the latest approximation to the eigenvector, and the
index of the eigenvector is returned in IFAIL.
If JOBZ = 'N', then Z is not referenced.
Note: the user must ensure that at least max(1,M) columns are
supplied in the array Z; if RANGE = 'V', the exact value of M
is not known in advance and an upper bound must be used.
See also SIDE EFFECTS
SIDE EFFECTS
A (input/output) REAL(DP) array, dimension (CPLEX, N, N)
On entry, the (real symmetric|complex Hermitian) matrix A. If UPLO = 'U', the
leading N-by-N upper triangular part of A contains the
upper triangular part of the matrix A. If UPLO = 'L',
the leading N-by-N lower triangular part of A contains
the lower triangular part of the matrix A.
On exit, the lower triangle (if UPLO='L') or the upper
triangle (if UPLO='U') of A, including the diagonal, is
destroyed.
SOURCE
1580 subroutine xheevx_cplex(jobz, range, uplo, cplex, n, a, vl, vu, il, iu, & 1581 abstol, m, w, z, ldz, msg, ierr, comm) 1582 1583 !Arguments ------------------------------------ 1584 !scalars 1585 integer,intent(in) :: il,iu,ldz,n,cplex 1586 integer,optional,intent(in) :: comm 1587 integer,intent(inout) :: m 1588 integer,intent(out) :: ierr 1589 real(dp),intent(in) :: abstol,vl,vu 1590 character(len=*),intent(in) :: jobz,range,uplo 1591 character(len=*),intent(out) :: msg 1592 !arrays 1593 real(dp),intent(out) :: w(n) 1594 !real(dp),intent(out) :: z(cplex,ldz,n) 1595 real(dp),intent(out) :: z(cplex,ldz,m) 1596 real(dp),intent(inout) :: a(cplex,n,n) 1597 1598 !Local variables ------------------------------ 1599 !scalars 1600 integer :: lwork,nprocs 1601 logical :: use_scalapack 1602 !arrays 1603 integer,allocatable :: ifail(:),iwork(:) 1604 real(dp),allocatable :: rwork(:) 1605 real(dp),allocatable :: work_real(:) 1606 complex(dpc),allocatable :: work_cplx(:) 1607 #ifdef HAVE_LINALG_SCALAPACK 1608 !integer :: istwf_k 1609 !logical :: want_eigenvectors 1610 !type(matrix_scalapack) :: Slk_mat,Slk_vec 1611 !type(processor_scalapack) :: Slk_processor 1612 #endif 1613 1614 !************************************************************************ 1615 1616 use_scalapack=.FALSE. 1617 if (PRESENT(comm)) then 1618 nprocs = xmpi_comm_size(comm) 1619 #ifdef HAVE_LINALG_SCALAPACK 1620 use_scalapack = (nprocs>1) 1621 #endif 1622 end if 1623 1624 if (ALL(cplex/=(/1,2/))) then 1625 write(msg,'(a,i0)')" Wrong value for cplex: ",cplex 1626 ierr = 1; return 1627 end if 1628 1629 SELECT CASE(use_scalapack) 1630 CASE (.FALSE.) 1631 ! Standard LAPACK call. 1632 1633 if (cplex==1) then 1634 ! Real symmetric case 1635 lwork = MAX(1,8*n) 1636 ABI_MALLOC(work_real,(lwork)) 1637 ABI_MALLOC(iwork,(5*n)) 1638 ABI_MALLOC(ifail,(n)) 1639 1640 call DSYEVX(jobz,range,uplo,n,a,n,vl,vu,il,iu,abstol,m,w,z,ldz,work_real,lwork,iwork,ifail,ierr) 1641 1642 if (ierr < 0) then 1643 write(msg,'(a,i0,a)')" The ",-ierr,"-th argument of DSYEVX had an illegal value." 1644 end if 1645 1646 if (ierr > 0) then 1647 write(msg,'(2a,i0,a)')& 1648 "DSYEVX: the algorithm failed to converge; ",ch10,ierr,"eigenvectors failed to converge. " 1649 end if 1650 1651 ABI_FREE(work_real) 1652 ABI_FREE(iwork) 1653 ABI_FREE(ifail) 1654 RETURN 1655 1656 else 1657 ! Complex Hermitian case. 1658 lwork = MAX(1,2*n) 1659 ABI_MALLOC(work_cplx,(lwork)) 1660 ABI_MALLOC(rwork,(7*n)) 1661 ABI_MALLOC(iwork,(5*n)) 1662 ABI_MALLOC(ifail,(n)) 1663 1664 call ZHEEVX(jobz,range,uplo,n,a,n,vl,vu,il,iu,abstol,m,w,z,ldz,work_cplx,lwork,rwork,iwork,ifail,ierr) 1665 1666 if (ierr < 0) then 1667 write(msg,'(a,i0,a)')" The ",-ierr,"-th argument of ZHEEVX had an illegal value." 1668 end if 1669 1670 if (ierr > 0) then 1671 write(msg,'(2a,i0,a)')& 1672 "ZHEEVX: the algorithm failed to converge; ",ch10,ierr,"eigenvectors failed to converge. " 1673 end if 1674 1675 ABI_FREE(iwork) 1676 ABI_FREE(ifail) 1677 ABI_FREE(rwork) 1678 ABI_FREE(work_cplx) 1679 RETURN 1680 end if 1681 1682 CASE (.TRUE.) 1683 1684 #ifdef HAVE_LINALG_SCALAPACK 1685 ABI_ERROR("Not coded yet") 1686 ! call Slk_processor%init(comm) 1687 ! istwf_k=1 1688 1689 ! ! Initialize and fill Scalapack matrix from the global one. 1690 ! call Slk_mat%init(n,n,Slk_processor,istwf_k) 1691 ! call slk_matrix_from_global_dpc_2D(Slk_mat,uplo,a) 1692 1693 ! want_eigenvectors = firstchar(jobz,(/"V","v"/)) 1694 ! if (want_eigenvectors) then ! Initialize the distributed vectors. 1695 ! call Slk_vec%init(n,n,Slk_processor,istwf_k) 1696 ! end if 1697 1698 ! ! Solve the problem. 1699 ! call slk_mat%pzheevx(jobz,range,uplo,vl,vu,il,iu,abstol,Slk_vec,m,w) 1700 ! call Slk_mat%free() 1701 ! 1702 ! if (want_eigenvectors) then ! A is overwritten with the eigenvectors 1703 ! z = czero 1704 ! call slk_matrix_to_global_dpc_2D(Slk_vec,"All",z) ! Fill the entries calculated by this node. 1705 ! call Slk_vec%free() 1706 ! call xmpi_sum(z,comm,ierr) ! Fill the remaing entries of the global matrix 1707 ! end if 1708 1709 ! call Slk_processor%free() 1710 1711 RETURN 1712 #endif 1713 1714 ABI_BUG("You should not be here!") 1715 END SELECT 1716 1717 end subroutine xheevx_cplex
m_hide_lapack/xhegv_cplex [ Functions ]
[ Top ] [ m_hide_lapack ] [ Functions ]
NAME
xhegv_cplex
FUNCTION
xhegv_cplex computes all the eigenvalues, and optionally, the eigenvectors of a
(real generalized symmetric-definite| complex generalized Hermitian-definite)
eigenproblem, of the form
A*x=(lambda)*B*x (1),
A*Bx=(lambda)*x, (2), or
B*A*x=(lambda)*x (3).
Here A and B are assumed to be (symmetric|Hermitian) and B is also positive definite.
INPUTS
ITYPE (input) INTEGER Specifies the problem type to be solved:
= 1: A*x = (lambda)*B*x
= 2: A*B*x = (lambda)*x
= 3: B*A*x = (lambda)*x
JOBZ (input) CHARACTER*1
= "N": Compute eigenvalues only;
= "V": Compute eigenvalues and eigenvectors.
UPLO (input) CHARACTER*1
= "U": Upper triangle of A and B are stored;
= "L": Lower triangle of A and B are stored.
CPLEX Size of the first dimension of the A and B matrices.
1 for a real symmetric matrix.
2 for a complex Hermitian matrix.
N (input) INTEGER
The order of the matrices A and B. N >= 0.
[comm]=MPI communicator for ScaLAPACK inversion. Only available if the code has been compiled with Scalapack support.
To avoid wasting CPU time the scalapack initialization is avoided if the number of processors in 1,
in this case the sequential LAPACK routine is called.
OUTPUT
W (output) REAL(DP) array, dimension (N)
If INFO = 0, the eigenvalues in ascending order.
See also SIDE EFFECTS
SIDE EFFECTS
A (input/output) REAL(DP) array, dimension (CPLEX,N, N)
On entry, the (real symmetric|Hermitian) matrix A. If UPLO = "U", the leading N-by-N upper triangular part of A
<S-F1>contains the upper triangular part of the matrix A.
If UPLO = "L", the leading N-by-N lower triangular part of A contains the lower triangular part of the matrix A.
On exit, if JOBZ = "V", then A contains the matrix Z of eigenvectors.
The eigenvectors are normalized as follows:
if ITYPE =1 or 2:
Z**T*B*Z = I if CPLEX=1
Z**H*B*Z = I if CPLEX=2.
if ITYPE = 3,
Z**T*inv(B)*Z = I if CPLEX=1
Z**H*inv(B)*Z = I if CPLEX=2
If JOBZ = "N", then on exit the upper triangle (if UPLO="U") or the lower triangle
(if UPLO="L") of A, including the diagonal, is destroyed.
B (input/output) REAL(DP) array, dimension (CPLEX,N, N)
On entry, the (real symmetric|Hermitian) positive definite matrix B.
If UPLO = "U", the leading N-by-N upper triangular part of B contains the upper triangular part of the matrix B.
If UPLO = "L", the leading N-by-N lower triangular part of B contains the lower triangular part of the matrix B.
On exit, if INFO <= N, the part of B containing the matrix is overwritten by the triangular
factor U or L from the Cholesky factorization
B = U**T*U or B = L*L**T if CPLEX=1
B = U**H*U or B = L*L**H if CPLEX=2
SOURCE
1099 subroutine xhegv_cplex(itype, jobz, uplo, cplex, n, a, b, w, msg, ierr, comm) 1100 1101 !Arguments ------------------------------------ 1102 !scalars 1103 integer,intent(in) :: n,itype,cplex 1104 character(len=*),intent(in) :: jobz, uplo 1105 character(len=*),intent(out) :: msg 1106 integer,intent(out) :: ierr 1107 integer,optional,intent(in) :: comm 1108 !arrays 1109 real(dp),intent(inout) :: a(cplex,n,n), b(cplex,n,n) 1110 real(dp),intent(out) :: w(n) 1111 1112 !Local variables ------------------------------ 1113 !scalars 1114 integer :: lwork,nprocs,ii 1115 logical :: use_scalapack 1116 !arrays 1117 real(dp),allocatable :: rwork(:), work_real(:) 1118 complex(dpc),allocatable :: work_cplx(:) 1119 #ifdef HAVE_LINALG_SCALAPACK 1120 !integer :: istwf_k 1121 !type(matrix_scalapack) :: Slk_matA,Slk_matB 1122 !type(processor_scalapack) :: Slk_processor 1123 #endif 1124 !************************************************************************ 1125 1126 use_scalapack = .FALSE. 1127 if (present(comm)) then 1128 nprocs = xmpi_comm_size(comm) 1129 #ifdef HAVE_LINALG_SCALAPACK 1130 use_scalapack = nprocs > 1 1131 #endif 1132 end if 1133 1134 if (all(cplex /= [1, 2])) then 1135 write(msg,'(a,i0)')"Wrong value for cplex: ",cplex 1136 ierr = 1; return 1137 end if 1138 1139 SELECT CASE(use_scalapack) 1140 CASE (.FALSE.) 1141 1142 if (cplex==1) then 1143 ! Real symmetric case. 1144 lwork = MAX(1,3*n-1) 1145 1146 ABI_MALLOC(work_real,(lwork)) 1147 call DSYGV(itype,jobz,uplo,n,a,n,b,n,w,work_real,lwork,ierr) 1148 1149 if (ierr < 0) then 1150 write(msg,'(a,i0,a)')" The ",-ierr,"-th argument of DSYGV had an illegal value." 1151 end if 1152 1153 if (ierr > 0) then 1154 if (ierr <= n) then 1155 write(msg,'(2a,i0,a)')& 1156 " DSYGV failed to converge: ",ch10,& 1157 ierr," off-diagonal elements of an intermediate tridiagonal form did not converge to zero. " 1158 else 1159 ii = ierr - n 1160 write(msg,'(3a,i0,3a)')& 1161 "DSYGV failed to converge: ",ch10,& 1162 "The leading minor of order ",ii," of B is not positive definite. ",ch10,& 1163 "The factorization of B could not be completed and no eigenvalues or eigenvectors were computed." 1164 end if 1165 end if 1166 1167 ABI_FREE(work_real) 1168 return 1169 1170 else 1171 ! complex Hermitian case 1172 lwork = MAX(1,2*n-1) 1173 1174 ABI_MALLOC(work_cplx,(lwork)) 1175 ABI_MALLOC(rwork,(MAX(1,3*n-2))) 1176 1177 call ZHEGV(itype,jobz,uplo,n,a,n,b,n,w,work_cplx,lwork,rwork,ierr) 1178 1179 if (ierr < 0) then 1180 write(msg,'(a,i0,a)')" The ",-ierr,"-th argument of ZHEGV had an illegal value." 1181 end if 1182 1183 if (ierr > 0) then 1184 if (ierr <= n) then 1185 write(msg,'(2a,i0,a)')& 1186 "ZHEGV failed to converge: ",ch10,& 1187 ierr," off-diagonal elements of an intermediate tridiagonal form did not converge to zero. " 1188 else 1189 ii = ierr -n 1190 write(msg,'(3a,i0,3a)')& 1191 "ZHEGV failed to converge: ",ch10,& 1192 "The leading minor of order ",ii," of B is not positive definite. ",ch10,& 1193 "The factorization of B could not be completed and no eigenvalues or eigenvectors were computed." 1194 end if 1195 end if 1196 1197 ABI_FREE(rwork) 1198 ABI_FREE(work_cplx) 1199 return 1200 end if ! cplex 1201 1202 CASE (.TRUE.) 1203 1204 #ifdef HAVE_LINALG_SCALAPACK 1205 1206 ABI_ERROR("Not coded yet") 1207 ! call Slk_processor%init(comm) 1208 ! istwf_k=1 1209 1210 ! ! Initialize and fill Scalapack matrix from the global one. 1211 ! call Slk_matA%init(n,n,Slk_processor,istwf_k) 1212 ! call slk_matrix_from_global_dpc_2D(Slk_matA,uplo,a) 1213 1214 ! call Slk_matB%init(n,n,Slk_processor,istwf_k) 1215 ! call slk_matrix_from_global_dpc_2D(Slk_matB,uplo,b) 1216 1217 ! ! Solve the problem with scaLAPACK. 1218 ! ABI_ERROR("slk_pZHEGV not yet coded") 1219 ! ! TODO 1220 ! call slk_pzhegv(itype,jobz,uplo,Slk_matA,Slk_matB,w) 1221 1222 ! call Slk_matB%free() 1223 ! 1224 ! if (firstchar(jobz,(/"V","v"/))) then ! A is overwritten with the eigenvectors 1225 ! a = czero 1226 ! call slk_matrix_to_global_dpc_2D(Slk_matA,"All",a) ! Fill the entries calculated by this node. 1227 ! call xmpi_sum(a,comm,ierr) ! Fill the remaing entries of the global matrix 1228 ! end if 1229 1230 ! call Slk_matA%free() 1231 1232 ! call Slk_processor%free() 1233 1234 RETURN 1235 #endif 1236 1237 ABI_BUG("You should not be here!") 1238 END SELECT 1239 1240 end subroutine xhegv_cplex
m_hide_lapack/xhegvx_cplex [ Functions ]
[ Top ] [ m_hide_lapack ] [ Functions ]
NAME
xhegvx_cplex
FUNCTION
xhegvx_cplex - compute selected eigenvalues, and optionally, eigenvectors of a (real symmetric-definite|complex generalized Hermitian-definite) eigenproblem, of the form A*x=(lambda)*B*x, A*Bx=(lambda)*x, or B*A*x=(lambda)*x. Here A and B are assumed to be (real symmetric|complex Hermitian) and B is also positive definite. Eigenvalues and eigenvectors can be selected by specifying either a range of values or a range of indices for the desired eigenvalues.
INPUTS
ITYPE (input) INTEGER Specifies the problem type to be solved:
= 1: A*x = (lambda)*B*x
= 2: A*B*x = (lambda)*x
= 3: B*A*x = (lambda)*x
JOBZ (input) CHARACTER*1
= 'N': Compute eigenvalues only;
= 'V': Compute eigenvalues and eigenvectors.
RANGE (input) CHARACTER*1
= 'A': all eigenvalues will be found.
= 'V': all eigenvalues in the half-open interval (VL,VU]
will be found.
= 'I': the IL-th through IU-th eigenvalues will be found.
UPLO (input) CHARACTER*1
= 'U': Upper triangle of A is stored;
= 'L': Lower triangle of A is stored.
CPLEX Size of the first dimension of the matrices A and B
1 for Real symmetric matrices
2 for complex Hermitianmatrices
N (input) INTEGER
The order of the matrices A and B. N >= 0.
LDA (input) INTEGER
The leading dimension of the array A. LDA >= max(1,N).
VL (input) REAL(DP)
VU (input) REAL(DP)
If RANGE='V', the lower and upper bounds of the interval to
be searched for eigenvalues. VL < VU.
Not referenced if RANGE = 'A' or 'I'.
IL (input) INTEGER
IU (input) INTEGER
If RANGE='I', the indices (in ascending order) of the
smallest and largest eigenvalues to be returned.
1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
Not referenced if RANGE = 'A' or 'V'.
ABSTOL (input) REAL(DP)
The absolute error tolerance for the eigenvalues.
An approximate eigenvalue is accepted as converged
when it is determined to lie in an interval [a,b]
of width less than or equal to
ABSTOL + EPS * max( |a|,|b| ) ,
where EPS is the machine precision. If ABSTOL is less than
or equal to zero, then EPS*|T| will be used in its place,
where |T| is the 1-norm of the tridiagonal matrix obtained
by reducing A to tridiagonal form.
Eigenvalues will be computed most accurately when ABSTOL is
set to twice the underflow threshold 2*DLAMCH('S'), not zero.
If this routine returns with INFO>0, indicating that some
eigenvectors did not converge, try setting ABSTOL to
2*DLAMCH('S').
LDZ (input) INTEGER
The leading dimension of the array Z. LDZ >= 1, and if
JOBZ = 'V', LDZ >= max(1,N).
[comm]=MPI communicator for ScaLAPACK inversion. Only available if the code has been compiled with Scalapack support.
To avoid wasting CPU time the scalapack initialization is avoided if the number of processors in 1,
in this case the sequential LAPACK routine is called.
OUTPUT
M (output) INTEGER
The total number of eigenvalues found. 0 <= M <= N.
If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
W (output) REAL(DP) array, dimension (N)
On normal exit, the first M elements contain the selected
eigenvalues in ascending order.
Z (output) REAL(DP) array, dimension (CPLEX ,LDZ, max(1,M))
If JOBZ = 'V', then if INFO = 0, the first M columns of Z
contain the orthonormal eigenvectors of the matrix A
corresponding to the selected eigenvalues, with the i-th
column of Z holding the eigenvector associated with W(i).
The eigenvectors are normalized as follows:
if ITYPE = 1 or 2, Z**T*B*Z = I; if ITYPE = 3, Z**T*inv(B)*Z = I.
If an eigenvector fails to converge, then that column of Z
contains the latest approximation to the eigenvector, and the
index of the eigenvector is returned in IFAIL.
If JOBZ = 'N', then Z is not referenced.
Note: the user must ensure that at least max(1,M) columns are
supplied in the array Z; if RANGE = 'V', the exact value of M
is not known in advance and an upper bound must be used.
See also SIDE EFFECTS
SIDE EFFECTS
A (input/output) REAL(DP) array, dimension (CPLEX, N, N)
On entry, the (real symmetric| complex Hermitian) matrix A. If UPLO = 'U', the
leading N-by-N upper triangular part of A contains the
upper triangular part of the matrix A. If UPLO = "L",
the leading N-by-N lower triangular part of A contains
the lower triangular part of the matrix A.
On exit, the lower triangle (if UPLO="L") or the upper
triangle (if UPLO="U") of A, including the diagonal, is
destroyed.
B (input/output) REAL(DP) array, dimension (CPLEX, LDB, N)
On entry, the (real symmetric| complex Hermitian) matrix B. If UPLO = "U", the leading N-by-N upper triangular part
of B contains the upper triangular part of the matrix B.
If UPLO = "L", the leading N-by-N lower triangular part of B contains the lower triangular part of the matrix B.
On exit, if INFO <= N, the part of B containing the matrix is overwritten by the triangular factor
U or L from the Cholesky factorization B = U**H*U or B = L*L**H.
SOURCE
2100 subroutine xhegvx_cplex(itype, jobz, range, uplo, cplex, n, a, b, & 2101 vl, vu, il, iu, abstol, m, w, z, ldz, msg, ierr, comm) 2102 2103 !Arguments ------------------------------------ 2104 !scalars 2105 integer,intent(in) :: il,iu,ldz,n,itype,cplex 2106 integer,optional,intent(in) :: comm 2107 integer,intent(inout) :: m 2108 integer,intent(out) :: ierr 2109 real(dp),intent(in) :: abstol,vl,vu 2110 character(len=*),intent(in) :: jobz,range,uplo 2111 character(len=*),intent(out) :: msg 2112 !arrays 2113 real(dp),intent(out) :: w(n) 2114 !real(dp),intent(out) :: z(cplex,ldz,n) 2115 real(dp),intent(out) :: z(cplex,ldz,m) 2116 real(dp),intent(inout) :: a(cplex,n,n),b(cplex,n,n) 2117 2118 !Local variables ------------------------------ 2119 !scalars 2120 integer :: lwork,nprocs,ii 2121 logical :: use_scalapack 2122 !arrays 2123 integer,allocatable :: ifail(:),iwork(:) 2124 real(dp),allocatable :: rwork(:) 2125 real(dp),allocatable :: work_real(:) 2126 complex(dpc),allocatable :: work_cplx(:) 2127 #ifdef HAVE_LINALG_SCALAPACK 2128 !integer :: istwf_k 2129 !logical :: want_eigenvectors 2130 !type(matrix_scalapack) :: Slk_matA,Slk_matB,Slk_vec 2131 !type(processor_scalapack) :: Slk_processor 2132 #endif 2133 2134 !************************************************************************ 2135 2136 use_scalapack=.FALSE. 2137 if (PRESENT(comm)) then 2138 nprocs = xmpi_comm_size(comm) 2139 #ifdef HAVE_LINALG_SCALAPACK 2140 use_scalapack = (nprocs>1) 2141 #endif 2142 end if 2143 2144 if (ALL(cplex/=(/1,2/))) then 2145 write(msg,'(a,i0)')" Wrong value for cplex: ",cplex 2146 ierr = 1; return 2147 end if 2148 2149 SELECT CASE(use_scalapack) 2150 2151 CASE (.FALSE.) 2152 ! Standard LAPACK call. 2153 if (cplex==1) then 2154 ! Real symmetric case 2155 lwork = MAX(1,8*n) 2156 2157 ABI_MALLOC(work_real,(lwork)) 2158 ABI_MALLOC(iwork,(5*n)) 2159 ABI_MALLOC(ifail,(n)) 2160 2161 call DSYGVX(itype,jobz,range,uplo,n,a,n,b,n,vl,vu,il,iu,abstol,m,w,z,ldz,work_real,lwork,iwork,ifail,ierr) 2162 2163 if (ierr < 0) then 2164 write(msg,'(a,i0,a)')" The ",-ierr,"-th argument of DSYGVX had an illegal value." 2165 end if 2166 2167 if (ierr > 0) then 2168 if (ierr<= n) then 2169 write(msg,'(a,i0,a)')" DSYGVX failed to converge: ",ierr," eigenvectors failed to converge. " 2170 else 2171 ii = ierr - n 2172 write(msg,'(3a,i0,3a)')& 2173 " DSYGVX failed to converge: ",ch10,& 2174 " The leading minor of order ",ii," of B is not positive definite. ",ch10,& 2175 " The factorization of B could not be completed and no eigenvalues or eigenvectors were computed." 2176 end if 2177 end if 2178 2179 ABI_FREE(iwork) 2180 ABI_FREE(ifail) 2181 ABI_FREE(work_real) 2182 RETURN 2183 2184 else 2185 ! Complex Hermitian case. 2186 lwork = MAX(1,2*n) 2187 2188 ABI_MALLOC(work_cplx,(lwork)) 2189 ABI_MALLOC(rwork,(7*n)) 2190 ABI_MALLOC(iwork,(5*n)) 2191 ABI_MALLOC(ifail,(n)) 2192 2193 !write(std_out,*)"Calling ZHEGVX" 2194 call ZHEGVX(itype,jobz,range,uplo,n,a,n,b,n,vl,vu,il,iu,abstol,m,w,z,ldz,work_cplx,lwork,rwork,iwork,ifail,ierr) 2195 2196 if (ierr < 0) then 2197 write(msg,'(a,i0,a)')"The ",-ierr,"-th argument of ZHEGVX had an illegal value." 2198 end if 2199 2200 if (ierr > 0) then 2201 if (ierr<= n) then 2202 write(msg,'(a,i0,a)')"ZHEGVX failed to converge: ",ierr," eigenvectors failed to converge. " 2203 else 2204 ii = ierr -n 2205 write(msg,'(3a,i0,3a)')& 2206 "ZHEEVX failed to converge: ",ch10,& 2207 "The leading minor of order ",ii," of B is not positive definite. ",ch10,& 2208 "The factorization of B could not be completed and no eigenvalues or eigenvectors were computed." 2209 end if 2210 end if 2211 2212 ABI_FREE(iwork) 2213 ABI_FREE(ifail) 2214 ABI_FREE(rwork) 2215 ABI_FREE(work_cplx) 2216 RETURN 2217 end if ! cplex 2218 2219 CASE (.TRUE.) 2220 2221 #ifdef HAVE_LINALG_SCALAPACK 2222 ABI_ERROR("not coded yet") 2223 ! call Slk_processor%init(comm) 2224 ! istwf_k=1 2225 2226 ! ! Initialize and fill Scalapack matrix from the global one. 2227 ! call Slk_matA%init(n,n,Slk_processor,istwf_k) 2228 ! call slk_matrix_from_global_dpc_2D(Slk_matA,uplo,a) 2229 2230 ! call Slk_matB%init(n,n,Slk_processor,istwf_k) 2231 ! call slk_matrix_from_global_dpc_2D(Slk_matB,uplo,b) 2232 2233 ! want_eigenvectors = firstchar(jobz,(/"V","v"/)) 2234 ! if (want_eigenvectors) then ! Initialize the distributed vectors. 2235 ! call Slk_vec%init(n,n,Slk_processor,istwf_k) 2236 ! end if 2237 2238 ! ! Solve the problem. 2239 ! ABI_ERROR("slk_pZHEGVX not coded yet") 2240 ! ! TODO write the scaLAPACK wrapper. 2241 ! call slk_pZHEGVX(itype,jobz,range,uplo,Slk_matA,Slk_matB,vl,vu,il,iu,abstol,Slk_vec,m,w) 2242 2243 ! call Slk_matA%free() 2244 ! call Slk_matB%free() 2245 ! 2246 ! if (want_eigenvectors) then ! A is overwritten with the eigenvectors 2247 ! z = czero 2248 ! call slk_matrix_to_global_dpc_2D(Slk_vec,"All",z) ! Fill the entries calculated by this node. 2249 ! call Slk_vec%free() 2250 ! call xmpi_sum(z,comm,ierr) ! Fill the remaing entries of the global matrix 2251 ! end if 2252 2253 ! call Slk_processor%free() 2254 2255 ! RETURN 2256 #endif 2257 2258 ABI_BUG("You should not be here!") 2259 END SELECT 2260 2261 end subroutine xhegvx_cplex
m_hide_lapack/zginv [ Functions ]
[ Top ] [ m_hide_lapack ] [ Functions ]
NAME
zginv
FUNCTION
Invert a general matrix of complex elements in double precision. ZGETRF computes an LU factorization of a general N-by-N matrix A using partial pivoting with row interchanges. ZGETRI computes the inverse of a matrix using the LU factorization computed by ZGETRF.
INPUTS
n=size of complex matrix a a=matrix of complex elements [comm]=MPI communicator for ScaLAPACK inversion. Only available if the code has been compiled with Scalapack support. To avoid wasting CPU time the scalapack initialization is avoided if the number of processors in 1. In this case the sequential LAPACK routine is called.
SIDE EFFECTS
a(n,n)= array of complex elements, input, inverted at output
SOURCE
2708 subroutine zginv(a, n, comm) 2709 2710 !Arguments ------------------------------------ 2711 !scalars 2712 integer,intent(in) :: n 2713 integer,optional,intent(in) :: comm 2714 !arrays 2715 complex(dpc),intent(inout) :: a(n,n) 2716 2717 !Local variables------------------------------- 2718 !scalars 2719 integer :: lwork,info,nprocs 2720 logical :: use_scalapack 2721 character(len=500) :: msg 2722 !arrays 2723 integer,allocatable :: ipiv(:) 2724 complex(dpc),allocatable :: work(:) 2725 #ifdef HAVE_LINALG_SCALAPACK 2726 integer :: istwf_k,ierr 2727 type(matrix_scalapack) :: Slk_mat 2728 type(processor_scalapack) :: Slk_processor 2729 #endif 2730 2731 ! ************************************************************************* 2732 2733 use_scalapack=.FALSE. 2734 if (PRESENT(comm)) then 2735 nprocs = xmpi_comm_size(comm) 2736 #ifdef HAVE_LINALG_SCALAPACK 2737 use_scalapack = (nprocs>1) 2738 #endif 2739 end if 2740 2741 SELECT CASE(use_scalapack) 2742 CASE (.FALSE.) 2743 ABI_MALLOC(ipiv, (n)) 2744 call ZGETRF(n,n,a,n,ipiv,info) ! P* L* U Factorization. 2745 2746 if (info < 0) then 2747 write(msg,'(a,i0,a)')" The ",-info,"-th argument of ZGETRF had an illegal value." 2748 ABI_ERROR(msg) 2749 end if 2750 2751 if (info > 0) then 2752 write(msg,'(3a,i0,4a)')& 2753 "The matrix that has been passed in argument is probably either singular or nearly singular.",ch10,& 2754 "U(i,i) in the P*L*U factorization is exactly zero for i = ",info,ch10,& 2755 "The factorization has been completed but the factor U is exactly singular.",ch10,& 2756 "Division by zero will occur if it is used to solve a system of equations." 2757 ABI_ERROR(msg) 2758 end if 2759 2760 lwork=MAX(1,n) 2761 ABI_MALLOC(work,(lwork)) 2762 2763 call ZGETRI(n,a,n,ipiv,work,lwork,info) ! Invert U and then compute inv(A) 2764 2765 if (info < 0) then 2766 write(msg,'(a,i0,a)')" The ",-info,"-th argument of ZGETRI had an illegal value." 2767 ABI_ERROR(msg) 2768 end if 2769 2770 if (info > 0) then 2771 write(msg,'(3a,i0,a)')& 2772 "The matrix that has been passed to this subroutine is probably either singular or nearly singular.",ch10,& 2773 "U(i,i) for i= ",info," is exactly zero; the matrix is singular and its inverse could not be computed." 2774 ABI_ERROR(msg) 2775 end if 2776 2777 ABI_FREE(ipiv) 2778 ABI_FREE(work) 2779 RETURN 2780 2781 CASE (.TRUE.) 2782 2783 #ifdef HAVE_LINALG_SCALAPACK 2784 call Slk_processor%init(comm) 2785 istwf_k=1 2786 2787 ! Initialize and fill Scalapack matrix from the global one. 2788 call Slk_mat%init(n,n,Slk_processor,istwf_k) 2789 call slk_matrix_from_global_dpc_2D(Slk_mat,"All",a) 2790 2791 ! Perform the calculation with scaLAPACK. 2792 call Slk_mat%invert() 2793 2794 ! Reconstruct the global matrix from the distributed one. 2795 a = czero 2796 call slk_matrix_to_global_dpc_2D(Slk_mat,"All",a) ! Fill the entries calculated by this node. 2797 call Slk_mat%free() 2798 2799 call xmpi_sum(a,comm,ierr) ! Fill the remaing entries of the global matrix 2800 call Slk_processor%free() 2801 2802 return 2803 #endif 2804 2805 ABI_BUG("You should not be here!") 2806 END SELECT 2807 2808 end subroutine zginv
m_hide_lapack/zhpd_invert [ Functions ]
[ Top ] [ m_hide_lapack ] [ Functions ]
NAME
zhpd_invert
FUNCTION
Invert a Hermitian positive definite matrix of complex elements in double precision.
INPUTS
uplo= 'U': Upper triangle of A is stored;
= 'L': Lower triangle of A is stored.
n=size of complex matrix a
a=matrix of complex elements
[comm]=MPI communicator for ScaLAPACK inversion. Only available if the code has been compiled with Scalapack support.
To avoid wasting CPU time the scalapack initialization is avoided if the number of processors in 1.
In this case the sequential LAPACK routine is called.
SIDE EFFECTS
a(n,n)=
On entry, the Hermitian matrix A. If UPLO = 'U', the leading
N-by-N upper triangular part of A contains the upper
triangular part of the matrix A, and the strictly lower
triangular part of A is not referenced. If UPLO = 'L', the
leading N-by-N lower triangular part of A contains the lower
triangular part of the matrix A, and the strictly upper
triangular part of A is not referenced.
On exit, the upper or lower triangle of the (Hermitian) inverse of A
SOURCE
2842 subroutine zhpd_invert(uplo, a, n, comm) 2843 2844 !Arguments ------------------------------------ 2845 !scalars 2846 character(len=*),intent(in) :: uplo 2847 integer,intent(in) :: n 2848 integer,optional,intent(in) :: comm 2849 !arrays 2850 complex(dpc),intent(inout) :: a(n,n) 2851 2852 !Local variables------------------------------- 2853 !scalars 2854 integer :: info,nprocs 2855 logical :: use_scalapack 2856 character(len=500) :: msg 2857 !arrays 2858 #ifdef HAVE_LINALG_SCALAPACK 2859 integer :: istwf_k,ierr 2860 type(matrix_scalapack) :: Slk_mat 2861 type(processor_scalapack) :: Slk_processor 2862 #endif 2863 2864 ! ************************************************************************* 2865 2866 use_scalapack=.FALSE. 2867 if (PRESENT(comm)) then 2868 nprocs = xmpi_comm_size(comm) 2869 #ifdef HAVE_LINALG_SCALAPACK 2870 use_scalapack = (nprocs>1) 2871 #endif 2872 end if 2873 2874 SELECT CASE(use_scalapack) 2875 CASE (.FALSE.) 2876 ! * ZPOTRF computes the Cholesky factorization of a complex Hermitian positive definite. 2877 ! * A = U**H * U, if UPLO = 'U', or 2878 ! * A = L * L**H, if UPLO = 'L', 2879 call ZPOTRF(uplo,n,a,n,info) 2880 2881 if (info < 0) then 2882 write(msg,'(a,i0,a)')" The ",-info,"-th argument of ZPOTRF had an illegal value." 2883 ABI_ERROR(msg) 2884 end if 2885 2886 if (info > 0) then 2887 write(msg,'(a,i0,3a)')& 2888 "The leading minor of order ",info," is not positive definite, ",ch10,& 2889 "and the factorization could not be completed." 2890 ABI_ERROR(msg) 2891 end if 2892 ! 2893 ! * ZPOTRI computes the inverse of a complex Hermitian positive definite 2894 ! * matrix A using the Cholesky factorization A = U**H*U or A = L*L**H 2895 ! * computed by ZPOTRF. 2896 ! * On exit, the upper or lower triangle of the (Hermitian) 2897 ! * inverse of A, overwriting the input factor U or L. 2898 call ZPOTRI(uplo,n,a,n,info) 2899 2900 if (info < 0) then 2901 write(msg,'(a,i0,a)')" The ",-info,"-th argument of ZPOTRI had an illegal value." 2902 ABI_ERROR(msg) 2903 end if 2904 2905 if (info > 0) then 2906 write(msg,'(a,2(1x,i0),a)')& 2907 "The ( ",info,info,")element of the factor U or L is zero, and the inverse could not be computed." 2908 ABI_ERROR(msg) 2909 end if 2910 2911 RETURN 2912 2913 CASE (.TRUE.) 2914 2915 #ifdef HAVE_LINALG_SCALAPACK 2916 call Slk_processor%init(comm) 2917 istwf_k=1 2918 2919 ! Initialize and fill Scalapack matrix from the global one. 2920 call Slk_mat%init(n,n,Slk_processor,istwf_k) 2921 call slk_matrix_from_global_dpc_2D(Slk_mat,uplo,a) 2922 2923 ! Perform the calculation with scaLAPACK. 2924 call Slk_mat%hpd_invert(uplo, full=.False.) 2925 2926 ! Reconstruct the global matrix from the distributed one. 2927 a = czero 2928 call slk_matrix_to_global_dpc_2D(Slk_mat,uplo,a) ! Fill the entries calculated by this node. 2929 call Slk_mat%free() 2930 2931 call xmpi_sum(a,comm,ierr) ! Fill the remaing entries of the global matrix 2932 call Slk_processor%free() 2933 2934 RETURN 2935 #endif 2936 2937 ABI_BUG("You should not be here!") 2938 END SELECT 2939 2940 end subroutine zhpd_invert