Branch data Line data Source code
1 : : !!****m* ABINIT/m_hide_lapack
2 : : !! NAME
3 : : !! m_hide_lapack
4 : : !!
5 : : !! FUNCTION
6 : : !! ABINIT Linear Algebra Subroutine Interfaces.
7 : : !!
8 : : !! This modules provides interfaces performing the overloading of commonly used Lapack routines.
9 : : !! The main purpose of this module is to create a layer between abinit routines and Lapack procedures.
10 : : !! This layer can be used to hide the parallel Scalapack version. In this case, only the MPI commutator
11 : : !! has to be provided in input as the wrapper will take care of the initialization of the Scalapack grid as
12 : : !! well as of the distribution of the matrix. Note that this allows one to reduce
13 : : !! the CPU time per processor but not the memory allocated since the entire matrix has to be provided in input.
14 : : !! The interfaces are very similar to the Lapack F77 version (neither F90 constructs nor
15 : : !! F90 assumed size arrays are used). The main simplification with respect to the F77 version
16 : : !! of Lapack is that the work arrays are allocated inside the wrapper with optimal size
17 : : !! thus reducing the number of input argcomm_scalapackuments that has to be passed.
18 : : !! Leading dimensions have been removed from the interface whenever possible.
19 : : !! In F90 one can pass the array descriptor if the routines should operate on a slice
20 : : !! of the local array (seldom done in abinit).
21 : : !! Using array descriptor is OK but will it likely slow-down the calculation as
22 : : !! some compilers perform a copy of the input-output data.
23 : : !! If efficiency is a concern, then the F77 call should be used
24 : : !!
25 : : !! COPYRIGHT
26 : : !! Copyright (C) 1992-2026 ABINIT group (MG, GMR, XG)
27 : : !! This file is distributed under the terms of the
28 : : !! GNU General Public License, see ~abinit/COPYING
29 : : !! or http://www.gnu.org/copyleft/gpl.txt .
30 : : !! For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
31 : : !!
32 : : !! TODO
33 : : !! 1) Use a function to define the size of the Scalapack block according to some heuristic method.
34 : : !! 2) Define a threshold below which Scalapack is not used although the MPI communicator is passed.
35 : : !! 3) Split MPI communicator for Scalapack (.i.e. use a max size for the Scalapack comm; see abi_linalg_init).
36 : : !! 4) On certain networks, xmpi_sum might crash due to the size of the MPI packet.
37 : : !! This problem should be solved in hide_mpi (Module containing a private global variable
38 : : !! defining a threshold above which the input array is split into smaller chunks.
39 : : !!
40 : : !! SOURCE
41 : :
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 : : use m_linalg_interfaces
56 : :
57 : : use m_time, only : cwtime
58 : : use m_fstrings, only : firstchar
59 : :
60 : : implicit none
61 : :
62 : : private
63 : :
64 : : ! ==========================
65 : : ! Complex Hermitian matrices
66 : : ! ==========================
67 : : ! The _cplex version receives matrices declared as arr(cplex,N,M)
68 : : ! and calls the complex/real version depending on cplex (used e.g. for istwf_k = 2)
69 : :
70 : : public :: xheev ! Computes all the eigenvalues and, optionally, eigenvectors of a complex Hermitian matrix.
71 : : public :: xheev_cplex
72 : :
73 : : public :: xhpev ! Computes all the eigenvalues and, optionally, eigenvectors of a complex Hermitian matrix
74 : : ! in packed storage (Scalapack version not available)
75 : :
76 : : public :: xhegv ! Compute all the eigenvalues, and optionally, the eigenvectors of a complex generalized
77 : : ! Hermitian-definite eigenproblem, of the form:
78 : : ! A*x=(lambda)*B*x, A*Bx=(lambda)*x, or B*A*x=(lambda)*x
79 : :
80 : : public :: xheevx ! Computes selected eigenvalues and, optionally, eigenvectors of a complex Hermitian matrix A.
81 : : ! Eigenvalues and eigenvectors can be selected by specifying either a range of values or a range of
82 : : ! indices for the desired eigenvalues.
83 : :
84 : : public :: xheevx_cplex
85 : :
86 : : public :: xhegvx ! Computes selected eigenvalues, and optionally, eigenvectors
87 : : ! of a complex generalized Hermitian-definite eigenproblem,
88 : : ! of the form A*x=(lambda)*B*x, A*Bx=(lambda)*x, or B*A*x=(lambda)*x.
89 : : ! Eigenvalues and eigenvectors can be selected by specifying either a range of values or a range of
90 : : ! indices for the desired eigenvalues.
91 : : public :: xhegvx_cplex
92 : :
93 : :
94 : : 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.
95 : :
96 : :
97 : : ! ==============================
98 : : ! Complex non-symmetric matrices
99 : : ! ==============================
100 : :
101 : : public :: xgeev ! Computes for a complex nonsymmetric matrix A, the eigenvalues and, optionally,
102 : : ! the left and/or right eigenvectors.
103 : :
104 : : public :: xginv ! Invert a general matrix of complex elements by means of LU factorization.
105 : :
106 : :
107 : : public :: xhdp_invert ! Invert a Hermitian positive definite matrix.
108 : :
109 : : interface xheev
110 : : module procedure wrap_CHEEV
111 : : module procedure wrap_ZHEEV
112 : : end interface xheev
113 : :
114 : : interface xhpev
115 : : module procedure wrap_CHPEV
116 : : module procedure wrap_ZHPEV
117 : : end interface xhpev
118 : :
119 : : interface xhegv
120 : : module procedure wrap_ZHEGV
121 : : end interface xhegv
122 : :
123 : : public :: xhegv_cplex
124 : :
125 : : interface xheevx
126 : : module procedure wrap_ZHEEVX
127 : : end interface xheevx
128 : :
129 : : interface xhegvx
130 : : module procedure wrap_ZHEGVX
131 : : end interface xhegvx
132 : :
133 : : interface xgeev
134 : : module procedure wrap_CGEEV
135 : : module procedure wrap_ZGEEV
136 : : end interface xgeev
137 : :
138 : : interface xginv
139 : : module procedure cginv
140 : : module procedure zginv
141 : : end interface xginv
142 : :
143 : : interface xhdp_invert
144 : : module procedure zhpd_invert
145 : : end interface xhdp_invert
146 : :
147 : : public :: matrginv ! Invert a general matrix of real*8 elements.
148 : : public :: matr3eigval ! Find the eigenvalues of a real symmetric 3x3 matrix, entered in full storage mode.
149 : :
150 : : !FIXME This procedures are deprecated, use lapack API
151 : : public :: jacobi ! Computes all eigenvalues and eigenvectors of a real symmetric matrix a,
152 : : public :: ludcmp
153 : : public :: lubksb
154 : : public :: dzgedi
155 : : public :: dzgefa
156 : :
157 : : !----------------------------------------------------------------------
158 : : ! support for unitary tests and profiling.
159 : :
160 : : type,public :: latime_t
161 : : character(len=500) :: testname
162 : : integer :: msize
163 : : real(dp) :: ctime
164 : : real(dp) :: wtime
165 : : real(dp) :: max_abserr=-one
166 : : real(dp) :: gflops
167 : : end type latime_t
168 : :
169 : : public :: test_xginv
170 : :
171 : : !----------------------------------------------------------------------
172 : : ! private variables
173 : :
174 : : integer,private,parameter :: SLK_BLOCK_SIZE = 24
175 : : ! Default block size for Scalapack distribution.
176 : : ! As recommended by Intel MKL, a more sensible default than the previous value of 40
177 : :
178 : : CONTAINS !=========================================================================================================================
179 : : !!***
180 : :
181 : : !!****f* m_hide_lapack/wrap_CHEEV
182 : : !! NAME
183 : : !! wrap_CHEEV
184 : : !!
185 : : !! FUNCTION
186 : : !! wrap_CHEEV computes the eigenvalues and, optionally, the eigenvectors of a
187 : : !! complex Hermitian matrix in single precision. [PRIVATE]
188 : : !!
189 : : !! INPUTS
190 : : !! JOBZ (input) CHARACTER*1
191 : : !! = 'N': Compute eigenvalues only;
192 : : !! = 'V': Compute eigenvalues and eigenvectors.
193 : : !!
194 : : !! UPLO (input) CHARACTER*1
195 : : !! = 'U': Upper triangle of A is stored;
196 : : !! = 'L': Lower triangle of A is stored.
197 : : !!
198 : : !! N (input) INTEGER
199 : : !! The order of the matrix A. N >= 0.
200 : : !!
201 : : !! OUTPUT
202 : : !! W (output) REAL(SP) array, dimension (N)
203 : : !! If INFO = 0, the eigenvalues in ascending order.
204 : : !!
205 : : !! See also SIDE EFFECTS
206 : : !!
207 : : !! SIDE EFFECTS
208 : : !! A (input/output) COMPLEX(sp) array, dimension (N, N)
209 : : !! On entry, the Hermitian matrix A. If UPLO = 'U', the
210 : : !! leading N-by-N upper triangular part of A contains the
211 : : !! upper triangular part of the matrix A. If UPLO = 'L',
212 : : !! the leading N-by-N lower triangular part of A contains
213 : : !! the lower triangular part of the matrix A.
214 : : !! On exit, if JOBZ = 'V', then if INFO = 0, A contains the
215 : : !! orthonormal eigenvectors of the matrix A.
216 : : !! If JOBZ = 'N', then on exit the lower triangle (if UPLO='L')
217 : : !! or the upper triangle (if UPLO='U') of A, including the
218 : : !! diagonal, is destroyed.
219 : : !!
220 : : !! SOURCE
221 : :
222 : 0 : subroutine wrap_CHEEV(jobz, uplo, n, a, w)
223 : :
224 : : !Arguments ------------------------------------
225 : : !scalars
226 : : integer,intent(in) :: n
227 : : character(len=*),intent(in) :: jobz,uplo
228 : : !scalars
229 : : real(sp),intent(out) :: w(n)
230 : : complex(sp),intent(inout) :: a(n,n)
231 : :
232 : : !Local variables ------------------------------
233 : : !scalars
234 : : integer :: lwork,info
235 : : character(len=500) :: msg
236 : : !arrays
237 : 0 : real(sp),allocatable :: rwork(:)
238 : 0 : complex(sp),allocatable :: work(:)
239 : : !************************************************************************
240 : :
241 : 0 : lwork = MAX(1,2*n-1)
242 : :
243 [ # # # # : 0 : ABI_MALLOC(work, (lwork))
# # # # #
# ]
244 [ # # # # : 0 : ABI_MALLOC(rwork, (MAX(1,3*n-2)))
# # # # #
# ]
245 : :
246 : 0 : call CHEEV(jobz,uplo,n,a,n,w,work,lwork,rwork,info)
247 : :
248 [ # # ]: 0 : if (info < 0) then
249 : 0 : write(msg,'(a,i0,a)')"The ",-info,"-th argument of CHEEV had an illegal value."
250 : 0 : ABI_ERROR(msg)
251 : : end if
252 : :
253 [ # # ]: 0 : if (info > 0) then
254 : : write(msg,'(2a,i0,a)')&
255 : 0 : "CHEEV: the algorithm failed to converge; ",ch10,&
256 : 0 : info," off-diagonal elements of an intermediate tridiagonal form did not converge to zero. "
257 : 0 : ABI_ERROR(msg)
258 : : end if
259 : :
260 [ # # ]: 0 : ABI_FREE(rwork)
261 [ # # ]: 0 : ABI_FREE(work)
262 : :
263 : : !TODO scaLAPACK version (complex single precision buffer is needed in slkmat_dp_t)
264 : :
265 : 0 : end subroutine wrap_CHEEV
266 : : !!***
267 : :
268 : : !----------------------------------------------------------------------
269 : :
270 : : !!****f* m_hide_lapack/wrap_ZHEEV
271 : : !! NAME
272 : : !! wrap_ZHEEV
273 : : !!
274 : : !! FUNCTION
275 : : !! wrap_ZHEEV computes the eigenvalues and, optionally, the eigenvectors of a
276 : : !! complex Hermitian matrix in double precision. [PRIVATE]
277 : : !!
278 : : !! INPUTS
279 : : !! JOBZ (input) CHARACTER*1
280 : : !! = 'N': Compute eigenvalues only;
281 : : !! = 'V': Compute eigenvalues and eigenvectors.
282 : : !!
283 : : !! UPLO (input) CHARACTER*1
284 : : !! = 'U': Upper triangle of A is stored;
285 : : !! = 'L': Lower triangle of A is stored.
286 : : !!
287 : : !! N (input) INTEGER
288 : : !! The order of the matrix A. N >= 0.
289 : : !!
290 : : !! [comm]=MPI communicator for ScaLAPACK inversion. Only available if the code has been compiled with Scalapack support.
291 : : !! To avoid wasting CPU time the scalapack initialization is avoided if the number of processors in 1,
292 : : !! in this case the sequential LAPACK routine is called.
293 : : !! OUTPUT
294 : : !! W (output) REAL(DP) array, dimension (N)
295 : : !! If INFO = 0, the eigenvalues in ascending order.
296 : : !!
297 : : !! See also SIDE EFFECTS
298 : : !!
299 : : !! SIDE EFFECTS
300 : : !! A (input/output) COMPLEX(DP) array, dimension (N, N)
301 : : !! On entry, the Hermitian matrix A. If UPLO = 'U', the
302 : : !! leading N-by-N upper triangular part of A contains the
303 : : !! upper triangular part of the matrix A. If UPLO = 'L',
304 : : !! the leading N-by-N lower triangular part of A contains
305 : : !! the lower triangular part of the matrix A.
306 : : !! On exit, if JOBZ = 'V', then if INFO = 0, A contains the
307 : : !! orthonormal eigenvectors of the matrix A.
308 : : !! If JOBZ = 'N', then on exit the lower triangle (if UPLO='L')
309 : : !! or the upper triangle (if UPLO='U') of A, including the
310 : : !! diagonal, is destroyed.
311 : : !!
312 : : !! SOURCE
313 : :
314 : 52400 : subroutine wrap_ZHEEV(jobz, uplo, n, a, w, comm)
315 : :
316 : : !Arguments ------------------------------------
317 : : !scalars
318 : : integer,intent(in) :: n
319 : : integer,optional,intent(in) :: comm
320 : : character(len=*),intent(in) :: jobz,uplo
321 : : !arrays
322 : : complex(dp),intent(inout) :: a(n,n)
323 : : real(dp),intent(out) :: w(n)
324 : :
325 : : !Local variables ------------------------------
326 : : !scalars
327 : : integer :: lwork,info,nprocs
328 : : logical :: use_scalapack
329 : : character(len=500) :: msg
330 : : !arrays
331 : 52400 : real(dp),allocatable :: rwork(:)
332 : 52400 : complex(dp),allocatable :: work(:)
333 : : #ifdef HAVE_LINALG_SCALAPACK
334 : : integer :: ierr,istwf_k
335 : : logical :: want_eigenvectors
336 : : type(slkmat_dp_t) :: Slk_mat,Slk_vec
337 : : type(slk_processor_t) :: Slk_processor
338 : : #endif
339 : : !************************************************************************
340 : :
341 : 52400 : use_scalapack=.FALSE.
342 [ - + ]: 52400 : if (PRESENT(comm)) then
343 : 0 : nprocs = xmpi_comm_size(comm)
344 : : #ifdef HAVE_LINALG_SCALAPACK
345 : : use_scalapack = (nprocs>1)
346 : : #endif
347 : : end if
348 : :
349 : : SELECT CASE(use_scalapack)
350 : : CASE (.FALSE.)
351 : :
352 : 52400 : lwork = MAX(1,2*n-1)
353 [ - + + - : 157200 : ABI_MALLOC(work, (lwork))
+ - - + -
+ ]
354 [ - + + - : 157200 : ABI_MALLOC(rwork, (MAX(1,3*n-2)))
+ - - + -
+ ]
355 : :
356 : 52400 : call ZHEEV(jobz,uplo,n,a,n,w,work,lwork,rwork,info)
357 : :
358 [ - + ]: 52400 : if (info < 0) then
359 : 0 : write(msg,'(a,i0,a)')" The ",-info,"-th argument of ZHEEV had an illegal value."
360 : 0 : ABI_ERROR(msg)
361 : : end if
362 : :
363 [ - + ]: 52400 : if (info > 0) then
364 : : write(msg,'(2a,i0,a)')&
365 : 0 : "ZHEEV: the algorithm failed to converge; ",ch10,&
366 : 0 : info," off-diagonal elements of an intermediate tridiagonal form did not converge to zero. "
367 : 0 : ABI_ERROR(msg)
368 : : end if
369 : :
370 [ - + ]: 52400 : ABI_FREE(rwork)
371 [ - + ]: 52400 : ABI_FREE(work)
372 : : RETURN
373 : :
374 : : CASE (.TRUE.)
375 : : #ifdef HAVE_LINALG_SCALAPACK
376 : : call Slk_processor%init(comm)
377 : : istwf_k=1
378 : :
379 : : ! Initialize and fill Scalapack matrix from the global one.
380 : : call Slk_mat%init(n,n,Slk_processor,istwf_k)
381 : : call slk_matrix_from_global_dpc_2D(Slk_mat,uplo,a)
382 : :
383 : : want_eigenvectors = firstchar(jobz,(/"V","v"/))
384 : : if (want_eigenvectors) then
385 : : ! Initialize the distributed vectors.
386 : : call Slk_vec%init(n,n,Slk_processor,istwf_k)
387 : : end if
388 : :
389 : : ! Solve the problem with scaLAPACK.
390 : : call slk_mat%heev(jobz, uplo, Slk_vec, w)
391 : : call Slk_mat%free()
392 : :
393 : : if (want_eigenvectors) then ! A is overwritten with the eigenvectors
394 : : a = czero
395 : : call slk_matrix_to_global_dpc_2D(Slk_vec,"All",a) ! Fill the entries calculated by this node.
396 : : call Slk_vec%free()
397 : : call xmpi_sum(a,comm,ierr) ! Fill the remaining entries of the global matrix
398 : : end if
399 : :
400 : : call Slk_processor%free()
401 : : RETURN
402 : : #endif
403 : :
404 : : ABI_BUG("You should not be here!")
405 : : END SELECT
406 : :
407 : : end subroutine wrap_ZHEEV
408 : : !!***
409 : :
410 : : !----------------------------------------------------------------------
411 : :
412 : : !!****f* m_hide_lapack/xheev_cplex
413 : : !! NAME
414 : : !! xheev_cplex
415 : : !!
416 : : !! FUNCTION
417 : : !! xheev_cplex computes the eigenvalues and, optionally, the eigenvectors of a
418 : : !! (complex Hermitian| real symmetric) matrix in double precision.
419 : : !!
420 : : !! INPUTS
421 : : !! JOBZ (input) CHARACTER*1
422 : : !! = 'N': Compute eigenvalues only;
423 : : !! = 'V': Compute eigenvalues and eigenvectors.
424 : : !!
425 : : !! UPLO (input) CHARACTER*1
426 : : !! = 'U': Upper triangle of A is stored;
427 : : !! = 'L': Lower triangle of A is stored.
428 : : !!
429 : : !! CPLEX= Size of the first dimension of the A matrix.
430 : : !! 1 for a real symmetric matrix.
431 : : !! 2 for complex Hermitian matrix stored in a real array with real and imaginary part.
432 : : !!
433 : : !! N (input) INTEGER
434 : : !! The order of the matrix A. N >= 0.
435 : : !!
436 : : !! [comm]=MPI communicator for ScaLAPACK inversion. Only available if the code has been compiled with Scalapack support.
437 : : !! To avoid wasting CPU time the scalapack initialization is avoided if the number of processors in 1,
438 : : !! in this case the sequential LAPACK routine is called.
439 : : !!
440 : : !! OUTPUT
441 : : !! W (output) REAL(DP) array, dimension (N)
442 : : !! If INFO = 0, the eigenvalues in ascending order.
443 : : !!
444 : : !! See also SIDE EFFECTS
445 : : !!
446 : : !! SIDE EFFECTS
447 : : !! A (input/output) REAL(DP) array, dimension (CPLEX, N, N)
448 : : !! On entry, the (complex Hermitian|Real symmetric) matrix A. If UPLO = 'U', the
449 : : !! leading N-by-N upper triangular part of A contains the
450 : : !! upper triangular part of the matrix A. If UPLO = 'L',
451 : : !! the leading N-by-N lower triangular part of A contains
452 : : !! the lower triangular part of the matrix A.
453 : : !! On exit, if JOBZ = 'V', then if INFO = 0, A contains the
454 : : !! orthonormal eigenvectors of the matrix A.
455 : : !! If JOBZ = 'N', then on exit the lower triangle (if UPLO='L')
456 : : !! or the upper triangle (if UPLO='U') of A, including the
457 : : !! diagonal, is destroyed.
458 : : !!
459 : : !! SOURCE
460 : :
461 : 0 : subroutine xheev_cplex(jobz, uplo, cplex, n, a, w, msg, ierr, comm)
462 : :
463 : : !Arguments ------------------------------------
464 : : !scalars
465 : : integer,intent(in) :: n,cplex
466 : : integer,optional,intent(in) :: comm
467 : : character(len=*),intent(in) :: jobz,uplo
468 : : integer,intent(out) :: ierr
469 : : character(len=*),intent(out) :: msg
470 : : !arrays
471 : : real(dp),intent(inout) :: a(cplex,n,n)
472 : : real(dp),intent(out) :: w(n)
473 : :
474 : : !Local variables ------------------------------
475 : : !scalars
476 : : integer :: lwork,nprocs
477 : : logical :: use_scalapack
478 : : !arrays
479 : 0 : real(dp),allocatable :: rwork(:)
480 : 0 : real(dp),allocatable :: work_real(:)
481 : 0 : complex(dp),allocatable :: work_cplx(:)
482 : : #ifdef HAVE_LINALG_SCALAPACK
483 : : !integer :: istwf_k
484 : : !logical :: want_eigenvectors
485 : : !type(slkmat_dp_t) :: Slk_mat,Slk_vec
486 : : !type(slk_processor_t) :: Slk_processor
487 : : #endif
488 : : !************************************************************************
489 : :
490 : 0 : use_scalapack=.FALSE.
491 [ # # ]: 0 : if (PRESENT(comm)) then
492 : 0 : nprocs = xmpi_comm_size(comm)
493 : : #ifdef HAVE_LINALG_SCALAPACK
494 : : use_scalapack = (nprocs>1)
495 : : #endif
496 : : end if
497 : :
498 [ # # ]: 0 : if (ALL(cplex/= [1, 2])) then
499 : 0 : write(msg,'(a,i0)')" Wrong value for cplex: ",cplex
500 : 0 : ierr = 1; return
501 : : end if
502 : :
503 : : SELECT CASE(use_scalapack)
504 : : CASE (.FALSE.)
505 [ # # ]: 0 : if (cplex==1) then
506 : : ! Real symmetric case.
507 : 0 : lwork = MAX(1,3*n-1)
508 [ # # # # : 0 : ABI_MALLOC(work_real,(lwork))
# # # # #
# ]
509 : :
510 : 0 : call DSYEV(jobz,uplo,n,a,n,w,work_real,lwork,ierr)
511 : :
512 [ # # ]: 0 : if (ierr < 0) then
513 : 0 : write(msg,'(a,i0,a)')" The ",-ierr,"-th argument of DSYEV had an illegal value."
514 : : end if
515 : :
516 [ # # ]: 0 : if (ierr > 0) then
517 : : write(msg,'(2a,i0,a)')&
518 : 0 : "DSYEV: the algorithm failed to converge; ",ch10,&
519 : 0 : ierr," off-diagonal elements of an intermediate tridiagonal form did not converge to zero. "
520 : : end if
521 : :
522 [ # # ]: 0 : ABI_FREE(work_real)
523 : 0 : RETURN
524 : :
525 : : else
526 : : ! Hermitian case.
527 : 0 : lwork = MAX(1,2*n-1)
528 : :
529 [ # # # # : 0 : ABI_MALLOC(work_cplx, (lwork))
# # # # #
# ]
530 [ # # # # : 0 : ABI_MALLOC(rwork, (MAX(1,3*n-2)))
# # # # #
# ]
531 : :
532 : 0 : call ZHEEV(jobz,uplo,n,a,n,w,work_cplx,lwork,rwork,ierr)
533 : :
534 [ # # ]: 0 : if (ierr < 0) then
535 : 0 : write(msg,'(a,i0,a)')" The ",-ierr,"-th argument of ZHEEV had an illegal value."
536 : : end if
537 : :
538 [ # # ]: 0 : if (ierr > 0) then
539 : : write(msg,'(2a,i0,a)')&
540 : 0 : "ZHEEV: the algorithm failed to converge; ",ch10,&
541 : 0 : ierr," off-diagonal elements of an intermediate tridiagonal form did not converge to zero. "
542 : : end if
543 : :
544 [ # # ]: 0 : ABI_FREE(rwork)
545 [ # # ]: 0 : ABI_FREE(work_cplx)
546 : 0 : RETURN
547 : : end if ! cplex
548 : :
549 : : CASE (.TRUE.)
550 : :
551 : : #ifdef HAVE_LINALG_SCALAPACK
552 : : ABI_ERROR("Not coded yet")
553 : :
554 : : !call Slk_processor%init(comm)
555 : : !istwf_k=1
556 : : !
557 : : !! Initialize and fill Scalapack matrix from the global one.
558 : : !call Slk_mat%init(n,n,Slk_processor,istwf_k)
559 : : !
560 : : !call slk_matrix_from_global_dpc_2D(Slk_mat,uplo,a)
561 : : !
562 : : !want_eigenvectors = firstchar(jobz,(/"V","v"/))
563 : : !if (want_eigenvectors) then ! Initialize the distributed vectors.
564 : : ! call Slk_vec%init(n,n,Slk_processor,istwf_k)
565 : : !end if
566 : : !
567 : : !! Solve the problem with scaLAPACK.
568 : : !call slk_mat%heev(jobz,uplo,Slk_vec,w)
569 : : !call Slk_mat%free()
570 : : !
571 : : !if (want_eigenvectors) then ! A is overwritten with the eigenvectors
572 : : ! a = czero
573 : : ! call slk_matrix_to_global_dpc_2D(Slk_vec,"All",a) ! Fill the entries calculated by this node.
574 : : ! call Slk_vec%free()
575 : : ! call xmpi_sum(a,comm,ierr) ! Fill the remaining entries of the global matrix
576 : : !end if
577 : : !
578 : : !call Slk_processor%free()
579 : :
580 : : RETURN
581 : : #endif
582 : :
583 : : ABI_BUG("You should not be here!")
584 : : END SELECT
585 : :
586 : : end subroutine xheev_cplex
587 : : !!***
588 : :
589 : : !----------------------------------------------------------------------
590 : :
591 : : !!****f* m_hide_lapack/wrap_CHPEV
592 : : !! NAME
593 : : !! wrap_CHPEV
594 : : !!
595 : : !! FUNCTION
596 : : !! wrap_CHPEV computes all the eigenvalues and, optionally, eigenvectors of a
597 : : !! complex Hermitian matrix in packed storage. Scalapack version is not available. [PRIVATE].
598 : : !!
599 : : !! INPUTS
600 : : !! JOBZ (input) CHARACTER*1
601 : : !! = 'N': Compute eigenvalues only;
602 : : !! = 'V': Compute eigenvalues and eigenvectors.
603 : : !!
604 : : !! UPLO (input) CHARACTER*1
605 : : !! = 'U': Upper triangle of A is stored;
606 : : !! = 'L': Lower triangle of A is stored.
607 : : !!
608 : : !! N (input) INTEGER
609 : : !! The order of the matrix A. N >= 0.
610 : : !!
611 : : !! LDZ (input) INTEGER
612 : : !! The leading dimension of the array Z. LDZ >= 1, and if
613 : : !! JOBZ = 'V', LDZ >= max(1,N).
614 : : !!
615 : : !! OUTPUT
616 : : !! W (output) REAL(SP) array, dimension (N)
617 : : !! If INFO = 0, the eigenvalues in ascending order.
618 : : !!
619 : : !! Z (output) COMPLEX(sp) array, dimension (LDZ, N)
620 : : !! If JOBZ = 'V', then if INFO = 0, Z contains the orthonormal
621 : : !! eigenvectors of the matrix A, with the i-th column of Z
622 : : !! holding the eigenvector associated with W(i).
623 : : !! If JOBZ = 'N', then Z is not referenced.
624 : : !!
625 : : !! See also SIDE EFFECTS
626 : : !!
627 : : !! SIDE EFFECTS
628 : : !!
629 : : !! AP (input/output) COMPLEX(sp) array, dimension (N*(N+1)/2)
630 : : !! On entry, the upper or lower triangle of the Hermitian matrix
631 : : !! A, packed columnwise in a linear array. The j-th column of A
632 : : !! is stored in the array AP as follows:
633 : : !! if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
634 : : !! if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n.
635 : : !!
636 : : !! On exit, AP is overwritten by values generated during the
637 : : !! reduction to tridiagonal form. If UPLO = 'U', the diagonal
638 : : !! and first superdiagonal of the tridiagonal matrix T overwrite
639 : : !! the corresponding elements of A, and if UPLO = 'L', the
640 : : !! diagonal and first subdiagonal of T overwrite the
641 : : !! corresponding elements of A.
642 : : !!
643 : : !! SOURCE
644 : :
645 : 0 : subroutine wrap_CHPEV(jobz, uplo, n, ap, w, z, ldz)
646 : :
647 : : !Arguments ------------------------------------
648 : : !scalars
649 : : integer,intent(in) :: n,ldz
650 : : character(len=*),intent(in) :: jobz,uplo
651 : : !arrays
652 : : real(sp),intent(out) :: w(n)
653 : : complex(sp),intent(inout) :: ap(n*(n+1)/2)
654 : : complex(sp),intent(out) :: z(ldz,n)
655 : :
656 : : !Local variables ------------------------------
657 : : !scalars
658 : : integer :: info
659 : : character(len=500) :: msg
660 : : !arrays
661 : 0 : real(sp),allocatable :: rwork(:)
662 : 0 : complex(sp),allocatable :: work(:)
663 : : !************************************************************************
664 : :
665 [ # # # # : 0 : ABI_MALLOC(work, (MAX(1,2*n-1)))
# # # # #
# ]
666 [ # # # # : 0 : ABI_MALLOC(rwork, (MAX(1,3*n-2)))
# # # # #
# ]
667 : :
668 : 0 : call CHPEV( JOBZ, UPLO, N, AP, W, Z, LDZ, WORK, RWORK, INFO )
669 : :
670 [ # # ]: 0 : if (info < 0) then
671 : 0 : write(msg,'(a,i0,a)')" The ",-info,"-th argument of ZHEEV had an illegal value."
672 : 0 : ABI_ERROR(msg)
673 : : end if
674 : :
675 [ # # ]: 0 : if (info > 0) then
676 : : write(msg,'(2a,i0,a)')&
677 : 0 : "ZHPEV: the algorithm failed to converge; ",ch10,&
678 : 0 : info," off-diagonal elements of an intermediate tridiagonal form did not converge to zero. "
679 : 0 : ABI_ERROR(msg)
680 : : end if
681 : :
682 [ # # ]: 0 : ABI_FREE(rwork)
683 [ # # ]: 0 : ABI_FREE(work)
684 : :
685 : 0 : end subroutine wrap_CHPEV
686 : : !!***
687 : :
688 : : !----------------------------------------------------------------------
689 : :
690 : : !!****f* m_hide_lapack/wrap_ZHPEV
691 : : !! NAME
692 : : !! wrap_ZHPEV
693 : : !!
694 : : !! FUNCTION
695 : : !! wrap_ZHPEV computes all the eigenvalues and, optionally, eigenvectors of a
696 : : !! complex Hermitian matrix in packed storage. Scalapack version is not available. [PRIVATE].
697 : : !!
698 : : !! INPUTS
699 : : !! JOBZ (input) CHARACTER*1
700 : : !! = 'N': Compute eigenvalues only;
701 : : !! = 'V': Compute eigenvalues and eigenvectors.
702 : : !!
703 : : !! UPLO (input) CHARACTER*1
704 : : !! = 'U': Upper triangle of A is stored;
705 : : !! = 'L': Lower triangle of A is stored.
706 : : !!
707 : : !! N (input) INTEGER
708 : : !! The order of the matrix A. N >= 0.
709 : : !!
710 : : !! LDZ (input) INTEGER
711 : : !! The leading dimension of the array Z. LDZ >= 1, and if
712 : : !! JOBZ = 'V', LDZ >= max(1,N).
713 : : !!
714 : : !! [comm]=MPI communicator for ScaLAPACK inversion. Only available if the code has been compiled with Scalapack support.
715 : : !! To avoid wasting CPU time the scalapack initialization is avoided if the number of processors in 1,
716 : : !! in this case the sequential LAPACK routine is called. Note that scalapack does not provide native
717 : : !! support for packed symmetric matrices. Therefore we have to distribute the full matrix among the nodes.
718 : : !! in order to perform the calculation in parallel.
719 : : !!
720 : : !! OUTPUT
721 : : !! W (output) REAL(DP) array, dimension (N)
722 : : !! If INFO = 0, the eigenvalues in ascending order.
723 : : !!
724 : : !! Z (output) COMPLEX(DP) array, dimension (LDZ, N)
725 : : !! If JOBZ = 'V', then if INFO = 0, Z contains the orthonormal
726 : : !! eigenvectors of the matrix A, with the i-th column of Z
727 : : !! holding the eigenvector associated with W(i).
728 : : !! If JOBZ = 'N', then Z is not referenced.
729 : : !!
730 : : !! See also SIDE EFFECTS
731 : : !!
732 : : !! SIDE EFFECTS
733 : : !!
734 : : !! AP (input/output) COMPLEX(DP) array, dimension (N*(N+1)/2)
735 : : !! On entry, the upper or lower triangle of the Hermitian matrix
736 : : !! A, packed columnwise in a linear array. The j-th column of A
737 : : !! is stored in the array AP as follows:
738 : : !! if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
739 : : !! if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n.
740 : : !!
741 : : !! On exit, AP is overwritten by values generated during the
742 : : !! reduction to tridiagonal form. If UPLO = 'U', the diagonal
743 : : !! and first superdiagonal of the tridiagonal matrix T overwrite
744 : : !! the corresponding elements of A, and if UPLO = 'L', the
745 : : !! diagonal and first subdiagonal of T overwrite the
746 : : !! corresponding elements of A. Unchanged if ScaLAPACK is used.
747 : : !!
748 : : !! SOURCE
749 : :
750 : 21 : subroutine wrap_ZHPEV(jobz, uplo, n, ap, w, z, ldz, comm)
751 : :
752 : : !Arguments ------------------------------------
753 : : !scalars
754 : : integer,intent(in) :: n,ldz
755 : : integer,optional,intent(in) :: comm
756 : : character(len=*),intent(in) :: jobz,uplo
757 : : !arrays
758 : : real(dp),intent(out) :: w(n)
759 : : complex(dp),intent(inout) :: ap(n*(n+1)/2)
760 : : complex(dp),intent(out) :: z(ldz,n)
761 : :
762 : : !Local variables ------------------------------
763 : : !scalars
764 : : integer :: info,nprocs
765 : : logical :: use_scalapack
766 : : character(len=500) :: msg
767 : : !arrays
768 : 21 : real(dp),allocatable :: rwork(:)
769 : 21 : complex(dp),allocatable :: work(:)
770 : : #ifdef HAVE_LINALG_SCALAPACK
771 : : integer :: ierr,istwf_k
772 : : logical :: want_eigenvectors
773 : : type(slkmat_dp_t) :: Slk_mat,Slk_vec
774 : : type(slk_processor_t) :: Slk_processor
775 : : #endif
776 : : !************************************************************************
777 : :
778 : 21 : use_scalapack=.FALSE.
779 [ - + ]: 21 : if (PRESENT(comm)) then
780 : 0 : nprocs = xmpi_comm_size(comm)
781 : : #ifdef HAVE_LINALG_SCALAPACK
782 : : use_scalapack = (nprocs>1)
783 : : #endif
784 : : end if
785 : :
786 : : SELECT CASE(use_scalapack)
787 : : CASE (.FALSE.)
788 [ - + + - : 63 : ABI_MALLOC(work, (MAX(1,2*n-1)))
+ - - + -
+ ]
789 [ - + + - : 63 : ABI_MALLOC(rwork, (MAX(1,3*n-2)))
+ - - + -
+ ]
790 : :
791 : 21 : call ZHPEV(jobz,uplo,n,ap,w,z,ldz,work,rwork,info)
792 : :
793 [ - + ]: 21 : if (info < 0) then
794 : 0 : write(msg,'(a,i0,a)')" The ",-info,"-th argument of ZHPEV had an illegal value."
795 : 0 : ABI_ERROR(msg)
796 : : end if
797 : :
798 [ - + ]: 21 : if (info > 0) then
799 : : write(msg,'(2a,i0,a)')&
800 : 0 : "ZHPEV: the algorithm failed to converge; ",ch10,&
801 : 0 : info," off-diagonal elements of an intermediate tridiagonal form did not converge to zero. "
802 : 0 : ABI_ERROR(msg)
803 : : end if
804 : :
805 [ - + ]: 21 : ABI_FREE(rwork)
806 [ - + ]: 21 : ABI_FREE(work)
807 : : RETURN
808 : :
809 : : CASE (.TRUE.)
810 : :
811 : : #ifdef HAVE_LINALG_SCALAPACK
812 : : call Slk_processor%init(comm)
813 : : istwf_k=1
814 : :
815 : : ! Initialize and fill Scalapack matrix from the global one.
816 : : call Slk_mat%init(n,n,Slk_processor,istwf_k)
817 : : call slk_matrix_from_global_dpc_1Dp(Slk_mat,uplo,ap)
818 : :
819 : : want_eigenvectors = firstchar(jobz,(/"V","v"/))
820 : : if (want_eigenvectors) then
821 : : ! Initialize the distributed vectors.
822 : : call Slk_vec%init(n,n,Slk_processor,istwf_k)
823 : : end if
824 : :
825 : : ! Solve the problem with scaLAPACK.
826 : : call slk_mat%heev(jobz,uplo,Slk_vec,w)
827 : : call Slk_mat%free()
828 : :
829 : : if (want_eigenvectors) then ! Collect the eigenvectors.
830 : : z = zero
831 : : call slk_matrix_to_global_dpc_2D(Slk_vec,"All",z) ! Fill the entries calculated by this node.
832 : : call Slk_vec%free()
833 : : call xmpi_sum(z,comm,ierr) ! Fill the remaining entries of the global matrix
834 : : end if
835 : :
836 : : call Slk_processor%free()
837 : :
838 : : RETURN
839 : : #endif
840 : :
841 : : ABI_BUG("You should not be here!")
842 : : END SELECT
843 : :
844 : : end subroutine wrap_ZHPEV
845 : : !!***
846 : :
847 : : !----------------------------------------------------------------------
848 : :
849 : : !!****f* m_hide_lapack/wrap_ZHEGV
850 : : !! NAME
851 : : !! wrap_ZHEGV
852 : : !!
853 : : !! FUNCTION
854 : : !! wrap_ZHEGV computes all the eigenvalues, and optionally, the eigenvectors of a complex generalized
855 : : !! Hermitian-definite eigenproblem, of the form
856 : : !! A*x=(lambda)*B*x (1),
857 : : !! A*Bx=(lambda)*x, (2), or
858 : : !! B*A*x=(lambda)*x (3).
859 : : !! Here A and B are assumed to be Hermitian and B is also positive definite.
860 : : !!
861 : : !! INPUTS
862 : : !! ITYPE (input) INTEGER Specifies the problem type to be solved:
863 : : !! = 1: A*x = (lambda)*B*x
864 : : !! = 2: A*B*x = (lambda)*x
865 : : !! = 3: B*A*x = (lambda)*x
866 : : !!
867 : : !! JOBZ (input) CHARACTER*1
868 : : !! = "N": Compute eigenvalues only;
869 : : !! = "V": Compute eigenvalues and eigenvectors.
870 : : !!
871 : : !! UPLO (input) CHARACTER*1
872 : : !! = "U": Upper triangle of A and B are stored;
873 : : !! = "L": Lower triangle of A and B are stored.
874 : : !!
875 : : !! N (input) INTEGER
876 : : !! The order of the matrices A and B. N >= 0.
877 : : !!
878 : : !! [comm]=MPI communicator for ScaLAPACK inversion. Only available if the code has been compiled with Scalapack support.
879 : : !! To avoid wasting CPU time the scalapack initialization is avoided if the number of processors in 1,
880 : : !! in this case the sequential LAPACK routine is called.
881 : : !!
882 : : !! OUTPUT
883 : : !! W (output) REAL(DP) array, dimension (N)
884 : : !! If INFO = 0, the eigenvalues in ascending order.
885 : : !!
886 : : !! See also SIDE EFFECTS
887 : : !!
888 : : !! SIDE EFFECTS
889 : : !! A (input/output) COMPLEX(DP) array, dimension (N, N)
890 : : !! On entry, the Hermitian matrix A. If UPLO = "U", the leading N-by-N upper triangular part of A
891 : : !! <S-F1>contains the upper triangular part of the matrix A.
892 : : !! If UPLO = "L", the leading N-by-N lower triangular part of A contains the lower triangular part of the matrix A.
893 : : !!
894 : : !! On exit, if JOBZ = "V", then A contains the matrix Z of eigenvectors.
895 : : !! The eigenvectors are normalized as follows: if ITYPE = 1 or 2,
896 : : !! Z**H*B*Z = I; if ITYPE = 3, Z**H*inv(B)*Z = I.
897 : : !! If JOBZ = "N", then on exit the upper triangle (if UPLO="U") or the lower triangle
898 : : !! (if UPLO="L") of A, including the diagonal, is destroyed.
899 : : !!
900 : : !! B (input/output) COMPLEX*16 array, dimension (LDB, N)
901 : : !! On entry, the Hermitian positive definite matrix B.
902 : : !! If UPLO = "U", the leading N-by-N upper triangular part of B contains the upper triangular part of the matrix B.
903 : : !! If UPLO = "L", the leading N-by-N lower triangular part of B contains the lower triangular part of the matrix B.
904 : : !!
905 : : !! On exit, if INFO <= N, the part of B containing the matrix is overwritten by the triangular
906 : : !! factor U or L from the Cholesky factorization B = U**H*U or B = L*L**H.
907 : : !!
908 : : !! SOURCE
909 : :
910 : 19 : subroutine wrap_ZHEGV(itype, jobz, uplo, n, a, b, w, comm)
911 : :
912 : : !Arguments ------------------------------------
913 : : !scalars
914 : : integer,intent(in) :: n,itype
915 : : integer,optional,intent(in) :: comm
916 : : character(len=*),intent(in) :: jobz,uplo
917 : : !arrays
918 : : complex(dp),intent(inout) :: a(n,n),b(n,n)
919 : : real(dp),intent(out) :: w(n)
920 : :
921 : : !Local variables ------------------------------
922 : : !scalars
923 : : integer :: lwork,info,nprocs,ii
924 : : logical :: use_scalapack
925 : : character(len=500) :: msg
926 : : !arrays
927 : 19 : real(dp),allocatable :: rwork(:)
928 : 19 : complex(dp),allocatable :: work(:)
929 : : #ifdef HAVE_LINALG_SCALAPACK
930 : : integer :: ierr,istwf_k
931 : : type(slkmat_dp_t) :: Slk_matA,Slk_matB
932 : : type(slk_processor_t) :: Slk_processor
933 : : #endif
934 : : !************************************************************************
935 : :
936 : 19 : use_scalapack=.FALSE.
937 [ - + ]: 19 : if (PRESENT(comm)) then
938 : 0 : nprocs = xmpi_comm_size(comm)
939 : : #ifdef HAVE_LINALG_SCALAPACK
940 : : use_scalapack = (nprocs>1)
941 : : #endif
942 : : end if
943 : :
944 : : SELECT CASE(use_scalapack)
945 : : CASE (.FALSE.)
946 : 19 : lwork = MAX(1,2*n-1)
947 : :
948 [ - + + - : 57 : ABI_MALLOC(work, (lwork))
+ - - + -
+ ]
949 [ - + + - : 57 : ABI_MALLOC(rwork, (MAX(1,3*n-2)))
+ - - + -
+ ]
950 : :
951 : 19 : call ZHEGV(itype,jobz,uplo,n,a,n,b,n,w,work,lwork,rwork,info)
952 : :
953 [ - + ]: 19 : if (info < 0) then
954 : 0 : write(msg,'(a,i0,a)')" The ",-info,"-th argument of ZHEGV had an illegal value."
955 : 0 : ABI_ERROR(msg)
956 : : end if
957 : :
958 [ - + ]: 19 : if (info > 0) then
959 [ # # ]: 0 : if (info<= n) then
960 : : write(msg,'(2a,i0,a)')&
961 : 0 : "ZHEGV failed to converge: ",ch10,&
962 : 0 : info," off-diagonal elements of an intermediate tridiagonal form did not converge to zero. "
963 : : else
964 : 0 : ii = info -n
965 : : write(msg,'(3a,i0,3a)')&
966 : 0 : "ZHEGV failed to converge: ",ch10,&
967 : 0 : "The leading minor of order ",ii," of B is not positive definite. ",ch10,&
968 : 0 : "The factorization of B could not be completed and no eigenvalues or eigenvectors were computed."
969 : : end if
970 : 0 : ABI_ERROR(msg)
971 : : end if
972 : :
973 [ - + ]: 19 : ABI_FREE(rwork)
974 [ - + ]: 19 : ABI_FREE(work)
975 : : RETURN
976 : :
977 : : CASE (.TRUE.)
978 : :
979 : : #ifdef HAVE_LINALG_SCALAPACK
980 : : call Slk_processor%init(comm)
981 : : istwf_k=1
982 : :
983 : : ! Initialize and fill Scalapack matrix from the global one.
984 : : call Slk_matA%init(n,n,Slk_processor,istwf_k)
985 : : call slk_matrix_from_global_dpc_2D(Slk_matA,uplo,a)
986 : :
987 : : call Slk_matB%init(n,n,Slk_processor,istwf_k)
988 : : call slk_matrix_from_global_dpc_2D(Slk_matB,uplo,b)
989 : :
990 : : ! Solve the problem with scaLAPACK.
991 : : ABI_ERROR("slk_pZHEGV not yet coded")
992 : : ! TODO
993 : : !% call slk_pzhegv(itype,jobz,uplo,Slk_matA,Slk_matB,w)
994 : :
995 : : call Slk_matB%free()
996 : :
997 : : if (firstchar(jobz,(/"V","v"/))) then ! A is overwritten with the eigenvectors
998 : : a = czero
999 : : call slk_matrix_to_global_dpc_2D(Slk_matA,"All",a) ! Fill the entries calculated by this node.
1000 : : call xmpi_sum(a,comm,ierr) ! Fill the remaining entries of the global matrix
1001 : : end if
1002 : :
1003 : : call Slk_matA%free()
1004 : : call Slk_processor%free()
1005 : : RETURN
1006 : : #endif
1007 : :
1008 : : ABI_BUG("You should not be here!")
1009 : : END SELECT
1010 : :
1011 : : end subroutine wrap_ZHEGV
1012 : : !!***
1013 : :
1014 : : !----------------------------------------------------------------------
1015 : :
1016 : : !!****f* m_hide_lapack/xhegv_cplex
1017 : : !! NAME
1018 : : !! xhegv_cplex
1019 : : !!
1020 : : !! FUNCTION
1021 : : !! xhegv_cplex computes all the eigenvalues, and optionally, the eigenvectors of a
1022 : : !! (real generalized symmetric-definite| complex generalized Hermitian-definite)
1023 : : !! eigenproblem, of the form
1024 : : !! A*x=(lambda)*B*x (1),
1025 : : !! A*Bx=(lambda)*x, (2), or
1026 : : !! B*A*x=(lambda)*x (3).
1027 : : !! Here A and B are assumed to be (symmetric|Hermitian) and B is also positive definite.
1028 : : !!
1029 : : !! INPUTS
1030 : : !! ITYPE (input) INTEGER Specifies the problem type to be solved:
1031 : : !! = 1: A*x = (lambda)*B*x
1032 : : !! = 2: A*B*x = (lambda)*x
1033 : : !! = 3: B*A*x = (lambda)*x
1034 : : !!
1035 : : !! JOBZ (input) CHARACTER*1
1036 : : !! = "N": Compute eigenvalues only;
1037 : : !! = "V": Compute eigenvalues and eigenvectors.
1038 : : !!
1039 : : !! UPLO (input) CHARACTER*1
1040 : : !! = "U": Upper triangle of A and B are stored;
1041 : : !! = "L": Lower triangle of A and B are stored.
1042 : : !!
1043 : : !! CPLEX Size of the first dimension of the A and B matrices.
1044 : : !! 1 for a real symmetric matrix.
1045 : : !! 2 for a complex Hermitian matrix.
1046 : : !!
1047 : : !! N (input) INTEGER
1048 : : !! The order of the matrices A and B. N >= 0.
1049 : : !!
1050 : : !! [comm]=MPI communicator for ScaLAPACK inversion. Only available if the code has been compiled with Scalapack support.
1051 : : !! To avoid wasting CPU time the scalapack initialization is avoided if the number of processors in 1,
1052 : : !! in this case the sequential LAPACK routine is called.
1053 : : !!
1054 : : !! OUTPUT
1055 : : !! W (output) REAL(DP) array, dimension (N)
1056 : : !! If INFO = 0, the eigenvalues in ascending order.
1057 : : !!
1058 : : !! See also SIDE EFFECTS
1059 : : !!
1060 : : !! SIDE EFFECTS
1061 : : !! A (input/output) REAL(DP) array, dimension (CPLEX,N, N)
1062 : : !! On entry, the (real symmetric|Hermitian) matrix A. If UPLO = "U", the leading N-by-N upper triangular part of A
1063 : : !! <S-F1>contains the upper triangular part of the matrix A.
1064 : : !! If UPLO = "L", the leading N-by-N lower triangular part of A contains the lower triangular part of the matrix A.
1065 : : !!
1066 : : !! On exit, if JOBZ = "V", then A contains the matrix Z of eigenvectors.
1067 : : !! The eigenvectors are normalized as follows:
1068 : : !! if ITYPE =1 or 2:
1069 : : !! Z**T*B*Z = I if CPLEX=1
1070 : : !! Z**H*B*Z = I if CPLEX=2.
1071 : : !! if ITYPE = 3,
1072 : : !! Z**T*inv(B)*Z = I if CPLEX=1
1073 : : !! Z**H*inv(B)*Z = I if CPLEX=2
1074 : : !!
1075 : : !! If JOBZ = "N", then on exit the upper triangle (if UPLO="U") or the lower triangle
1076 : : !! (if UPLO="L") of A, including the diagonal, is destroyed.
1077 : : !!
1078 : : !! B (input/output) REAL(DP) array, dimension (CPLEX,N, N)
1079 : : !! On entry, the (real symmetric|Hermitian) positive definite matrix B.
1080 : : !! If UPLO = "U", the leading N-by-N upper triangular part of B contains the upper triangular part of the matrix B.
1081 : : !! If UPLO = "L", the leading N-by-N lower triangular part of B contains the lower triangular part of the matrix B.
1082 : : !!
1083 : : !! On exit, if INFO <= N, the part of B containing the matrix is overwritten by the triangular
1084 : : !! factor U or L from the Cholesky factorization
1085 : : !! B = U**T*U or B = L*L**T if CPLEX=1
1086 : : !! B = U**H*U or B = L*L**H if CPLEX=2
1087 : : !!
1088 : : !! SOURCE
1089 : :
1090 : 0 : subroutine xhegv_cplex(itype, jobz, uplo, cplex, n, a, b, w, msg, ierr, comm)
1091 : :
1092 : : !Arguments ------------------------------------
1093 : : !scalars
1094 : : integer,intent(in) :: n,itype,cplex
1095 : : character(len=*),intent(in) :: jobz, uplo
1096 : : character(len=*),intent(out) :: msg
1097 : : integer,intent(out) :: ierr
1098 : : integer,optional,intent(in) :: comm
1099 : : !arrays
1100 : : real(dp),intent(inout) :: a(cplex,n,n), b(cplex,n,n)
1101 : : real(dp),intent(out) :: w(n)
1102 : :
1103 : : !Local variables ------------------------------
1104 : : !scalars
1105 : : integer :: lwork,nprocs,ii
1106 : : logical :: use_scalapack
1107 : : !arrays
1108 : 0 : real(dp),allocatable :: rwork(:), work_real(:)
1109 : 0 : complex(dp),allocatable :: work_cplx(:)
1110 : : #ifdef HAVE_LINALG_SCALAPACK
1111 : : !integer :: istwf_k
1112 : : !type(slkmat_dp_t) :: Slk_matA,Slk_matB
1113 : : !type(slk_processor_t) :: Slk_processor
1114 : : #endif
1115 : : !************************************************************************
1116 : :
1117 : 0 : use_scalapack = .FALSE.
1118 [ # # ]: 0 : if (present(comm)) then
1119 : 0 : nprocs = xmpi_comm_size(comm)
1120 : : #ifdef HAVE_LINALG_SCALAPACK
1121 : : use_scalapack = nprocs > 1
1122 : : #endif
1123 : : end if
1124 : :
1125 [ # # ]: 0 : if (all(cplex /= [1, 2])) then
1126 : 0 : write(msg,'(a,i0)')"Wrong value for cplex: ",cplex
1127 : 0 : ierr = 1; return
1128 : : end if
1129 : :
1130 : : SELECT CASE(use_scalapack)
1131 : : CASE (.FALSE.)
1132 : :
1133 [ # # ]: 0 : if (cplex==1) then
1134 : : ! Real symmetric case.
1135 : 0 : lwork = MAX(1,3*n-1)
1136 : :
1137 [ # # # # : 0 : ABI_MALLOC(work_real,(lwork))
# # # # #
# ]
1138 : 0 : call DSYGV(itype,jobz,uplo,n,a,n,b,n,w,work_real,lwork,ierr)
1139 : :
1140 [ # # ]: 0 : if (ierr < 0) then
1141 : 0 : write(msg,'(a,i0,a)')" The ",-ierr,"-th argument of DSYGV had an illegal value."
1142 : : end if
1143 : :
1144 [ # # ]: 0 : if (ierr > 0) then
1145 [ # # ]: 0 : if (ierr <= n) then
1146 : : write(msg,'(2a,i0,a)')&
1147 : 0 : " DSYGV failed to converge: ",ch10,&
1148 : 0 : ierr," off-diagonal elements of an intermediate tridiagonal form did not converge to zero. "
1149 : : else
1150 : 0 : ii = ierr - n
1151 : : write(msg,'(3a,i0,3a)')&
1152 : 0 : "DSYGV failed to converge: ",ch10,&
1153 : 0 : "The leading minor of order ",ii," of B is not positive definite. ",ch10,&
1154 : 0 : "The factorization of B could not be completed and no eigenvalues or eigenvectors were computed."
1155 : : end if
1156 : : end if
1157 : :
1158 [ # # ]: 0 : ABI_FREE(work_real)
1159 : 0 : return
1160 : :
1161 : : else
1162 : : ! complex Hermitian case
1163 : 0 : lwork = MAX(1,2*n-1)
1164 : :
1165 [ # # # # : 0 : ABI_MALLOC(work_cplx,(lwork))
# # # # #
# ]
1166 [ # # # # : 0 : ABI_MALLOC(rwork,(MAX(1,3*n-2)))
# # # # #
# ]
1167 : :
1168 : 0 : call ZHEGV(itype,jobz,uplo,n,a,n,b,n,w,work_cplx,lwork,rwork,ierr)
1169 : :
1170 [ # # ]: 0 : if (ierr < 0) then
1171 : 0 : write(msg,'(a,i0,a)')" The ",-ierr,"-th argument of ZHEGV had an illegal value."
1172 : : end if
1173 : :
1174 [ # # ]: 0 : if (ierr > 0) then
1175 [ # # ]: 0 : if (ierr <= n) then
1176 : : write(msg,'(2a,i0,a)')&
1177 : 0 : "ZHEGV failed to converge: ",ch10,&
1178 : 0 : ierr," off-diagonal elements of an intermediate tridiagonal form did not converge to zero. "
1179 : : else
1180 : 0 : ii = ierr -n
1181 : : write(msg,'(3a,i0,3a)')&
1182 : 0 : "ZHEGV failed to converge: ",ch10,&
1183 : 0 : "The leading minor of order ",ii," of B is not positive definite. ",ch10,&
1184 : 0 : "The factorization of B could not be completed and no eigenvalues or eigenvectors were computed."
1185 : : end if
1186 : : end if
1187 : :
1188 [ # # ]: 0 : ABI_FREE(rwork)
1189 [ # # ]: 0 : ABI_FREE(work_cplx)
1190 : 0 : return
1191 : : end if ! cplex
1192 : :
1193 : : CASE (.TRUE.)
1194 : :
1195 : : #ifdef HAVE_LINALG_SCALAPACK
1196 : :
1197 : : ABI_ERROR("Not coded yet")
1198 : : ! call Slk_processor%init(comm)
1199 : : ! istwf_k=1
1200 : :
1201 : : ! ! Initialize and fill Scalapack matrix from the global one.
1202 : : ! call Slk_matA%init(n,n,Slk_processor,istwf_k)
1203 : : ! call slk_matrix_from_global_dpc_2D(Slk_matA,uplo,a)
1204 : :
1205 : : ! call Slk_matB%init(n,n,Slk_processor,istwf_k)
1206 : : ! call slk_matrix_from_global_dpc_2D(Slk_matB,uplo,b)
1207 : :
1208 : : ! ! Solve the problem with scaLAPACK.
1209 : : ! ABI_ERROR("slk_pZHEGV not yet coded")
1210 : : ! ! TODO
1211 : : ! call slk_pzhegv(itype,jobz,uplo,Slk_matA,Slk_matB,w)
1212 : :
1213 : : ! call Slk_matB%free()
1214 : : !
1215 : : ! if (firstchar(jobz,(/"V","v"/))) then ! A is overwritten with the eigenvectors
1216 : : ! a = czero
1217 : : ! call slk_matrix_to_global_dpc_2D(Slk_matA,"All",a) ! Fill the entries calculated by this node.
1218 : : ! call xmpi_sum(a,comm,ierr) ! Fill the remaining entries of the global matrix
1219 : : ! end if
1220 : :
1221 : : ! call Slk_matA%free()
1222 : :
1223 : : ! call Slk_processor%free()
1224 : :
1225 : : RETURN
1226 : : #endif
1227 : :
1228 : : ABI_BUG("You should not be here!")
1229 : : END SELECT
1230 : :
1231 : : end subroutine xhegv_cplex
1232 : : !!***
1233 : :
1234 : : !----------------------------------------------------------------------
1235 : :
1236 : : !!****f* m_hide_lapack/wrap_ZHEEVX
1237 : : !! NAME
1238 : : !! wrap_ZHEEVX
1239 : : !!
1240 : : !! FUNCTION
1241 : : !! wrap_ZHEEVX computes selected eigenvalues and, optionally, eigenvectors
1242 : : !! of a complex Hermitian matrix A. Eigenvalues and eigenvectors can
1243 : : !! be selected by specifying either a range of values or a range of
1244 : : !! indices for the desired eigenvalues.
1245 : : !!
1246 : : !! INPUTS
1247 : : !! JOBZ (input) CHARACTER*1
1248 : : !! = 'N': Compute eigenvalues only;
1249 : : !! = 'V': Compute eigenvalues and eigenvectors.
1250 : : !!
1251 : : !! RANGE (input) CHARACTER*1
1252 : : !! = 'A': all eigenvalues will be found.
1253 : : !! = 'V': all eigenvalues in the half-open interval (VL,VU]
1254 : : !! will be found.
1255 : : !! = 'I': the IL-th through IU-th eigenvalues will be found.
1256 : : !!
1257 : : !! UPLO (input) CHARACTER*1
1258 : : !! = 'U': Upper triangle of A is stored;
1259 : : !! = 'L': Lower triangle of A is stored.
1260 : : !!
1261 : : !! N (input) INTEGER
1262 : : !! The order of the matrix A. N >= 0.
1263 : : !!
1264 : : !! LDA (input) INTEGER
1265 : : !! The leading dimension of the array A. LDA >= max(1,N).
1266 : : !!
1267 : : !! VL (input) REAL(DP)
1268 : : !! VU (input) REAL(DP)
1269 : : !! If RANGE='V', the lower and upper bounds of the interval to
1270 : : !! be searched for eigenvalues. VL < VU.
1271 : : !! Not referenced if RANGE = 'A' or 'I'.
1272 : : !!
1273 : : !! IL (input) INTEGER
1274 : : !! IU (input) INTEGER
1275 : : !! If RANGE='I', the indices (in ascending order) of the
1276 : : !! smallest and largest eigenvalues to be returned.
1277 : : !! 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
1278 : : !! Not referenced if RANGE = 'A' or 'V'.
1279 : : !!
1280 : : !! ABSTOL (input) REAL(DP)
1281 : : !! The absolute error tolerance for the eigenvalues.
1282 : : !! An approximate eigenvalue is accepted as converged
1283 : : !! when it is determined to lie in an interval [a,b]
1284 : : !! of width less than or equal to
1285 : : !!
1286 : : !! ABSTOL + EPS * max( |a|,|b| ) ,
1287 : : !!
1288 : : !! where EPS is the machine precision. If ABSTOL is less than
1289 : : !! or equal to zero, then EPS*|T| will be used in its place,
1290 : : !! where |T| is the 1-norm of the tridiagonal matrix obtained
1291 : : !! by reducing A to tridiagonal form.
1292 : : !!
1293 : : !! Eigenvalues will be computed most accurately when ABSTOL is
1294 : : !! set to twice the underflow threshold 2*DLAMCH('S'), not zero.
1295 : : !! If this routine returns with INFO>0, indicating that some
1296 : : !! eigenvectors did not converge, try setting ABSTOL to
1297 : : !! 2*DLAMCH('S').
1298 : : !!
1299 : : !! See "Computing Small Singular Values of Bidiagonal Matrices
1300 : : !! with Guaranteed High Relative Accuracy," by Demmel and
1301 : : !! Kahan, LAPACK Working Note #3.
1302 : : !!
1303 : : !! LDZ (input) INTEGER
1304 : : !! The leading dimension of the array Z. LDZ >= 1, and if
1305 : : !! JOBZ = 'V', LDZ >= max(1,N).
1306 : : !!
1307 : : !! [comm]=MPI communicator for ScaLAPACK inversion. Only available if the code has been compiled with Scalapack support.
1308 : : !! To avoid wasting CPU time the scalapack initialization is avoided if the number of processors in 1,
1309 : : !! in this case the sequential LAPACK routine is called.
1310 : : !!
1311 : : !! OUTPUT
1312 : : !! M (output) INTEGER
1313 : : !! The total number of eigenvalues found. 0 <= M <= N.
1314 : : !! If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
1315 : : !!
1316 : : !! W (output) REAL(DP) array, dimension (N)
1317 : : !! On normal exit, the first M elements contain the selected
1318 : : !! eigenvalues in ascending order.
1319 : : !!
1320 : : !! Z (output) COMPLEX(DP) array, dimension (LDZ, max(1,M))
1321 : : !! If JOBZ = 'V', then if INFO = 0, the first M columns of Z
1322 : : !! contain the orthonormal eigenvectors of the matrix A
1323 : : !! corresponding to the selected eigenvalues, with the i-th
1324 : : !! column of Z holding the eigenvector associated with W(i).
1325 : : !! If an eigenvector fails to converge, then that column of Z
1326 : : !! contains the latest approximation to the eigenvector, and the
1327 : : !! index of the eigenvector is returned in IFAIL.
1328 : : !! If JOBZ = 'N', then Z is not referenced.
1329 : : !! Note: the user must ensure that at least max(1,M) columns are
1330 : : !! supplied in the array Z; if RANGE = 'V', the exact value of M
1331 : : !! is not known in advance and an upper bound must be used.
1332 : : !!
1333 : : !! See also SIDE EFFECTS
1334 : : !!
1335 : : !! SIDE EFFECTS
1336 : : !! A (input/output) COMPLEX(DP) array, dimension (N, N)
1337 : : !! On entry, the Hermitian matrix A. If UPLO = 'U', the
1338 : : !! leading N-by-N upper triangular part of A contains the
1339 : : !! upper triangular part of the matrix A. If UPLO = 'L',
1340 : : !! the leading N-by-N lower triangular part of A contains
1341 : : !! the lower triangular part of the matrix A.
1342 : : !! On exit, the lower triangle (if UPLO='L') or the upper
1343 : : !! triangle (if UPLO='U') of A, including the diagonal, is
1344 : : !! destroyed.
1345 : : !!
1346 : : !! SOURCE
1347 : :
1348 : 1 : subroutine wrap_ZHEEVX(jobz,range,uplo,n,a,vl,vu,il,iu,abstol,m,w,z,ldz,comm)
1349 : :
1350 : : !Arguments ------------------------------------
1351 : : !scalars
1352 : : integer,intent(in) :: il,iu,ldz,n
1353 : : integer,optional,intent(in) :: comm
1354 : : integer,intent(inout) :: m
1355 : : real(dp),intent(in) :: abstol,vl,vu
1356 : : character(len=*),intent(in) :: jobz,range,uplo
1357 : : !arrays
1358 : : real(dp),intent(out) :: w(n)
1359 : : complex(dp),intent(out) :: z(ldz,m)
1360 : : complex(dp),intent(inout) :: a(n,n)
1361 : :
1362 : : !Local variables ------------------------------
1363 : : !scalars
1364 : : integer :: lwork,info,nprocs
1365 : : logical :: use_scalapack
1366 : : character(len=500) :: msg
1367 : : !arrays
1368 : 1 : integer,allocatable :: ifail(:),iwork(:)
1369 : 1 : real(dp),allocatable :: rwork(:)
1370 : 1 : complex(dp),allocatable :: work(:)
1371 : : #ifdef HAVE_LINALG_SCALAPACK
1372 : : integer :: ierr,istwf_k
1373 : : logical :: want_eigenvectors
1374 : : type(slkmat_dp_t) :: Slk_mat,Slk_vec
1375 : : type(slk_processor_t) :: Slk_processor
1376 : : #endif
1377 : : !************************************************************************
1378 : :
1379 : 1 : use_scalapack=.FALSE.
1380 [ - + ]: 1 : if (PRESENT(comm)) then
1381 : 0 : nprocs = xmpi_comm_size(comm)
1382 : : #ifdef HAVE_LINALG_SCALAPACK
1383 : : use_scalapack = (nprocs>1)
1384 : : #endif
1385 : : end if
1386 : :
1387 : : SELECT CASE(use_scalapack)
1388 : : CASE (.FALSE.) ! Standard LAPACK call.
1389 : :
1390 : 1 : lwork = MAX(1,2*n)
1391 [ - + + - : 3 : ABI_MALLOC(work,(lwork))
+ - - + -
+ ]
1392 [ + - - + : 3 : ABI_MALLOC(rwork,(7*n))
+ - + - +
- - + -
+ ]
1393 [ + - - + : 3 : ABI_MALLOC(iwork,(5*n))
+ - + - +
- - + -
+ ]
1394 [ - + + - : 3 : ABI_MALLOC(ifail,(n))
+ - + - -
+ - + ]
1395 : :
1396 : 1 : call ZHEEVX(jobz,range,uplo,n,a,n,vl,vu,il,iu,abstol,m,w,z,ldz,work,lwork,rwork,iwork,ifail,info)
1397 : :
1398 [ - + ]: 1 : if (info < 0) then
1399 : 0 : write(msg,'(a,i0,a)')" The ",-info,"-th argument of ZHEEVX had an illegal value."
1400 : 0 : ABI_ERROR(msg)
1401 : : end if
1402 : :
1403 [ - + ]: 1 : if (info > 0) then
1404 : 0 : write(msg,'(2a,i0,a)')"ZHEEVX: the algorithm failed to converge; ",ch10,&
1405 : 0 : info,"eigenvectors failed to converge. "
1406 : 0 : ABI_ERROR(msg)
1407 : : end if
1408 : :
1409 [ - + ]: 1 : ABI_FREE(iwork)
1410 [ - + ]: 1 : ABI_FREE(ifail)
1411 [ - + ]: 1 : ABI_FREE(rwork)
1412 [ - + ]: 1 : ABI_FREE(work)
1413 : : RETURN
1414 : :
1415 : : CASE (.TRUE.)
1416 : :
1417 : : #ifdef HAVE_LINALG_SCALAPACK
1418 : : call Slk_processor%init(comm)
1419 : : istwf_k=1
1420 : :
1421 : : ! Initialize and fill Scalapack matrix from the global one.
1422 : : call Slk_mat%init(n,n,Slk_processor,istwf_k)
1423 : : call slk_matrix_from_global_dpc_2D(Slk_mat,uplo,a)
1424 : :
1425 : : want_eigenvectors = firstchar(jobz,(/"V","v"/))
1426 : : if (want_eigenvectors) then
1427 : : ! Initialize the distributed vectors.
1428 : : call Slk_vec%init(n,n,Slk_processor,istwf_k)
1429 : : end if
1430 : :
1431 : : ! Solve the problem.
1432 : : call slk_mat%pzheevx(jobz,range,uplo,vl,vu,il,iu,abstol,Slk_vec,m,w)
1433 : : call Slk_mat%free()
1434 : :
1435 : : if (want_eigenvectors) then ! A is overwritten with the eigenvectors
1436 : : z = czero
1437 : : call slk_matrix_to_global_dpc_2D(Slk_vec,"All",z) ! Fill the entries calculated by this node.
1438 : : call Slk_vec%free()
1439 : : call xmpi_sum(z,comm,ierr) ! Fill the remaining entries of the global matrix
1440 : : end if
1441 : :
1442 : : call Slk_processor%free()
1443 : : RETURN
1444 : : #endif
1445 : :
1446 : : ABI_BUG("You should not be here!")
1447 : : END SELECT
1448 : :
1449 : : end subroutine wrap_ZHEEVX
1450 : : !!***
1451 : :
1452 : : !----------------------------------------------------------------------
1453 : :
1454 : : !!****f* m_hide_lapack/xheevx_cplex
1455 : : !! NAME
1456 : : !! xheevx_cplex
1457 : : !!
1458 : : !! FUNCTION
1459 : : !! xheevx_cplex computes selected eigenvalues and, optionally, eigenvectors
1460 : : !! of a (real symmetric|complex Hermitian) matrix A. Eigenvalues and eigenvectors can
1461 : : !! be selected by specifying either a range of values or a range of
1462 : : !! indices for the desired eigenvalues.
1463 : : !!
1464 : : !! INPUTS
1465 : : !! JOBZ (input) CHARACTER*1
1466 : : !! = 'N': Compute eigenvalues only;
1467 : : !! = 'V': Compute eigenvalues and eigenvectors.
1468 : : !!
1469 : : !! RANGE (input) CHARACTER*1
1470 : : !! = 'A': all eigenvalues will be found.
1471 : : !! = 'V': all eigenvalues in the half-open interval (VL,VU]
1472 : : !! will be found.
1473 : : !! = 'I': the IL-th through IU-th eigenvalues will be found.
1474 : : !!
1475 : : !! UPLO (input) CHARACTER*1
1476 : : !! = 'U': Upper triangle of A is stored;
1477 : : !! = 'L': Lower triangle of A is stored.
1478 : : !!
1479 : : !! CPLEX Size of the first dimension of the matrix A.
1480 : : !! 1 for real symmetric matrix
1481 : : !! 2 for complex Hermitian matrix.
1482 : : !!
1483 : : !! N (input) INTEGER
1484 : : !! The order of the matrix A. N >= 0.
1485 : : !!
1486 : : !! LDA (input) INTEGER
1487 : : !! The leading dimension of the array A. LDA >= max(1,N).
1488 : : !!
1489 : : !! VL (input) REAL(DP)
1490 : : !! VU (input) REAL(DP)
1491 : : !! If RANGE='V', the lower and upper bounds of the interval to
1492 : : !! be searched for eigenvalues. VL < VU.
1493 : : !! Not referenced if RANGE = 'A' or 'I'.
1494 : : !!
1495 : : !! IL (input) INTEGER
1496 : : !! IU (input) INTEGER
1497 : : !! If RANGE='I', the indices (in ascending order) of the
1498 : : !! smallest and largest eigenvalues to be returned.
1499 : : !! 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
1500 : : !! Not referenced if RANGE = 'A' or 'V'.
1501 : : !!
1502 : : !! ABSTOL (input) REAL(DP)
1503 : : !! The absolute error tolerance for the eigenvalues.
1504 : : !! An approximate eigenvalue is accepted as converged
1505 : : !! when it is determined to lie in an interval [a,b]
1506 : : !! of width less than or equal to
1507 : : !!
1508 : : !! ABSTOL + EPS * max( |a|,|b| ) ,
1509 : : !!
1510 : : !! where EPS is the machine precision. If ABSTOL is less than
1511 : : !! or equal to zero, then EPS*|T| will be used in its place,
1512 : : !! where |T| is the 1-norm of the tridiagonal matrix obtained
1513 : : !! by reducing A to tridiagonal form.
1514 : : !!
1515 : : !! Eigenvalues will be computed most accurately when ABSTOL is
1516 : : !! set to twice the underflow threshold 2*DLAMCH('S'), not zero.
1517 : : !! If this routine returns with INFO>0, indicating that some
1518 : : !! eigenvectors did not converge, try setting ABSTOL to
1519 : : !! 2*DLAMCH('S').
1520 : : !!
1521 : : !! See "Computing Small Singular Values of Bidiagonal Matrices
1522 : : !! with Guaranteed High Relative Accuracy," by Demmel and
1523 : : !! Kahan, LAPACK Working Note #3.
1524 : : !!
1525 : : !! LDZ (input) INTEGER
1526 : : !! The leading dimension of the array Z. LDZ >= 1, and if
1527 : : !! JOBZ = 'V', LDZ >= max(1,N).
1528 : : !!
1529 : : !! [comm]=MPI communicator for ScaLAPACK inversion. Only available if the code has been compiled with Scalapack support.
1530 : : !! To avoid wasting CPU time the scalapack initialization is avoided if the number of processors in 1,
1531 : : !! in this case the sequential LAPACK routine is called.
1532 : : !!
1533 : : !! OUTPUT
1534 : : !! M (output) INTEGER
1535 : : !! The total number of eigenvalues found. 0 <= M <= N.
1536 : : !! If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
1537 : : !!
1538 : : !! W (output) REAL(DP) array, dimension (N)
1539 : : !! On normal exit, the first M elements contain the selected
1540 : : !! eigenvalues in ascending order.
1541 : : !!
1542 : : !! Z (output) REAL(DP) array, dimension (CPLEX, LDZ, max(1,M))
1543 : : !! If JOBZ = 'V', then if INFO = 0, the first M columns of Z
1544 : : !! contain the orthonormal eigenvectors of the matrix A
1545 : : !! corresponding to the selected eigenvalues, with the i-th
1546 : : !! column of Z holding the eigenvector associated with W(i).
1547 : : !! If an eigenvector fails to converge, then that column of Z
1548 : : !! contains the latest approximation to the eigenvector, and the
1549 : : !! index of the eigenvector is returned in IFAIL.
1550 : : !! If JOBZ = 'N', then Z is not referenced.
1551 : : !! Note: the user must ensure that at least max(1,M) columns are
1552 : : !! supplied in the array Z; if RANGE = 'V', the exact value of M
1553 : : !! is not known in advance and an upper bound must be used.
1554 : : !!
1555 : : !! See also SIDE EFFECTS
1556 : : !!
1557 : : !! SIDE EFFECTS
1558 : : !! A (input/output) REAL(DP) array, dimension (CPLEX, N, N)
1559 : : !! On entry, the (real symmetric|complex Hermitian) matrix A. If UPLO = 'U', the
1560 : : !! leading N-by-N upper triangular part of A contains the
1561 : : !! upper triangular part of the matrix A. If UPLO = 'L',
1562 : : !! the leading N-by-N lower triangular part of A contains
1563 : : !! the lower triangular part of the matrix A.
1564 : : !! On exit, the lower triangle (if UPLO='L') or the upper
1565 : : !! triangle (if UPLO='U') of A, including the diagonal, is
1566 : : !! destroyed.
1567 : : !!
1568 : : !! SOURCE
1569 : :
1570 : 10 : subroutine xheevx_cplex(jobz, range, uplo, cplex, n, a, vl, vu, il, iu, &
1571 : 10 : abstol, m, w, z, ldz, msg, ierr, comm)
1572 : :
1573 : : !Arguments ------------------------------------
1574 : : !scalars
1575 : : integer,intent(in) :: il,iu,ldz,n,cplex
1576 : : integer,optional,intent(in) :: comm
1577 : : integer,intent(inout) :: m
1578 : : integer,intent(out) :: ierr
1579 : : real(dp),intent(in) :: abstol,vl,vu
1580 : : character(len=*),intent(in) :: jobz,range,uplo
1581 : : character(len=*),intent(out) :: msg
1582 : : !arrays
1583 : : real(dp),intent(out) :: w(n)
1584 : : !real(dp),intent(out) :: z(cplex,ldz,n)
1585 : : real(dp),intent(out) :: z(cplex,ldz,m)
1586 : : real(dp),intent(inout) :: a(cplex,n,n)
1587 : :
1588 : : !Local variables ------------------------------
1589 : : !scalars
1590 : : integer :: lwork,nprocs
1591 : : logical :: use_scalapack
1592 : : !arrays
1593 : 10 : integer,allocatable :: ifail(:),iwork(:)
1594 : 10 : real(dp),allocatable :: rwork(:)
1595 : 10 : real(dp),allocatable :: work_real(:)
1596 : 10 : complex(dp),allocatable :: work_cplx(:)
1597 : : #ifdef HAVE_LINALG_SCALAPACK
1598 : : !integer :: istwf_k
1599 : : !logical :: want_eigenvectors
1600 : : !type(slkmat_dp_t) :: Slk_mat,Slk_vec
1601 : : !type(slk_processor_t) :: Slk_processor
1602 : : #endif
1603 : : !************************************************************************
1604 : :
1605 : 10 : use_scalapack=.FALSE.
1606 [ - + ]: 10 : if (PRESENT(comm)) then
1607 : 0 : nprocs = xmpi_comm_size(comm)
1608 : : #ifdef HAVE_LINALG_SCALAPACK
1609 : : use_scalapack = (nprocs>1)
1610 : : #endif
1611 : : end if
1612 : :
1613 [ - + ]: 10 : if (ALL(cplex/=(/1,2/))) then
1614 : 0 : write(msg,'(a,i0)')" Wrong value for cplex: ",cplex
1615 : 0 : ierr = 1; return
1616 : : end if
1617 : :
1618 : : SELECT CASE(use_scalapack)
1619 : : CASE (.FALSE.)
1620 : : ! Standard LAPACK call.
1621 : :
1622 [ - + ]: 10 : if (cplex==1) then
1623 : : ! Real symmetric case
1624 : 0 : lwork = MAX(1,8*n)
1625 [ # # # # : 0 : ABI_MALLOC(work_real,(lwork))
# # # # #
# ]
1626 [ # # # # : 0 : ABI_MALLOC(iwork,(5*n))
# # # # #
# # # #
# ]
1627 [ # # # # : 0 : ABI_MALLOC(ifail,(n))
# # # # #
# # # ]
1628 : :
1629 : 0 : call DSYEVX(jobz,range,uplo,n,a,n,vl,vu,il,iu,abstol,m,w,z,ldz,work_real,lwork,iwork,ifail,ierr)
1630 : :
1631 [ # # ]: 0 : if (ierr < 0) then
1632 : 0 : write(msg,'(a,i0,a)')" The ",-ierr,"-th argument of DSYEVX had an illegal value."
1633 : : end if
1634 : :
1635 [ # # ]: 0 : if (ierr > 0) then
1636 : : write(msg,'(2a,i0,a)')&
1637 : 0 : "DSYEVX: the algorithm failed to converge; ",ch10,ierr,"eigenvectors failed to converge. "
1638 : : end if
1639 : :
1640 [ # # ]: 0 : ABI_FREE(work_real)
1641 [ # # ]: 0 : ABI_FREE(iwork)
1642 [ # # ]: 0 : ABI_FREE(ifail)
1643 : 0 : RETURN
1644 : :
1645 : : else
1646 : : ! Complex Hermitian case.
1647 : 10 : lwork = MAX(1,2*n)
1648 [ - + + - : 30 : ABI_MALLOC(work_cplx,(lwork))
+ - - + -
+ ]
1649 [ + - - + : 30 : ABI_MALLOC(rwork,(7*n))
+ - + - +
- - + -
+ ]
1650 [ + - - + : 30 : ABI_MALLOC(iwork,(5*n))
+ - + - +
- - + -
+ ]
1651 [ - + + - : 30 : ABI_MALLOC(ifail,(n))
+ - + - -
+ - + ]
1652 : :
1653 : 10 : call ZHEEVX(jobz,range,uplo,n,a,n,vl,vu,il,iu,abstol,m,w,z,ldz,work_cplx,lwork,rwork,iwork,ifail,ierr)
1654 : :
1655 [ - + ]: 10 : if (ierr < 0) then
1656 : 0 : write(msg,'(a,i0,a)')" The ",-ierr,"-th argument of ZHEEVX had an illegal value."
1657 : : end if
1658 : :
1659 [ - + ]: 10 : if (ierr > 0) then
1660 : : write(msg,'(2a,i0,a)')&
1661 : 0 : "ZHEEVX: the algorithm failed to converge; ",ch10,ierr,"eigenvectors failed to converge. "
1662 : : end if
1663 : :
1664 [ - + ]: 10 : ABI_FREE(iwork)
1665 [ - + ]: 10 : ABI_FREE(ifail)
1666 [ - + ]: 10 : ABI_FREE(rwork)
1667 [ - + ]: 10 : ABI_FREE(work_cplx)
1668 : 10 : RETURN
1669 : : end if
1670 : :
1671 : : CASE (.TRUE.)
1672 : :
1673 : : #ifdef HAVE_LINALG_SCALAPACK
1674 : : ABI_ERROR("Not coded yet")
1675 : : ! call Slk_processor%init(comm)
1676 : : ! istwf_k=1
1677 : :
1678 : : ! ! Initialize and fill Scalapack matrix from the global one.
1679 : : ! call Slk_mat%init(n,n,Slk_processor,istwf_k)
1680 : : ! call slk_matrix_from_global_dpc_2D(Slk_mat,uplo,a)
1681 : :
1682 : : ! want_eigenvectors = firstchar(jobz,(/"V","v"/))
1683 : : ! if (want_eigenvectors) then ! Initialize the distributed vectors.
1684 : : ! call Slk_vec%init(n,n,Slk_processor,istwf_k)
1685 : : ! end if
1686 : :
1687 : : ! ! Solve the problem.
1688 : : ! call slk_mat%pzheevx(jobz,range,uplo,vl,vu,il,iu,abstol,Slk_vec,m,w)
1689 : : ! call Slk_mat%free()
1690 : : !
1691 : : ! if (want_eigenvectors) then ! A is overwritten with the eigenvectors
1692 : : ! z = czero
1693 : : ! call slk_matrix_to_global_dpc_2D(Slk_vec,"All",z) ! Fill the entries calculated by this node.
1694 : : ! call Slk_vec%free()
1695 : : ! call xmpi_sum(z,comm,ierr) ! Fill the remaining entries of the global matrix
1696 : : ! end if
1697 : :
1698 : : ! call Slk_processor%free()
1699 : :
1700 : : RETURN
1701 : : #endif
1702 : :
1703 : : ABI_BUG("You should not be here!")
1704 : : END SELECT
1705 : :
1706 : : end subroutine xheevx_cplex
1707 : : !!***
1708 : :
1709 : : !----------------------------------------------------------------------
1710 : :
1711 : : !!****f* m_hide_lapack/wrap_ZHEGVX
1712 : : !! NAME
1713 : : !! wrap_ZHEGVX
1714 : : !!
1715 : : !! FUNCTION
1716 : : !! wrap_ZHEGVX - compute selected eigenvalues, and optionally, eigenvectors of a
1717 : : !! complex generalized Hermitian-definite eigenproblem, of the form A*x=(lambda)*B*x, A*Bx=(lambda)*x, or B*A*x=(lambda)*x.
1718 : : !! Eigenvalues and eigenvectors can be selected by specifying either a range of values or a range of
1719 : : !! indices for the desired eigenvalues.
1720 : : !!
1721 : : !! INPUTS
1722 : : !!
1723 : : !! ITYPE (input) INTEGER Specifies the problem type to be solved:
1724 : : !! = 1: A*x = (lambda)*B*x
1725 : : !! = 2: A*B*x = (lambda)*x
1726 : : !! = 3: B*A*x = (lambda)*x
1727 : : !!
1728 : : !! JOBZ (input) CHARACTER*1
1729 : : !! = 'N': Compute eigenvalues only;
1730 : : !! = 'V': Compute eigenvalues and eigenvectors.
1731 : : !!
1732 : : !! RANGE (input) CHARACTER*1
1733 : : !! = 'A': all eigenvalues will be found.
1734 : : !! = 'V': all eigenvalues in the half-open interval (VL,VU]
1735 : : !! will be found.
1736 : : !! = 'I': the IL-th through IU-th eigenvalues will be found.
1737 : : !!
1738 : : !! UPLO (input) CHARACTER*1
1739 : : !! = 'U': Upper triangle of A is stored;
1740 : : !! = 'L': Lower triangle of A is stored.
1741 : : !!
1742 : : !! N (input) INTEGER
1743 : : !! The order of the matrices A and B. N >= 0.
1744 : : !!
1745 : : !! LDA (input) INTEGER
1746 : : !! The leading dimension of the array A. LDA >= max(1,N).
1747 : : !!
1748 : : !! VL (input) REAL(DP)
1749 : : !! VU (input) REAL(DP)
1750 : : !! If RANGE='V', the lower and upper bounds of the interval to
1751 : : !! be searched for eigenvalues. VL < VU.
1752 : : !! Not referenced if RANGE = 'A' or 'I'.
1753 : : !!
1754 : : !! IL (input) INTEGER
1755 : : !! IU (input) INTEGER
1756 : : !! If RANGE='I', the indices (in ascending order) of the
1757 : : !! smallest and largest eigenvalues to be returned.
1758 : : !! 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
1759 : : !! Not referenced if RANGE = 'A' or 'V'.
1760 : : !!
1761 : : !! ABSTOL (input) REAL(DP)
1762 : : !! The absolute error tolerance for the eigenvalues.
1763 : : !! An approximate eigenvalue is accepted as converged
1764 : : !! when it is determined to lie in an interval [a,b]
1765 : : !! of width less than or equal to
1766 : : !!
1767 : : !! ABSTOL + EPS * max( |a|,|b| ) ,
1768 : : !!
1769 : : !! where EPS is the machine precision. If ABSTOL is less than
1770 : : !! or equal to zero, then EPS*|T| will be used in its place,
1771 : : !! where |T| is the 1-norm of the tridiagonal matrix obtained
1772 : : !! by reducing A to tridiagonal form.
1773 : : !!
1774 : : !! Eigenvalues will be computed most accurately when ABSTOL is
1775 : : !! set to twice the underflow threshold 2*DLAMCH('S'), not zero.
1776 : : !! If this routine returns with INFO>0, indicating that some
1777 : : !! eigenvectors did not converge, try setting ABSTOL to
1778 : : !! 2*DLAMCH('S').
1779 : : !!
1780 : : !! LDZ (input) INTEGER
1781 : : !! The leading dimension of the array Z. LDZ >= 1, and if
1782 : : !! JOBZ = 'V', LDZ >= max(1,N).
1783 : : !!
1784 : : !! [comm]=MPI communicator for ScaLAPACK inversion. Only available if the code has been compiled with Scalapack support.
1785 : : !! To avoid wasting CPU time the scalapack initialization is avoided if the number of processors in 1,
1786 : : !! in this case the sequential LAPACK routine is called.
1787 : : !!
1788 : : !! OUTPUT
1789 : : !! M (output) INTEGER
1790 : : !! The total number of eigenvalues found. 0 <= M <= N.
1791 : : !! If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
1792 : : !!
1793 : : !! W (output) REAL(DP) array, dimension (N)
1794 : : !! On normal exit, the first M elements contain the selected
1795 : : !! eigenvalues in ascending order.
1796 : : !!
1797 : : !! Z (output) COMPLEX(DP) array, dimension (LDZ, max(1,M))
1798 : : !! If JOBZ = 'V', then if INFO = 0, the first M columns of Z
1799 : : !! contain the orthonormal eigenvectors of the matrix A
1800 : : !! corresponding to the selected eigenvalues, with the i-th
1801 : : !! column of Z holding the eigenvector associated with W(i).
1802 : : !! 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.
1803 : : !! If an eigenvector fails to converge, then that column of Z
1804 : : !! contains the latest approximation to the eigenvector, and the
1805 : : !! index of the eigenvector is returned in IFAIL.
1806 : : !! If JOBZ = 'N', then Z is not referenced.
1807 : : !! Note: the user must ensure that at least max(1,M) columns are
1808 : : !! supplied in the array Z; if RANGE = 'V', the exact value of M
1809 : : !! is not known in advance and an upper bound must be used.
1810 : : !!
1811 : : !! See also SIDE EFFECTS
1812 : : !!
1813 : : !! SIDE EFFECTS
1814 : : !! A (input/output) COMPLEX(DP) array, dimension (N, N)
1815 : : !! On entry, the Hermitian matrix A. If UPLO = 'U', the
1816 : : !! leading N-by-N upper triangular part of A contains the
1817 : : !! upper triangular part of the matrix A. If UPLO = "L",
1818 : : !! the leading N-by-N lower triangular part of A contains
1819 : : !! the lower triangular part of the matrix A.
1820 : : !!
1821 : : !! On exit, the lower triangle (if UPLO="L") or the upper
1822 : : !! triangle (if UPLO="U") of A, including the diagonal, is
1823 : : !! destroyed.
1824 : : !!
1825 : : !! B (input/output) COMPLEX(DP) array, dimension (LDB, N)
1826 : : !! On entry, the Hermitian matrix B. If UPLO = "U", the leading N-by-N upper triangular part
1827 : : !! of B contains the upper triangular part of the matrix B.
1828 : : !! If UPLO = "L", the leading N-by-N lower triangular part of B contains the lower triangular part of the matrix B.
1829 : : !!
1830 : : !! On exit, if INFO <= N, the part of B containing the matrix is overwritten by the triangular factor
1831 : : !! U or L from the Cholesky factorization B = U**H*U or B = L*L**H.
1832 : : !!
1833 : : !! SOURCE
1834 : :
1835 : 0 : subroutine wrap_ZHEGVX(itype,jobz,range,uplo,n,a,b,vl,vu,il,iu,abstol,m,w,z,ldz,comm)
1836 : :
1837 : : !Arguments ------------------------------------
1838 : : !scalars
1839 : : integer,intent(in) :: il,iu,ldz,n,itype
1840 : : integer,optional,intent(in) :: comm
1841 : : integer,intent(inout) :: m
1842 : : real(dp),intent(in) :: abstol,vl,vu
1843 : : character(len=*),intent(in) :: jobz,range,uplo
1844 : : !arrays
1845 : : real(dp),intent(out) :: w(n)
1846 : : !complex(dp),intent(out) :: z(ldz,n)
1847 : : complex(dp),intent(out) :: z(ldz,m)
1848 : : complex(dp),intent(inout) :: a(n,n),b(n,n)
1849 : :
1850 : : !Local variables ------------------------------
1851 : : !scalars
1852 : : integer :: lwork,info,nprocs,ii
1853 : : logical :: use_scalapack
1854 : : character(len=500) :: msg
1855 : : !arrays
1856 : 0 : integer,allocatable :: ifail(:),iwork(:)
1857 : 0 : real(dp),allocatable :: rwork(:)
1858 : 0 : complex(dp),allocatable :: work(:)
1859 : : #ifdef HAVE_LINALG_SCALAPACK
1860 : : integer :: ierr,istwf_k
1861 : : logical :: want_eigenvectors
1862 : : type(slkmat_dp_t) :: Slk_matA,Slk_matB,Slk_vec
1863 : : type(slk_processor_t) :: Slk_processor
1864 : : #endif
1865 : : !************************************************************************
1866 : :
1867 : 0 : use_scalapack=.FALSE.
1868 [ # # ]: 0 : if (PRESENT(comm)) then
1869 : 0 : nprocs = xmpi_comm_size(comm)
1870 : : #ifdef HAVE_LINALG_SCALAPACK
1871 : : use_scalapack = (nprocs>1)
1872 : : #endif
1873 : : end if
1874 : :
1875 : : SELECT CASE(use_scalapack)
1876 : : CASE (.FALSE.)
1877 : : ! Standard LAPACK call.
1878 : 0 : lwork = MAX(1,2*n)
1879 [ # # # # : 0 : ABI_MALLOC(work,(lwork))
# # # # #
# ]
1880 [ # # # # : 0 : ABI_MALLOC(rwork,(7*n))
# # # # #
# # # #
# ]
1881 [ # # # # : 0 : ABI_MALLOC(iwork,(5*n))
# # # # #
# # # #
# ]
1882 [ # # # # : 0 : ABI_MALLOC(ifail,(n))
# # # # #
# # # ]
1883 : :
1884 : 0 : 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)
1885 : :
1886 [ # # ]: 0 : if (info < 0) then
1887 : 0 : write(msg,'(a,i0,a)')" The ",-info,"-th argument of ZHEGVX had an illegal value."
1888 : 0 : ABI_ERROR(msg)
1889 : : end if
1890 : :
1891 [ # # ]: 0 : if (info > 0) then
1892 [ # # ]: 0 : if (info<= n) then
1893 : 0 : write(msg,'(a,i0,a)')"ZHEGVX failed to converge: ",info," eigenvectors failed to converge. "
1894 : : else
1895 : 0 : ii = info -n
1896 : : write(msg,'(3a,i0,3a)')&
1897 : 0 : "ZHEEVX failed to converge: ",ch10,&
1898 : 0 : "The leading minor of order ",ii," of B is not positive definite. ",ch10,&
1899 : 0 : "The factorization of B could not be completed and no eigenvalues or eigenvectors were computed."
1900 : : end if
1901 : 0 : ABI_ERROR(msg)
1902 : : end if
1903 : :
1904 [ # # ]: 0 : ABI_FREE(iwork)
1905 [ # # ]: 0 : ABI_FREE(ifail)
1906 [ # # ]: 0 : ABI_FREE(rwork)
1907 [ # # ]: 0 : ABI_FREE(work)
1908 : : RETURN
1909 : :
1910 : : CASE (.TRUE.)
1911 : :
1912 : : #ifdef HAVE_LINALG_SCALAPACK
1913 : : call Slk_processor%init(comm)
1914 : : istwf_k=1
1915 : :
1916 : : ! Initialize and fill Scalapack matrix from the global one.
1917 : : call Slk_matA%init(n,n,Slk_processor,istwf_k)
1918 : : call slk_matrix_from_global_dpc_2D(Slk_matA,uplo,a)
1919 : :
1920 : : call Slk_matB%init(n,n,Slk_processor,istwf_k)
1921 : : call slk_matrix_from_global_dpc_2D(Slk_matB,uplo,b)
1922 : :
1923 : : want_eigenvectors = firstchar(jobz,(/"V","v"/))
1924 : : if (want_eigenvectors) then ! Initialize the distributed vectors.
1925 : : call Slk_vec%init(n,n,Slk_processor,istwf_k)
1926 : : end if
1927 : :
1928 : : ! Solve the problem.
1929 : : ABI_ERROR("slk_pZHEGVX not coded yet")
1930 : : ! TODO write the scaLAPACK wrapper.
1931 : : !call slk_pZHEGVX(itype,jobz,range,uplo,Slk_matA,Slk_matB,vl,vu,il,iu,abstol,Slk_vec,m,w)
1932 : :
1933 : : call Slk_matA%free()
1934 : : call Slk_matB%free()
1935 : :
1936 : : if (want_eigenvectors) then ! A is overwritten with the eigenvectors
1937 : : z = czero
1938 : : call slk_matrix_to_global_dpc_2D(Slk_vec,"All",z) ! Fill the entries calculated by this node.
1939 : : call Slk_vec%free()
1940 : : call xmpi_sum(z,comm,ierr) ! Fill the remaining entries of the global matrix
1941 : : end if
1942 : :
1943 : : call Slk_processor%free()
1944 : :
1945 : : RETURN
1946 : : #endif
1947 : :
1948 : : ABI_BUG("You should not be here!")
1949 : : END SELECT
1950 : :
1951 : : end subroutine wrap_ZHEGVX
1952 : : !!***
1953 : :
1954 : : !----------------------------------------------------------------------
1955 : :
1956 : : !!****f* m_hide_lapack/xhegvx_cplex
1957 : : !! NAME
1958 : : !! xhegvx_cplex
1959 : : !!
1960 : : !! FUNCTION
1961 : : !! xhegvx_cplex - compute selected eigenvalues, and optionally, eigenvectors of a
1962 : : !! (real symmetric-definite|complex generalized Hermitian-definite) eigenproblem, of the form
1963 : : !! A*x=(lambda)*B*x, A*Bx=(lambda)*x, or B*A*x=(lambda)*x.
1964 : : !! Here A and B are assumed to be (real symmetric|complex Hermitian) and B is also positive definite.
1965 : : !! Eigenvalues and eigenvectors can be selected by specifying either a range of values or a range of
1966 : : !! indices for the desired eigenvalues.
1967 : : !!
1968 : : !! INPUTS
1969 : : !!
1970 : : !! ITYPE (input) INTEGER Specifies the problem type to be solved:
1971 : : !! = 1: A*x = (lambda)*B*x
1972 : : !! = 2: A*B*x = (lambda)*x
1973 : : !! = 3: B*A*x = (lambda)*x
1974 : : !!
1975 : : !! JOBZ (input) CHARACTER*1
1976 : : !! = 'N': Compute eigenvalues only;
1977 : : !! = 'V': Compute eigenvalues and eigenvectors.
1978 : : !!
1979 : : !! RANGE (input) CHARACTER*1
1980 : : !! = 'A': all eigenvalues will be found.
1981 : : !! = 'V': all eigenvalues in the half-open interval (VL,VU]
1982 : : !! will be found.
1983 : : !! = 'I': the IL-th through IU-th eigenvalues will be found.
1984 : : !!
1985 : : !! UPLO (input) CHARACTER*1
1986 : : !! = 'U': Upper triangle of A is stored;
1987 : : !! = 'L': Lower triangle of A is stored.
1988 : : !!
1989 : : !! CPLEX Size of the first dimension of the matrices A and B
1990 : : !! 1 for Real symmetric matrices
1991 : : !! 2 for complex Hermitianmatrices
1992 : : !!
1993 : : !! N (input) INTEGER
1994 : : !! The order of the matrices A and B. N >= 0.
1995 : : !!
1996 : : !! LDA (input) INTEGER
1997 : : !! The leading dimension of the array A. LDA >= max(1,N).
1998 : : !!
1999 : : !! VL (input) REAL(DP)
2000 : : !! VU (input) REAL(DP)
2001 : : !! If RANGE='V', the lower and upper bounds of the interval to
2002 : : !! be searched for eigenvalues. VL < VU.
2003 : : !! Not referenced if RANGE = 'A' or 'I'.
2004 : : !!
2005 : : !! IL (input) INTEGER
2006 : : !! IU (input) INTEGER
2007 : : !! If RANGE='I', the indices (in ascending order) of the
2008 : : !! smallest and largest eigenvalues to be returned.
2009 : : !! 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
2010 : : !! Not referenced if RANGE = 'A' or 'V'.
2011 : : !!
2012 : : !! ABSTOL (input) REAL(DP)
2013 : : !! The absolute error tolerance for the eigenvalues.
2014 : : !! An approximate eigenvalue is accepted as converged
2015 : : !! when it is determined to lie in an interval [a,b]
2016 : : !! of width less than or equal to
2017 : : !!
2018 : : !! ABSTOL + EPS * max( |a|,|b| ) ,
2019 : : !!
2020 : : !! where EPS is the machine precision. If ABSTOL is less than
2021 : : !! or equal to zero, then EPS*|T| will be used in its place,
2022 : : !! where |T| is the 1-norm of the tridiagonal matrix obtained
2023 : : !! by reducing A to tridiagonal form.
2024 : : !!
2025 : : !! Eigenvalues will be computed most accurately when ABSTOL is
2026 : : !! set to twice the underflow threshold 2*DLAMCH('S'), not zero.
2027 : : !! If this routine returns with INFO>0, indicating that some
2028 : : !! eigenvectors did not converge, try setting ABSTOL to
2029 : : !! 2*DLAMCH('S').
2030 : : !!
2031 : : !! LDZ (input) INTEGER
2032 : : !! The leading dimension of the array Z. LDZ >= 1, and if
2033 : : !! JOBZ = 'V', LDZ >= max(1,N).
2034 : : !!
2035 : : !! [comm]=MPI communicator for ScaLAPACK inversion. Only available if the code has been compiled with Scalapack support.
2036 : : !! To avoid wasting CPU time the scalapack initialization is avoided if the number of processors in 1,
2037 : : !! in this case the sequential LAPACK routine is called.
2038 : : !!
2039 : : !! OUTPUT
2040 : : !! M (output) INTEGER
2041 : : !! The total number of eigenvalues found. 0 <= M <= N.
2042 : : !! If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
2043 : : !!
2044 : : !! W (output) REAL(DP) array, dimension (N)
2045 : : !! On normal exit, the first M elements contain the selected
2046 : : !! eigenvalues in ascending order.
2047 : : !!
2048 : : !! Z (output) REAL(DP) array, dimension (CPLEX ,LDZ, max(1,M))
2049 : : !! If JOBZ = 'V', then if INFO = 0, the first M columns of Z
2050 : : !! contain the orthonormal eigenvectors of the matrix A
2051 : : !! corresponding to the selected eigenvalues, with the i-th
2052 : : !! column of Z holding the eigenvector associated with W(i).
2053 : : !! The eigenvectors are normalized as follows:
2054 : : !! if ITYPE = 1 or 2, Z**T*B*Z = I; if ITYPE = 3, Z**T*inv(B)*Z = I.
2055 : : !!
2056 : : !! If an eigenvector fails to converge, then that column of Z
2057 : : !! contains the latest approximation to the eigenvector, and the
2058 : : !! index of the eigenvector is returned in IFAIL.
2059 : : !! If JOBZ = 'N', then Z is not referenced.
2060 : : !! Note: the user must ensure that at least max(1,M) columns are
2061 : : !! supplied in the array Z; if RANGE = 'V', the exact value of M
2062 : : !! is not known in advance and an upper bound must be used.
2063 : : !!
2064 : : !! See also SIDE EFFECTS
2065 : : !!
2066 : : !! SIDE EFFECTS
2067 : : !! A (input/output) REAL(DP) array, dimension (CPLEX, N, N)
2068 : : !! On entry, the (real symmetric| complex Hermitian) matrix A. If UPLO = 'U', the
2069 : : !! leading N-by-N upper triangular part of A contains the
2070 : : !! upper triangular part of the matrix A. If UPLO = "L",
2071 : : !! the leading N-by-N lower triangular part of A contains
2072 : : !! the lower triangular part of the matrix A.
2073 : : !!
2074 : : !! On exit, the lower triangle (if UPLO="L") or the upper
2075 : : !! triangle (if UPLO="U") of A, including the diagonal, is
2076 : : !! destroyed.
2077 : : !!
2078 : : !! B (input/output) REAL(DP) array, dimension (CPLEX, LDB, N)
2079 : : !! On entry, the (real symmetric| complex Hermitian) matrix B. If UPLO = "U", the leading N-by-N upper triangular part
2080 : : !! of B contains the upper triangular part of the matrix B.
2081 : : !! If UPLO = "L", the leading N-by-N lower triangular part of B contains the lower triangular part of the matrix B.
2082 : : !!
2083 : : !! On exit, if INFO <= N, the part of B containing the matrix is overwritten by the triangular factor
2084 : : !! U or L from the Cholesky factorization B = U**H*U or B = L*L**H.
2085 : : !!
2086 : : !! SOURCE
2087 : :
2088 : 0 : subroutine xhegvx_cplex(itype, jobz, range, uplo, cplex, n, a, b, &
2089 : 0 : vl, vu, il, iu, abstol, m, w, z, ldz, msg, ierr, comm)
2090 : :
2091 : : !Arguments ------------------------------------
2092 : : !scalars
2093 : : integer,intent(in) :: il,iu,ldz,n,itype,cplex
2094 : : integer,optional,intent(in) :: comm
2095 : : integer,intent(inout) :: m
2096 : : integer,intent(out) :: ierr
2097 : : real(dp),intent(in) :: abstol,vl,vu
2098 : : character(len=*),intent(in) :: jobz,range,uplo
2099 : : character(len=*),intent(out) :: msg
2100 : : !arrays
2101 : : real(dp),intent(out) :: w(n)
2102 : : !real(dp),intent(out) :: z(cplex,ldz,n)
2103 : : real(dp),intent(out) :: z(cplex,ldz,m)
2104 : : real(dp),intent(inout) :: a(cplex,n,n),b(cplex,n,n)
2105 : :
2106 : : !Local variables ------------------------------
2107 : : !scalars
2108 : : integer :: lwork,nprocs,ii
2109 : : logical :: use_scalapack
2110 : : !arrays
2111 : 0 : integer,allocatable :: ifail(:),iwork(:)
2112 : 0 : real(dp),allocatable :: rwork(:)
2113 : 0 : real(dp),allocatable :: work_real(:)
2114 : 0 : complex(dp),allocatable :: work_cplx(:)
2115 : : #ifdef HAVE_LINALG_SCALAPACK
2116 : : !integer :: istwf_k
2117 : : !logical :: want_eigenvectors
2118 : : !type(slkmat_dp_t) :: Slk_matA,Slk_matB,Slk_vec
2119 : : !type(slk_processor_t) :: Slk_processor
2120 : : #endif
2121 : : !************************************************************************
2122 : :
2123 : 0 : use_scalapack=.FALSE.
2124 [ # # ]: 0 : if (PRESENT(comm)) then
2125 : 0 : nprocs = xmpi_comm_size(comm)
2126 : : #ifdef HAVE_LINALG_SCALAPACK
2127 : : use_scalapack = (nprocs>1)
2128 : : #endif
2129 : : end if
2130 : :
2131 [ # # ]: 0 : if (ALL(cplex/=(/1,2/))) then
2132 : 0 : write(msg,'(a,i0)')" Wrong value for cplex: ",cplex
2133 : 0 : ierr = 1; return
2134 : : end if
2135 : :
2136 : : SELECT CASE(use_scalapack)
2137 : :
2138 : : CASE (.FALSE.)
2139 : : ! Standard LAPACK call.
2140 [ # # ]: 0 : if (cplex==1) then
2141 : : ! Real symmetric case
2142 : 0 : lwork = MAX(1,8*n)
2143 : :
2144 [ # # # # : 0 : ABI_MALLOC(work_real,(lwork))
# # # # #
# ]
2145 [ # # # # : 0 : ABI_MALLOC(iwork,(5*n))
# # # # #
# # # #
# ]
2146 [ # # # # : 0 : ABI_MALLOC(ifail,(n))
# # # # #
# # # ]
2147 : :
2148 : 0 : 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)
2149 : :
2150 [ # # ]: 0 : if (ierr < 0) then
2151 : 0 : write(msg,'(a,i0,a)')" The ",-ierr,"-th argument of DSYGVX had an illegal value."
2152 : : end if
2153 : :
2154 [ # # ]: 0 : if (ierr > 0) then
2155 [ # # ]: 0 : if (ierr<= n) then
2156 : 0 : write(msg,'(a,i0,a)')" DSYGVX failed to converge: ",ierr," eigenvectors failed to converge. "
2157 : : else
2158 : 0 : ii = ierr - n
2159 : : write(msg,'(3a,i0,3a)')&
2160 : 0 : " DSYGVX failed to converge: ",ch10,&
2161 : 0 : " The leading minor of order ",ii," of B is not positive definite. ",ch10,&
2162 : 0 : " The factorization of B could not be completed and no eigenvalues or eigenvectors were computed."
2163 : : end if
2164 : : end if
2165 : :
2166 [ # # ]: 0 : ABI_FREE(iwork)
2167 [ # # ]: 0 : ABI_FREE(ifail)
2168 [ # # ]: 0 : ABI_FREE(work_real)
2169 : 0 : RETURN
2170 : :
2171 : : else
2172 : : ! Complex Hermitian case.
2173 : 0 : lwork = MAX(1,2*n)
2174 : :
2175 [ # # # # : 0 : ABI_MALLOC(work_cplx,(lwork))
# # # # #
# ]
2176 [ # # # # : 0 : ABI_MALLOC(rwork,(7*n))
# # # # #
# # # #
# ]
2177 [ # # # # : 0 : ABI_MALLOC(iwork,(5*n))
# # # # #
# # # #
# ]
2178 [ # # # # : 0 : ABI_MALLOC(ifail,(n))
# # # # #
# # # ]
2179 : :
2180 : : !write(std_out,*)"Calling ZHEGVX"
2181 : 0 : 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)
2182 : :
2183 [ # # ]: 0 : if (ierr < 0) then
2184 : 0 : write(msg,'(a,i0,a)')"The ",-ierr,"-th argument of ZHEGVX had an illegal value."
2185 : : end if
2186 : :
2187 [ # # ]: 0 : if (ierr > 0) then
2188 [ # # ]: 0 : if (ierr<= n) then
2189 : 0 : write(msg,'(a,i0,a)')"ZHEGVX failed to converge: ",ierr," eigenvectors failed to converge. "
2190 : : else
2191 : 0 : ii = ierr -n
2192 : : write(msg,'(3a,i0,3a)')&
2193 : 0 : "ZHEEVX failed to converge: ",ch10,&
2194 : 0 : "The leading minor of order ",ii," of B is not positive definite. ",ch10,&
2195 : 0 : "The factorization of B could not be completed and no eigenvalues or eigenvectors were computed."
2196 : : end if
2197 : : end if
2198 : :
2199 [ # # ]: 0 : ABI_FREE(iwork)
2200 [ # # ]: 0 : ABI_FREE(ifail)
2201 [ # # ]: 0 : ABI_FREE(rwork)
2202 [ # # ]: 0 : ABI_FREE(work_cplx)
2203 : 0 : RETURN
2204 : : end if ! cplex
2205 : :
2206 : : CASE (.TRUE.)
2207 : :
2208 : : #ifdef HAVE_LINALG_SCALAPACK
2209 : : ABI_ERROR("not coded yet")
2210 : : ! call Slk_processor%init(comm)
2211 : : ! istwf_k=1
2212 : :
2213 : : ! ! Initialize and fill Scalapack matrix from the global one.
2214 : : ! call Slk_matA%init(n,n,Slk_processor,istwf_k)
2215 : : ! call slk_matrix_from_global_dpc_2D(Slk_matA,uplo,a)
2216 : :
2217 : : ! call Slk_matB%init(n,n,Slk_processor,istwf_k)
2218 : : ! call slk_matrix_from_global_dpc_2D(Slk_matB,uplo,b)
2219 : :
2220 : : ! want_eigenvectors = firstchar(jobz,(/"V","v"/))
2221 : : ! if (want_eigenvectors) then ! Initialize the distributed vectors.
2222 : : ! call Slk_vec%init(n,n,Slk_processor,istwf_k)
2223 : : ! end if
2224 : :
2225 : : ! ! Solve the problem.
2226 : : ! ABI_ERROR("slk_pZHEGVX not coded yet")
2227 : : ! ! TODO write the scaLAPACK wrapper.
2228 : : ! call slk_pZHEGVX(itype,jobz,range,uplo,Slk_matA,Slk_matB,vl,vu,il,iu,abstol,Slk_vec,m,w)
2229 : :
2230 : : ! call Slk_matA%free()
2231 : : ! call Slk_matB%free()
2232 : : !
2233 : : ! if (want_eigenvectors) then ! A is overwritten with the eigenvectors
2234 : : ! z = czero
2235 : : ! call slk_matrix_to_global_dpc_2D(Slk_vec,"All",z) ! Fill the entries calculated by this node.
2236 : : ! call Slk_vec%free()
2237 : : ! call xmpi_sum(z,comm,ierr) ! Fill the remaining entries of the global matrix
2238 : : ! end if
2239 : :
2240 : : ! call Slk_processor%free()
2241 : :
2242 : : ! RETURN
2243 : : #endif
2244 : :
2245 : : ABI_BUG("You should not be here!")
2246 : : END SELECT
2247 : :
2248 : : end subroutine xhegvx_cplex
2249 : : !!***
2250 : :
2251 : : !----------------------------------------------------------------------
2252 : :
2253 : : !!****f* m_hide_lapack/wrap_CGEEV
2254 : : !! NAME
2255 : : !! wrap_CGEEV
2256 : : !!
2257 : : !! FUNCTION
2258 : : !! wrap_CGEEV computes for an N-by-N complex nonsymmetric matrix A, the
2259 : : !! eigenvalues and, optionally, the left and/or right eigenvectors using single precision arithmetic. [PRIVATE]
2260 : : !!
2261 : : !! The right eigenvector v(j) of A satisfies: A * v(j) = lambda(j) * v(j)
2262 : : !! where lambda(j) is its eigenvalue.
2263 : : !! The left eigenvector u(j) of A satisfies u(j)**H * A = lambda(j) * u(j)**H
2264 : : !! where u(j)**H denotes the conjugate transpose of u(j).
2265 : : !!
2266 : : !! The computed eigenvectors are normalized to have Euclidean norm
2267 : : !! equal to 1 and largest component real.
2268 : : !!
2269 : : !! INPUTS
2270 : : !! JOBVL (input) CHARACTER*1
2271 : : !! = 'N': left eigenvectors of A are not computed;
2272 : : !! = 'V': left eigenvectors of are computed.
2273 : : !!
2274 : : !! JOBVR (input) CHARACTER*1
2275 : : !! = 'N': right eigenvectors of A are not computed;
2276 : : !! = 'V': right eigenvectors of A are computed.
2277 : : !!
2278 : : !! N (input) INTEGER
2279 : : !! The order of the matrix A. N >= 0.
2280 : : !!
2281 : : !! LDA (input) INTEGER
2282 : : !! The leading dimension of the array A. LDA >= max(1,N).
2283 : : !!
2284 : : !! LDVL (input) INTEGER
2285 : : !! The leading dimension of the array VL. LDVL >= 1; if
2286 : : !! JOBVL = 'V', LDVL >= N.
2287 : : !!
2288 : : !! LDVR (input) INTEGER
2289 : : !! The leading dimension of the array VR. LDVR >= 1; if
2290 : : !! JOBVR = 'V', LDVR >= N.
2291 : : !!
2292 : : !! OUTPUT
2293 : : !! W (output) COMPLEX(sp) array, dimension (N)
2294 : : !! W contains the computed eigenvalues.
2295 : : !! VL (output) COMPLEX(SCP) array, dimension (LDVL,N)
2296 : : !! If JOBVL = 'V', the left eigenvectors u(j) are stored one
2297 : : !! after another in the columns of VL, in the same order
2298 : : !! as their eigenvalues.
2299 : : !! If JOBVL = 'N', VL is not referenced.
2300 : : !! u(j) = VL(:,j), the j-th column of VL.
2301 : : !! VR (output) COMPLEX(sp) array, dimension (LDVR,N)
2302 : : !! If JOBVR = 'V', the right eigenvectors v(j) are stored one
2303 : : !! after another in the columns of VR, in the same order
2304 : : !! as their eigenvalues.
2305 : : !! If JOBVR = 'N', VR is not referenced.
2306 : : !! v(j) = VR(:,j), the j-th column of VR.
2307 : : !!
2308 : : !! See also SIDE EFFECTS
2309 : : !!
2310 : : !! SIDE EFFECTS
2311 : : !! A (input/output) COMPLEX(sp) array, dimension (LDA,N)
2312 : : !! On entry, the N-by-N matrix A.
2313 : : !! On exit, A has been overwritten.
2314 : : !!
2315 : : !! SOURCE
2316 : :
2317 : 0 : subroutine wrap_CGEEV(jobvl, jobvr, n, a, lda, w, vl, ldvl, vr, ldvr)
2318 : :
2319 : : !Arguments ------------------------------------
2320 : : !scalars
2321 : : integer,intent(in) :: n,lda,ldvl,ldvr
2322 : : character(len=*),intent(in) :: jobvl,jobvr
2323 : : !arrays
2324 : : complex(sp),intent(inout) :: a(lda,n)
2325 : : complex(sp),intent(out) :: w(n)
2326 : : complex(sp),intent(out) :: vl(ldvl,n)
2327 : : complex(sp),intent(out) :: vr(ldvr,n)
2328 : :
2329 : : !Local variables ------------------------------
2330 : : !scalars
2331 : : integer :: info,lwork
2332 : : character(len=500) :: msg
2333 : : !arrays
2334 : 0 : real(sp),allocatable :: rwork(:)
2335 : 0 : complex(sp),allocatable :: work(:)
2336 : : !************************************************************************
2337 : :
2338 : 0 : lwork = MAX(1,2*n)
2339 : :
2340 [ # # # # : 0 : ABI_MALLOC(work,(lwork))
# # # # #
# ]
2341 [ # # # # : 0 : ABI_MALLOC(rwork,(2*n))
# # # # #
# # # #
# ]
2342 : :
2343 : 0 : call CGEEV(jobvl,jobvr,n,a,lda,w,vl,ldvl,vr,ldvr,work,lwork,rwork,info)
2344 : :
2345 [ # # ]: 0 : if (info < 0) then
2346 : 0 : write(msg,'(a,i0,a)')" The ",-info,"-th argument of CGEEV had an illegal value."
2347 : 0 : ABI_ERROR(msg)
2348 : : end if
2349 : :
2350 [ # # ]: 0 : if (info > 0) then
2351 : : write(msg,'(3a,i0,a,i0,a)')&
2352 : 0 : "CGEEV: The QR algorithm failed to compute all the eigenvalues, and no eigenvectors have been computed;",ch10,&
2353 : 0 : "Elements ",info+1,":",n," of W contain eigenvalues which have converged. "
2354 : 0 : ABI_ERROR(msg)
2355 : : end if
2356 : :
2357 [ # # ]: 0 : ABI_FREE(work)
2358 [ # # ]: 0 : ABI_FREE(rwork)
2359 : :
2360 : 0 : end subroutine wrap_CGEEV
2361 : : !!***
2362 : :
2363 : : !----------------------------------------------------------------------
2364 : :
2365 : : !!****f* m_hide_lapack/wrap_ZGEEV
2366 : : !! NAME
2367 : : !! wrap_ZGEEV
2368 : : !!
2369 : : !! FUNCTION
2370 : : !! wrap_ZGEEV computes for an N-by-N complex nonsymmetric matrix A, the
2371 : : !! eigenvalues and, optionally, the left and/or right eigenvectors using double precision arithmetic. [PRIVATE]
2372 : : !!
2373 : : !! The right eigenvector v(j) of A satisfies: A * v(j) = lambda(j) * v(j)
2374 : : !! where lambda(j) is its eigenvalue.
2375 : : !! The left eigenvector u(j) of A satisfies u(j)**H * A = lambda(j) * u(j)**H
2376 : : !! where u(j)**H denotes the conjugate transpose of u(j).
2377 : : !!
2378 : : !! The computed eigenvectors are normalized to have Euclidean norm
2379 : : !! equal to 1 and largest component real.
2380 : : !! No scalapack version is available (PZGEEV is not provided by the Scalapack team)
2381 : : !!
2382 : : !! INPUTS
2383 : : !! JOBVL (input) CHARACTER*1
2384 : : !! = 'N': left eigenvectors of A are not computed;
2385 : : !! = 'V': left eigenvectors of are computed.
2386 : : !!
2387 : : !! JOBVR (input) CHARACTER*1
2388 : : !! = 'N': right eigenvectors of A are not computed;
2389 : : !! = 'V': right eigenvectors of A are computed.
2390 : : !!
2391 : : !! N (input) INTEGER
2392 : : !! The order of the matrix A. N >= 0.
2393 : : !!
2394 : : !! LDA (input) INTEGER
2395 : : !! The leading dimension of the array A. LDA >= max(1,N).
2396 : : !!
2397 : : !! LDVL (input) INTEGER
2398 : : !! The leading dimension of the array VL. LDVL >= 1; if
2399 : : !! JOBVL = 'V', LDVL >= N.
2400 : : !!
2401 : : !! LDVR (input) INTEGER
2402 : : !! The leading dimension of the array VR. LDVR >= 1; if
2403 : : !! JOBVR = 'V', LDVR >= N.
2404 : : !!
2405 : : !! OUTPUT
2406 : : !! W (output) COMPLEX(DP) array, dimension (N)
2407 : : !! W contains the computed eigenvalues.
2408 : : !! VL (output) COMPLEX(DP) array, dimension (LDVL,N)
2409 : : !! If JOBVL = 'V', the left eigenvectors u(j) are stored one
2410 : : !! after another in the columns of VL, in the same order
2411 : : !! as their eigenvalues.
2412 : : !! If JOBVL = 'N', VL is not referenced.
2413 : : !! u(j) = VL(:,j), the j-th column of VL.
2414 : : !! VR (output) COMPLEX(DP) array, dimension (LDVR,N)
2415 : : !! If JOBVR = 'V', the right eigenvectors v(j) are stored one
2416 : : !! after another in the columns of VR, in the same order
2417 : : !! as their eigenvalues.
2418 : : !! If JOBVR = 'N', VR is not referenced.
2419 : : !! v(j) = VR(:,j), the j-th column of VR.
2420 : : !!
2421 : : !! See also SIDE EFFECTS
2422 : : !!
2423 : : !! SIDE EFFECTS
2424 : : !! A (input/output) COMPLEX(DP) array, dimension (LDA,N)
2425 : : !! On entry, the N-by-N matrix A.
2426 : : !! On exit, A has been overwritten.
2427 : : !!
2428 : : !! SOURCE
2429 : :
2430 : 10 : subroutine wrap_ZGEEV(jobvl,jobvr,n,a,lda,w,vl,ldvl,vr,ldvr)
2431 : :
2432 : : !Arguments ------------------------------------
2433 : : !scalars
2434 : : integer,intent(in) :: n,lda,ldvl,ldvr
2435 : : character(len=*),intent(in) :: jobvl,jobvr
2436 : : !arrays
2437 : : complex(dp),intent(inout) :: a(lda,n)
2438 : : complex(dp),intent(out) :: w(n)
2439 : : complex(dp),intent(out) :: vl(ldvl,n)
2440 : : complex(dp),intent(out) :: vr(ldvr,n)
2441 : :
2442 : : !Local variables ------------------------------
2443 : : !scalars
2444 : : integer :: info,lwork
2445 : : logical :: use_scalapack
2446 : : character(len=500) :: msg
2447 : : !arrays
2448 : 10 : real(dp),allocatable :: rwork(:)
2449 : 10 : complex(dp),allocatable :: work(:)
2450 : : !************************************************************************
2451 : :
2452 : 10 : use_scalapack=.FALSE.
2453 : :
2454 : : SELECT CASE(use_scalapack)
2455 : : CASE (.FALSE.)
2456 : :
2457 : 10 : lwork = MAX(1,2*n)
2458 [ - + + - : 30 : ABI_MALLOC(work,(lwork))
+ - - + -
+ ]
2459 [ + - - + : 40 : ABI_MALLOC(rwork,(2*n))
+ - + - +
- - + -
+ ]
2460 : :
2461 : 10 : call ZGEEV(jobvl,jobvr,n,a,lda,w,vl,ldvl,vr,ldvr,work,lwork,rwork,info)
2462 : :
2463 [ - + ]: 10 : if (info < 0) then
2464 : 0 : write(msg,'(a,i0,a)')" The ",-info,"-th argument of ZGEEV had an illegal value."
2465 : 0 : ABI_ERROR(msg)
2466 : : end if
2467 : :
2468 [ - + ]: 10 : if (info > 0) then
2469 : : write(msg,'(3a,i0,a,i0,a)')&
2470 : 0 : "ZGEEV: The QR algorithm failed to compute all the eigenvalues, and no eigenvectors have been computed;",ch10,&
2471 : 0 : "Elements ",info+1,":",n," of W contain eigenvalues which have converged. "
2472 : 0 : ABI_ERROR(msg)
2473 : : end if
2474 : :
2475 [ - + ]: 10 : ABI_FREE(work)
2476 [ - + ]: 10 : ABI_FREE(rwork)
2477 : : RETURN
2478 : :
2479 : : CASE (.TRUE.)
2480 : : ABI_BUG("You should not be here!")
2481 : : END SELECT
2482 : :
2483 : : end subroutine wrap_ZGEEV
2484 : : !!***
2485 : :
2486 : : !----------------------------------------------------------------------
2487 : :
2488 : : !!****f* m_hide_lapack/cginv
2489 : : !! NAME
2490 : : !! cginv
2491 : : !!
2492 : : !! FUNCTION
2493 : : !! Invert a general matrix of complex elements in single precision.
2494 : : !! CGETRF computes an LU factorization of a general N-by-N matrix A using partial pivoting with row interchanges.
2495 : : !! CGETRI computes the inverse of a matrix using the LU factorization computed by CGETRF.
2496 : : !!
2497 : : !! INPUTS
2498 : : !! n=size of complex matrix a
2499 : : !! a=matrix of complex elements
2500 : : !! [comm]=MPI communicator for ScaLAPACK inversion. Only available if the code has been compiled with Scalapack support.
2501 : : !! To avoid wasting CPU time the scalapack initialization is avoided if the number of processors in 1,
2502 : : !! in this case the sequential LAPACK routine is called.
2503 : : !!
2504 : : !! SIDE EFFECTS
2505 : : !! a(n,n)= array of complex elements, input, inverted at output
2506 : : !!
2507 : : !! TODO
2508 : : !! Add Scalapack version
2509 : : !!
2510 : : !! SOURCE
2511 : :
2512 : 0 : subroutine cginv(a, n, comm)
2513 : :
2514 : : !Arguments ------------------------------------
2515 : : !scalars
2516 : : integer,intent(in) :: n
2517 : : integer,optional,intent(in) :: comm
2518 : : !arrays
2519 : : complex(sp),intent(inout) :: a(n,n)
2520 : :
2521 : : !Local variables-------------------------------
2522 : : !scalars
2523 : : integer :: lwork,info,nprocs
2524 : : logical :: use_scalapack
2525 : : character(len=500) :: msg
2526 : : !arrays
2527 : 0 : integer,allocatable :: ipiv(:)
2528 : 0 : complex(sp),allocatable :: work(:)
2529 : : #ifdef HAVE_LINALG_SCALAPACK
2530 : : !integer :: ierr,istwf_k,ipiv_size,liwork
2531 : : !integer,allocatable :: iwork(:)
2532 : : !type(slkmat_dp_t) :: Slk_mat
2533 : : !type(slk_processor_t) :: Slk_processor
2534 : : #endif
2535 : : ! *************************************************************************
2536 : :
2537 : 0 : use_scalapack=.FALSE.
2538 [ # # ]: 0 : if (PRESENT(comm)) then
2539 : 0 : nprocs = xmpi_comm_size(comm)
2540 : : ! TODO
2541 : : !#ifdef HAVE_LINALG_SCALAPACK
2542 : : ! use_scalapack = (nprocs>1)
2543 : : !#endif
2544 : : end if
2545 : :
2546 : : SELECT CASE(use_scalapack)
2547 : :
2548 : : CASE (.FALSE.)
2549 [ # # # # : 0 : ABI_MALLOC(ipiv, (n))
# # # # #
# # # ]
2550 : :
2551 : 0 : call CGETRF(n,n,a,n,ipiv,info) ! P* L* U Factorization.
2552 : :
2553 [ # # ]: 0 : if (info < 0) then
2554 : 0 : write(msg,'(a,i0,a)')" The ",-info,"-th argument of CGETRF had an illegal value."
2555 : 0 : ABI_ERROR(msg)
2556 : : end if
2557 : :
2558 [ # # ]: 0 : if (info > 0) then
2559 : : write(msg,'(3a,i0,4a)')&
2560 : 0 : "The matrix that has been passed in argument is probably either singular or nearly singular.",ch10,&
2561 : 0 : "U(i,i) in the P*L*U factorization is exactly zero for i = ",info,ch10,&
2562 : 0 : "The factorization has been completed but the factor U is exactly singular.",ch10,&
2563 : 0 : "Division by zero will occur if it is used to solve a system of equations."
2564 : 0 : ABI_ERROR(msg)
2565 : : end if
2566 : :
2567 : 0 : lwork=MAX(1,n)
2568 [ # # # # : 0 : ABI_MALLOC(work,(lwork))
# # # # #
# ]
2569 : :
2570 : 0 : call CGETRI(n,a,n,ipiv,work,lwork,info) ! Inverts U and the computes inv(A)
2571 : :
2572 [ # # ]: 0 : if (info < 0) then
2573 : 0 : write(msg,'(a,i0,a)')" The ",-info,"-th argument of CGETRI had an illegal value."
2574 : 0 : ABI_ERROR(msg)
2575 : : end if
2576 : :
2577 [ # # ]: 0 : if (info > 0) then
2578 : : write(msg,'(3a,i0,a)')&
2579 : 0 : "The matrix that has been passed to this subroutine is probably either singular or nearly singular.",ch10,&
2580 : 0 : "U(i,i) for i= ",info," is exactly zero; the matrix is singular and its inverse could not be computed."
2581 : 0 : ABI_ERROR(msg)
2582 : : end if
2583 : :
2584 [ # # ]: 0 : ABI_FREE(ipiv)
2585 [ # # ]: 0 : ABI_FREE(work)
2586 : : RETURN
2587 : :
2588 : : CASE (.TRUE.)
2589 : :
2590 : : #if 0
2591 : : ! FIXME slkmat_dp_t does not have a single precision complex buffer
2592 : :
2593 : : #ifdef HAVE_LINALG_SCALAPACK
2594 : : call Slk_processor%init(comm)
2595 : : istwf_k=1
2596 : :
2597 : : ! Initialize and fill Scalapack matrix from the global one.
2598 : : call Slk_mat%init(n,n,Slk_processor,istwf_k)
2599 : :
2600 : : ! IMPORTANT NOTE: PZGETRF requires square block decomposition i.e., MB_A = NB_A.
2601 : : if ( Slk_mat%descript%tab(MB_)/=Slk_mat%descript%tab(NB_) ) then
2602 : : msg ="PZGETRF requires square block decomposition i.e., MB_A = NB_A."
2603 : : ABI_ERROR(msg)
2604 : : end if
2605 : :
2606 : : !!call slk_matrix_from_global_dpc_2D(Slk_mat,"All",a)
2607 : :
2608 : : ipiv_size = my_locr(Slk_mat) + Slk_mat%descript%tab(MB_)
2609 : : ABI_MALLOC(ipiv,(ipiv_size))
2610 : :
2611 : : call PCGETRF(Slk_mat%sizeb_global(1),Slk_mat%sizeb_global(2),Slk_mat%buffer_cplx_sp,&
2612 : : & 1,1,Slk_mat%descript%tab,ipiv,info) ! P * L * U Factorization.
2613 : :
2614 : : if (info/=0) then
2615 : : write(msg,'(a,i0)')"PCGETRF returned info= ",info
2616 : : ABI_ERROR(msg)
2617 : : end if
2618 : :
2619 : : ! Get optimal size of workspace for PCGETRI.
2620 : : lwork=-1; liwork=-1
2621 : : ABI_MALLOC(work,(1))
2622 : : ABI_MALLOC(iwork,(1))
2623 : :
2624 : : call PCGETRI(Slk_mat%sizeb_global(1),Slk_mat%buffer_cplx_sp,1,1,Slk_mat%descript%tab,ipiv,&
2625 : : & work,lwork,iwork,liwork,info)
2626 : :
2627 : : ABI_CHECK(info==0,"PZGETRI: Error during computation of workspace size")
2628 : :
2629 : : lwork = NINT(DBLE(work(1))); liwork=iwork(1)
2630 : : ABI_FREE(work)
2631 : : ABI_FREE(iwork)
2632 : :
2633 : : ! Solve the problem.
2634 : : ABI_MALLOC(work,(lwork))
2635 : : ABI_MALLOC(iwork,(liwork))
2636 : :
2637 : : call PCGETRI(Slk_mat%sizeb_global(1),Slk_mat%buffer_cplx_sp,1,1,Slk_mat%descript%tab,ipiv,&
2638 : : & work,lwork,iwork,liwork,info)
2639 : :
2640 : : if (info/=0) then
2641 : : write(msg,'(a,i0)')"PZGETRI returned info= ",info
2642 : : ABI_ERROR(msg)
2643 : : end if
2644 : :
2645 : : ABI_FREE(work)
2646 : : ABI_FREE(iwork)
2647 : : ABI_FREE(ipiv)
2648 : :
2649 : : ! Reconstruct the global matrix from the distributed one.
2650 : : a = czero
2651 : : !! call slk_matrix_to_global_dpc_2D(Slk_mat,"All",a) ! Fill the entries calculated by this node.
2652 : : call Slk_mat%free()
2653 : :
2654 : : call xmpi_sum(a,comm,ierr) ! Fill the remaining entries of the global matrix
2655 : : call Slk_processor%free()
2656 : :
2657 : : RETURN
2658 : : #endif
2659 : :
2660 : : #endif
2661 : :
2662 : : ABI_BUG("You should not be here!")
2663 : :
2664 : : END SELECT
2665 : :
2666 : : end subroutine cginv
2667 : : !!***
2668 : :
2669 : : !----------------------------------------------------------------------
2670 : :
2671 : : !!****f* m_hide_lapack/zginv
2672 : : !! NAME
2673 : : !! zginv
2674 : : !!
2675 : : !! FUNCTION
2676 : : !! Invert a general matrix of complex elements in double precision.
2677 : : !! ZGETRF computes an LU factorization of a general N-by-N matrix A using partial pivoting with row interchanges.
2678 : : !! ZGETRI computes the inverse of a matrix using the LU factorization computed by ZGETRF.
2679 : : !!
2680 : : !! INPUTS
2681 : : !! n=size of complex matrix a
2682 : : !! a=matrix of complex elements
2683 : : !! [comm]=MPI communicator for ScaLAPACK inversion. Only available if the code has been compiled with Scalapack support.
2684 : : !! To avoid wasting CPU time the scalapack initialization is avoided if the number of processors in 1.
2685 : : !! In this case the sequential LAPACK routine is called.
2686 : : !!
2687 : : !! SIDE EFFECTS
2688 : : !! a(n,n)= array of complex elements, input, inverted at output
2689 : : !!
2690 : : !! SOURCE
2691 : :
2692 : 8463501 : subroutine zginv(a, n, comm)
2693 : :
2694 : : !Arguments ------------------------------------
2695 : : !scalars
2696 : : integer,intent(in) :: n
2697 : : integer,optional,intent(in) :: comm
2698 : : !arrays
2699 : : complex(dp),intent(inout) :: a(n,n)
2700 : :
2701 : : !Local variables-------------------------------
2702 : : !scalars
2703 : : integer :: lwork,info,nprocs
2704 : : logical :: use_scalapack
2705 : : character(len=500) :: msg
2706 : : !arrays
2707 : 8463501 : integer,allocatable :: ipiv(:)
2708 : 8463501 : complex(dp),allocatable :: work(:)
2709 : : #ifdef HAVE_LINALG_SCALAPACK
2710 : : integer :: istwf_k,ierr
2711 : : type(slkmat_dp_t) :: Slk_mat
2712 : : type(slk_processor_t) :: Slk_processor
2713 : : #endif
2714 : : ! *************************************************************************
2715 : :
2716 : 8463501 : use_scalapack=.FALSE.
2717 [ + + ]: 8463501 : if (PRESENT(comm)) then
2718 : 5995 : nprocs = xmpi_comm_size(comm)
2719 : : #ifdef HAVE_LINALG_SCALAPACK
2720 : : use_scalapack = (nprocs>1)
2721 : : #endif
2722 : : end if
2723 : :
2724 : : SELECT CASE(use_scalapack)
2725 : : CASE (.FALSE.)
2726 [ - + + - : 25390503 : ABI_MALLOC(ipiv, (n))
+ - + - -
+ - + ]
2727 : 8463501 : call ZGETRF(n,n,a,n,ipiv,info) ! P* L* U Factorization.
2728 : :
2729 [ - + ]: 8463501 : if (info < 0) then
2730 : 0 : write(msg,'(a,i0,a)')" The ",-info,"-th argument of ZGETRF had an illegal value."
2731 : 0 : ABI_ERROR(msg)
2732 : : end if
2733 : :
2734 [ - + ]: 8463501 : if (info > 0) then
2735 : : write(msg,'(3a,i0,4a)')&
2736 : 0 : "The matrix that has been passed in argument is probably either singular or nearly singular.",ch10,&
2737 : 0 : "U(i,i) in the P*L*U factorization is exactly zero for i = ",info,ch10,&
2738 : 0 : "The factorization has been completed but the factor U is exactly singular.",ch10,&
2739 : 0 : "Division by zero will occur if it is used to solve a system of equations."
2740 : 0 : ABI_ERROR(msg)
2741 : : end if
2742 : :
2743 : 8463501 : lwork=MAX(1,n)
2744 [ - + + - : 25390503 : ABI_MALLOC(work,(lwork))
+ - - + -
+ ]
2745 : :
2746 : 8463501 : call ZGETRI(n,a,n,ipiv,work,lwork,info) ! Invert U and then compute inv(A)
2747 : :
2748 [ - + ]: 8463501 : if (info < 0) then
2749 : 0 : write(msg,'(a,i0,a)')" The ",-info,"-th argument of ZGETRI had an illegal value."
2750 : 0 : ABI_ERROR(msg)
2751 : : end if
2752 : :
2753 [ - + ]: 8463501 : if (info > 0) then
2754 : : write(msg,'(3a,i0,a)')&
2755 : 0 : "The matrix that has been passed to this subroutine is probably either singular or nearly singular.",ch10,&
2756 : 0 : "U(i,i) for i= ",info," is exactly zero; the matrix is singular and its inverse could not be computed."
2757 : 0 : ABI_ERROR(msg)
2758 : : end if
2759 : :
2760 [ - + ]: 8463501 : ABI_FREE(ipiv)
2761 [ - + ]: 8463501 : ABI_FREE(work)
2762 : : RETURN
2763 : :
2764 : : CASE (.TRUE.)
2765 : :
2766 : : #ifdef HAVE_LINALG_SCALAPACK
2767 : : call Slk_processor%init(comm)
2768 : : istwf_k=1
2769 : :
2770 : : ! Initialize and fill Scalapack matrix from the global one.
2771 : : call Slk_mat%init(n,n,Slk_processor,istwf_k)
2772 : : call slk_matrix_from_global_dpc_2D(Slk_mat,"All",a)
2773 : :
2774 : : ! Perform the calculation with scaLAPACK.
2775 : : call Slk_mat%invert()
2776 : :
2777 : : ! Reconstruct the global matrix from the distributed one.
2778 : : a = czero
2779 : : call slk_matrix_to_global_dpc_2D(Slk_mat,"All",a) ! Fill the entries calculated by this node.
2780 : : call Slk_mat%free()
2781 : :
2782 : : call xmpi_sum(a,comm,ierr) ! Fill the remaining entries of the global matrix
2783 : : call Slk_processor%free()
2784 : :
2785 : : return
2786 : : #endif
2787 : :
2788 : : ABI_BUG("You should not be here!")
2789 : : END SELECT
2790 : :
2791 : : end subroutine zginv
2792 : : !!***
2793 : :
2794 : : !----------------------------------------------------------------------
2795 : :
2796 : : !!****f* m_hide_lapack/zhpd_invert
2797 : : !! NAME
2798 : : !! zhpd_invert
2799 : : !!
2800 : : !! FUNCTION
2801 : : !! Invert a Hermitian positive definite matrix of complex elements in double precision.
2802 : : !!
2803 : : !! INPUTS
2804 : : !! uplo= 'U': Upper triangle of A is stored;
2805 : : !! = 'L': Lower triangle of A is stored.
2806 : : !! n=size of complex matrix a
2807 : : !! a=matrix of complex elements
2808 : : !! [comm]=MPI communicator for ScaLAPACK inversion. Only available if the code has been compiled with Scalapack support.
2809 : : !! To avoid wasting CPU time the scalapack initialization is avoided if the number of processors in 1.
2810 : : !! In this case the sequential LAPACK routine is called.
2811 : : !!
2812 : : !! SIDE EFFECTS
2813 : : !! a(n,n)=
2814 : : !! On entry, the Hermitian matrix A. If UPLO = 'U', the leading
2815 : : !! N-by-N upper triangular part of A contains the upper
2816 : : !! triangular part of the matrix A, and the strictly lower
2817 : : !! triangular part of A is not referenced. If UPLO = 'L', the
2818 : : !! leading N-by-N lower triangular part of A contains the lower
2819 : : !! triangular part of the matrix A, and the strictly upper
2820 : : !! triangular part of A is not referenced.
2821 : : !! On exit, the upper or lower triangle of the (Hermitian) inverse of A
2822 : : !!
2823 : : !! SOURCE
2824 : :
2825 : 1 : subroutine zhpd_invert(uplo, a, n, comm)
2826 : :
2827 : : !Arguments ------------------------------------
2828 : : !scalars
2829 : : character(len=*),intent(in) :: uplo
2830 : : integer,intent(in) :: n
2831 : : integer,optional,intent(in) :: comm
2832 : : !arrays
2833 : : complex(dp),intent(inout) :: a(n,n)
2834 : :
2835 : : !Local variables-------------------------------
2836 : : !scalars
2837 : : integer :: info,nprocs
2838 : : logical :: use_scalapack
2839 : : character(len=500) :: msg
2840 : : !arrays
2841 : : #ifdef HAVE_LINALG_SCALAPACK
2842 : : integer :: istwf_k,ierr
2843 : : type(slkmat_dp_t) :: Slk_mat
2844 : : type(slk_processor_t) :: Slk_processor
2845 : : #endif
2846 : : ! *************************************************************************
2847 : :
2848 : 1 : use_scalapack=.FALSE.
2849 [ - + ]: 1 : if (PRESENT(comm)) then
2850 : 0 : nprocs = xmpi_comm_size(comm)
2851 : : #ifdef HAVE_LINALG_SCALAPACK
2852 : : use_scalapack = (nprocs>1)
2853 : : #endif
2854 : : end if
2855 : :
2856 : : SELECT CASE(use_scalapack)
2857 : : CASE (.FALSE.)
2858 : : ! * ZPOTRF computes the Cholesky factorization of a complex Hermitian positive definite.
2859 : : ! * A = U**H * U, if UPLO = 'U', or
2860 : : ! * A = L * L**H, if UPLO = 'L',
2861 : 1 : call ZPOTRF(uplo,n,a,n,info)
2862 : :
2863 [ - + ]: 1 : if (info < 0) then
2864 : 0 : write(msg,'(a,i0,a)')" The ",-info,"-th argument of ZPOTRF had an illegal value."
2865 : 0 : ABI_ERROR(msg)
2866 : : end if
2867 : :
2868 [ - + ]: 1 : if (info > 0) then
2869 : : write(msg,'(a,i0,3a)')&
2870 : 0 : "The leading minor of order ",info," is not positive definite, ",ch10,&
2871 : 0 : "and the factorization could not be completed."
2872 : 0 : ABI_ERROR(msg)
2873 : : end if
2874 : : !
2875 : : ! * ZPOTRI computes the inverse of a complex Hermitian positive definite
2876 : : ! * matrix A using the Cholesky factorization A = U**H*U or A = L*L**H
2877 : : ! * computed by ZPOTRF.
2878 : : ! * On exit, the upper or lower triangle of the (Hermitian)
2879 : : ! * inverse of A, overwriting the input factor U or L.
2880 : 1 : call ZPOTRI(uplo,n,a,n,info)
2881 : :
2882 [ - + ]: 1 : if (info < 0) then
2883 : 0 : write(msg,'(a,i0,a)')" The ",-info,"-th argument of ZPOTRI had an illegal value."
2884 : 0 : ABI_ERROR(msg)
2885 : : end if
2886 : :
2887 [ - + ]: 1 : if (info > 0) then
2888 : : write(msg,'(a,2(1x,i0),a)')&
2889 : 0 : "The ( ",info,info,")element of the factor U or L is zero, and the inverse could not be computed."
2890 : 0 : ABI_ERROR(msg)
2891 : : end if
2892 : :
2893 : 1 : RETURN
2894 : :
2895 : : CASE (.TRUE.)
2896 : :
2897 : : #ifdef HAVE_LINALG_SCALAPACK
2898 : : call Slk_processor%init(comm)
2899 : : istwf_k=1
2900 : :
2901 : : ! Initialize and fill Scalapack matrix from the global one.
2902 : : call Slk_mat%init(n,n,Slk_processor,istwf_k)
2903 : : call slk_matrix_from_global_dpc_2D(Slk_mat,uplo,a)
2904 : :
2905 : : ! Perform the calculation with scaLAPACK.
2906 : : call Slk_mat%hpd_invert(uplo, full=.False.)
2907 : :
2908 : : ! Reconstruct the global matrix from the distributed one.
2909 : : a = czero
2910 : : call slk_matrix_to_global_dpc_2D(Slk_mat,uplo,a) ! Fill the entries calculated by this node.
2911 : : call Slk_mat%free()
2912 : :
2913 : : call xmpi_sum(a,comm,ierr) ! Fill the remaining entries of the global matrix
2914 : : call Slk_processor%free()
2915 : :
2916 : : RETURN
2917 : : #endif
2918 : :
2919 : : ABI_BUG("You should not be here!")
2920 : : END SELECT
2921 : :
2922 : : end subroutine zhpd_invert
2923 : : !!***
2924 : :
2925 : : !----------------------------------------------------------------------
2926 : :
2927 : : !!****f* m_hide_lapack/matrginv
2928 : : !! NAME
2929 : : !! matrginv
2930 : : !!
2931 : : !! FUNCTION
2932 : : !! Invert a general matrix of real*8 elements.
2933 : : !!
2934 : : !! INPUTS
2935 : : !! lda=leading dimension of complex matrix a
2936 : : !! n=size of complex matrix a
2937 : : !! a=matrix of real elements
2938 : : !! OUTPUT
2939 : : !! a=inverse of a input matrix
2940 : : !!
2941 : : !! SIDE EFFECTS
2942 : : !! a(lda,n)= array of real elements, input, inverted at output
2943 : : !!
2944 : : !! SOURCE
2945 : :
2946 : 9242 : subroutine matrginv(a,lda,n)
2947 : :
2948 : : !Arguments ------------------------------------
2949 : : !scalars
2950 : : integer,intent(in) :: lda,n
2951 : : !arrays
2952 : : real(dp),intent(inout) :: a(lda,n)
2953 : :
2954 : : !Local variables-------------------------------
2955 : : !scalars
2956 : : integer :: ierr,nwork
2957 : : #if defined HAVE_LINALG_ESSL
2958 : : real(dp) :: rcond
2959 : : #endif
2960 : : character(len=500) :: message
2961 : : !arrays
2962 : 9242 : integer,allocatable :: ipvt(:)
2963 : : #if defined HAVE_LINALG_ESSL
2964 : : real(dp) :: det(2)
2965 : : #elif defined HAVE_LINALG_ASL
2966 : : real(dp) :: det(2)
2967 : : #endif
2968 : 9242 : real(dp),allocatable :: work(:)
2969 : : ! *************************************************************************
2970 : :
2971 : : #if defined HAVE_LINALG_ESSL
2972 : : nwork=200*n
2973 : : #else
2974 : 9242 : nwork=n
2975 : : #endif
2976 : :
2977 [ - + + - : 27726 : ABI_MALLOC(work,(nwork))
+ - + - -
+ - + ]
2978 [ - + + - : 27726 : ABI_MALLOC(ipvt,(n))
+ - + - -
+ - + ]
2979 : :
2980 : : #if defined HAVE_LINALG_ESSL
2981 : :
2982 : : call dgeicd(a,lda,n,0,rcond,det,work,nwork)
2983 : : if(abs(rcond)==zero) then
2984 : : write(message, '(7a)' )&
2985 : : ' The matrix that has been passed in argument of this subroutine',ch10,&
2986 : : ' is probably either singular or nearly singular.',ch10,&
2987 : : ' The ESSL routine dgeicd failed.',ch10,&
2988 : : ' Action: Contact ABINIT group '
2989 : : ABI_ERROR(message)
2990 : : end if
2991 : :
2992 : : #elif defined HAVE_LINALG_ASL
2993 : :
2994 : : call dbgmlu(a,lda,n,ipvt,ierr)
2995 : : if(ierr /= 0) then
2996 : : write(message, '(7a)' ) ch10,&
2997 : : ' The matrix that has been passed in argument of this subroutine',ch10,&
2998 : : ' is probably either singular or nearly singular.',ch10,&
2999 : : ' The ASL routine dbgmlu failed.',ch10,&
3000 : : ' Action: Contact ABINIT group '
3001 : : ABI_ERROR(message)
3002 : : end if
3003 : :
3004 : : call dbgmdi(a,lda,n,ipvt,det,-1,work,ierr)
3005 : :
3006 : : if(ierr /= 0) then
3007 : : write(message, '(7a)' ) &
3008 : : ' The matrix that has been passed in argument of this subroutine',ch10,&
3009 : : ' is probably either singular or nearly singular.',ch10,&
3010 : : ' The ASL routine dbgmdi failed.',ch10,&
3011 : : ' Action: Contact ABINIT group '
3012 : : ABI_ERROR(message)
3013 : : end if
3014 : :
3015 : : #else
3016 : :
3017 : 9242 : call dgetrf(n,n,a,lda,ipvt,ierr)
3018 [ - + ]: 9242 : if(ierr /= 0) then
3019 : : write(message, '(7a)' ) &
3020 : 0 : ' The matrix that has been passed in argument of this subroutine',ch10,&
3021 : 0 : ' is probably either singular or nearly singular.',ch10,&
3022 : 0 : ' The LAPACK routine dgetrf failed.',ch10,&
3023 : 0 : ' Action: Contact ABINIT group '
3024 : 0 : ABI_ERROR(message)
3025 : : end if
3026 : :
3027 : 9242 : call dgetri(n,a,lda,ipvt,work,n,ierr)
3028 : :
3029 [ - + ]: 9242 : if(ierr /= 0) then
3030 : : write(message, '(7a)' ) &
3031 : 0 : ' The matrix that has been passed in argument of this subroutine',ch10,&
3032 : 0 : ' is probably either singular or nearly singular.',ch10,&
3033 : 0 : ' The LAPACK routine dgetri failed.',ch10,&
3034 : 0 : ' Action: Contact ABINIT group '
3035 : 0 : ABI_ERROR(message)
3036 : : end if
3037 : :
3038 : : #endif
3039 : :
3040 [ - + ]: 9242 : ABI_FREE(work)
3041 [ - + ]: 9242 : ABI_FREE(ipvt)
3042 : :
3043 : 9242 : end subroutine matrginv
3044 : : !!***
3045 : :
3046 : : !!****f* m_hide_lapack/matr3eigval
3047 : : !! NAME
3048 : : !! matr3eigval
3049 : : !!
3050 : : !! FUNCTION
3051 : : !! Find the eigenvalues of a real symmetric 3x3 matrix, entered in full storage mode.
3052 : : !!
3053 : : !! INPUTS
3054 : : !! matr(3,3)=real symmetric 3x3 matrix
3055 : : !!
3056 : : !! OUTPUT
3057 : : !! eigval(3)=three eigenvalues
3058 : : !!
3059 : : !! SOURCE
3060 : :
3061 : 11540 : subroutine matr3eigval(eigval,matr)
3062 : :
3063 : : !Arguments ------------------------------------
3064 : : !arrays
3065 : : real(dp),intent(in) :: matr(3,3)
3066 : : real(dp),intent(out) :: eigval(3)
3067 : :
3068 : : !Local variables-------------------------------
3069 : : !scalars
3070 : : integer :: ier
3071 : : !arrays
3072 : : real(dp) :: eigvec(2,3,3),matrx(2,6),zhpev1(2,2*3-1),zhpev2(3*3-2)
3073 : : ! *************************************************************************
3074 : :
3075 : 11540 : matrx(1,1)=matr(1,1)
3076 : 11540 : matrx(1,2)=matr(1,2)
3077 : 11540 : matrx(1,3)=matr(2,2)
3078 : 11540 : matrx(1,4)=matr(1,3)
3079 : 11540 : matrx(1,5)=matr(2,3)
3080 : 11540 : matrx(1,6)=matr(3,3)
3081 [ + + ]: 80780 : matrx(2,:)=zero
3082 : :
3083 : 11540 : call ZHPEV ('V','U',3,matrx,eigval,eigvec,3,zhpev1,zhpev2,ier)
3084 : : !write(std_out,*)' eigval=',eigval
3085 : :
3086 : 11540 : end subroutine matr3eigval
3087 : : !!***
3088 : :
3089 : : !!****f* ABINIT/jacobi
3090 : : !! NAME
3091 : : !! jacobi
3092 : : !!
3093 : : !! FUNCTION
3094 : : !! Computes all eigenvalues and eigenvectors of a real symmetric matrix a,
3095 : : !! which is of size n by n, stored in a physical np by np array. On output,
3096 : : !! elements of a above the diagonal are destroyed. d returns the
3097 : : !! eigenvalues of a in its first n elements. v is a matrix with the same
3098 : : !! logical and physical dimensions as a, whose columns contain, on output,
3099 : : !! the normalized eigenvectors of a. nrot returns the number of Jacobi
3100 : : !! rotations that were required.
3101 : : !!
3102 : : !! INPUTS
3103 : : !!
3104 : : !! OUTPUT
3105 : : !!
3106 : : !! NOTES
3107 : : !! This routine is deprecated, use Lapack API
3108 : : !!
3109 : : !! SOURCE
3110 : :
3111 : 119591 : subroutine jacobi(a,n,np,d,v,nrot)
3112 : :
3113 : : !Arguments
3114 : : integer :: n,np,nrot
3115 : : real*8 :: a(np,np),d(np),v(np,np)
3116 : : !Local variables
3117 : : integer, parameter :: NMAX=500
3118 : : integer i,ip,iq,j
3119 : : real*8 c,g,h,s,sm,t,tau,theta,tresh,b(NMAX),z(NMAX)
3120 [ + + ]: 484092 : do ip=1,n
3121 [ + + ]: 1480916 : do iq=1,n
3122 : 1480916 : v(ip,iq)=0.
3123 : : enddo
3124 : 484092 : v(ip,ip)=1.
3125 : : enddo
3126 [ + + ]: 484092 : do ip=1,n
3127 : 364501 : b(ip)=a(ip,ip)
3128 : 364501 : d(ip)=b(ip)
3129 : 484092 : z(ip)=0.
3130 : : enddo
3131 : 119591 : nrot=0
3132 [ + - ]: 724069 : do i=1,50
3133 : : sm=0.
3134 [ + + ]: 2209846 : do ip=1,n-1
3135 [ + + ]: 4494970 : do iq=ip+1,n
3136 : 3770901 : sm=sm+abs(a(ip,iq))
3137 : : enddo
3138 : : enddo
3139 [ + + ]: 724069 : if(sm.eq.0.)return
3140 [ + + ]: 604478 : if(i.lt.4)then
3141 : 358692 : tresh=0.2*sm/n**2
3142 : : else
3143 : : tresh=0.
3144 : : endif
3145 [ + + ]: 1845345 : do ip=1,n-1
3146 [ + + ]: 3754512 : do iq=ip+1,n
3147 : 1909167 : g=100.*abs(a(ip,iq))
3148 : : if((i.gt.4).and.(abs(d(ip))+g.eq.abs(d(ip))) &
3149 [ + + + + : 3150034 : & .and.(abs(d(iq))+g.eq.abs(d(iq))))then
+ + ]
3150 : 365245 : a(ip,iq)=0.
3151 [ + + ]: 1543922 : else if(abs(a(ip,iq)).gt.tresh)then
3152 : 1114296 : h=d(iq)-d(ip)
3153 [ + + ]: 1114296 : if(abs(h)+g.eq.abs(h))then
3154 : 122227 : t=a(ip,iq)/h
3155 : : else
3156 : 992069 : theta=0.5*h/a(ip,iq)
3157 : 992069 : t=1./(abs(theta)+sqrt(1.+theta**2))
3158 [ + + ]: 992069 : if(theta.lt.0.)t=-t
3159 : : endif
3160 : 1114296 : c=1./sqrt(1+t**2)
3161 : 1114296 : s=t*c
3162 : 1114296 : tau=s/(1.+c)
3163 : 1114296 : h=t*a(ip,iq)
3164 : 1114296 : z(ip)=z(ip)-h
3165 : 1114296 : z(iq)=z(iq)+h
3166 : 1114296 : d(ip)=d(ip)-h
3167 : 1114296 : d(iq)=d(iq)+h
3168 : 1114296 : a(ip,iq)=0.
3169 [ + + ]: 1497377 : do j=1,ip-1
3170 : 383081 : g=a(j,ip)
3171 : 383081 : h=a(j,iq)
3172 : 383081 : a(j,ip)=g-s*(h+g*tau)
3173 : 1497377 : a(j,iq)=h+s*(g-h*tau)
3174 : : enddo
3175 [ + + ]: 1526917 : do j=ip+1,iq-1
3176 : 412621 : g=a(ip,j)
3177 : 412621 : h=a(j,iq)
3178 : 412621 : a(ip,j)=g-s*(h+g*tau)
3179 : 1526917 : a(j,iq)=h+s*(g-h*tau)
3180 : : enddo
3181 [ + + ]: 1537000 : do j=iq+1,n
3182 : 422704 : g=a(ip,j)
3183 : 422704 : h=a(iq,j)
3184 : 422704 : a(ip,j)=g-s*(h+g*tau)
3185 : 1537000 : a(iq,j)=h+s*(g-h*tau)
3186 : : enddo
3187 [ + + ]: 4561294 : do j=1,n
3188 : 3446998 : g=v(j,ip)
3189 : 3446998 : h=v(j,iq)
3190 : 3446998 : v(j,ip)=g-s*(h+g*tau)
3191 : 4561294 : v(j,iq)=h+s*(g-h*tau)
3192 : : enddo
3193 : 1114296 : nrot=nrot+1
3194 : : endif
3195 : : enddo
3196 : : enddo
3197 [ + + ]: 2449823 : do ip=1,n
3198 : 1845345 : b(ip)=b(ip)+z(ip)
3199 : 1845345 : d(ip)=b(ip)
3200 : 2449823 : z(ip)=0.
3201 : : enddo
3202 : : enddo
3203 : 0 : write(std_out,*) 'too many iterations in jacobi'
3204 : :
3205 : : end subroutine jacobi
3206 : : !!***
3207 : :
3208 : : !!****f* m_hide_lapack/ludcmp
3209 : : !! NAME
3210 : : !! ludcmp
3211 : : !!
3212 : : !! FUNCTION
3213 : : !! Given a matrix a(1:n,1:n), with physical dimension np by np, this
3214 : : !! routine replaces it by the LU decomposition of a rowwise permutation of
3215 : : !! itself. a and n are input. a is output, arranged as in equation (2.3.14)
3216 : : !! above; indx(1:n) is an output vector that records the row permutation
3217 : : !! effected by the partial pivoting; id is output as +- 1 depending on
3218 : : !! whether the number of row interchanges was even or odd,
3219 : : !! respectively. This routine is used in combination with lubksb to solve
3220 : : !! linear equations or invert a matrix.
3221 : : !!
3222 : : !! NOTES
3223 : : !! This routine is deprecated, use lapack API
3224 : : !!
3225 : : !! SOURCE
3226 : :
3227 : 2071 : SUBROUTINE ludcmp(a,n,np,indx,id,info)
3228 : :
3229 : : INTEGER n,np,indx(n),NMAX,id,info
3230 : : REAL*8 a(np,np),TINY
3231 : : PARAMETER (NMAX=500,TINY=1.0e-20)
3232 : :
3233 : : INTEGER i,imax,j,k
3234 : : REAL*8 aamax,dum,sum,vv(NMAX)
3235 : :
3236 : : ! write(std_out,*) 'ENTERING LUDCMP...'
3237 : : ! write(std_out,*) 'in ludcmp n=',n,' np=',np
3238 : : ! write(std_out,201) ((a(i,j),j=1,n),i=1,n)
3239 : : ! 201 FORMAT('A in ludcmp ',/,3F16.8,/,3F16.8,/,3F16.8)
3240 : 2071 : id=1
3241 : 2071 : info=0
3242 [ + + ]: 8284 : do i=1,n
3243 : : aamax=0.
3244 [ + + ]: 24852 : do j=1,n
3245 : 24852 : if (abs(a(i,j)).gt.aamax) aamax=abs(a(i,j))
3246 : : enddo
3247 [ - + ]: 6213 : if (aamax.eq.0.) then
3248 : 0 : write(std_out,*) 'LUDCMP: singular matrix !!!'
3249 [ # # ]: 0 : do j=1,3
3250 : 0 : write(std_out,*) (a(j,k),k=1,3)
3251 : : enddo
3252 : 0 : info=1
3253 : 0 : return
3254 : : ! stop 'singular matrix in ludcmp'
3255 : : endif
3256 : 8284 : vv(i)=1./aamax
3257 : : enddo
3258 [ + + ]: 8284 : do j=1,n
3259 [ + + ]: 12426 : do i=1,j-1
3260 : 6213 : sum=a(i,j)
3261 [ + + ]: 8284 : do k=1,i-1
3262 : 8284 : sum=sum-a(i,k)*a(k,j)
3263 : : enddo
3264 : 12426 : a(i,j)=sum
3265 : : enddo
3266 : : aamax=0.
3267 [ + + ]: 18639 : do i=j,n
3268 : 12426 : sum=a(i,j)
3269 [ + + ]: 20710 : do k=1,j-1
3270 : 20710 : sum=sum-a(i,k)*a(k,j)
3271 : : enddo
3272 : 12426 : a(i,j)=sum
3273 : 12426 : dum=vv(i)*abs(sum)
3274 [ + + ]: 18639 : if (dum.ge.aamax) then
3275 : 6940 : imax=i
3276 : 6940 : aamax=dum
3277 : : endif
3278 : : enddo
3279 [ + + ]: 6213 : if (j.ne.imax)then
3280 [ + + ]: 2768 : do k=1,n
3281 : 2076 : dum=a(imax,k)
3282 : 2076 : a(imax,k)=a(j,k)
3283 : 2768 : a(j,k)=dum
3284 : : enddo
3285 : 692 : id=-id
3286 : 692 : vv(imax)=vv(j)
3287 : : endif
3288 : 6213 : indx(j)=imax
3289 [ - + ]: 6213 : if(a(j,j).eq.0.)a(j,j)=TINY
3290 [ + + ]: 8284 : if(j.ne.n)then
3291 : 4142 : dum=1./a(j,j)
3292 [ + + ]: 10355 : do i=j+1,n
3293 : 10355 : a(i,j)=a(i,j)*dum
3294 : : enddo
3295 : : endif
3296 : : enddo
3297 : : ! write(std_out,*) 'LEAVING LUDCMP...'
3298 : : return
3299 : : END SUBROUTINE ludcmp
3300 : : !!***
3301 : :
3302 : : !!****f* m_hide_lapack/lubksb
3303 : : !! NAME
3304 : : !! lubksb
3305 : : !!
3306 : : !! FUNCTION
3307 : : !! Solves the set of n linear equations A . X = B. Here a is input, not as
3308 : : !! the matrix A but rather as its LU decomposition, determined by the
3309 : : !! routine ludcmp. indx is input as the permutation vector returned by
3310 : : !! ludcmp. b(1:n) is input as the right-hand side vector B, and returns
3311 : : !! with the solution vector X. a, n, np, and indx are not modified by this
3312 : : !! routine and can be left in place for successive calls with different
3313 : : !! right-hand sides b. This routine takes into account the possibility that
3314 : : !! b will begin with many zero elements, so it is efficient for use in
3315 : : !! matrix inversion.
3316 : : !!
3317 : : !! NOTES
3318 : : !! This routine is deprecated, use lapack API
3319 : : !!
3320 : : !! SOURCE
3321 : :
3322 : 6213 : SUBROUTINE lubksb(a,n,np,indx,b)
3323 : :
3324 : : INTEGER n,np,indx(n)
3325 : : REAL*8 a(np,np),b(n)
3326 : :
3327 : : INTEGER i,ii,j,ll
3328 : : REAL*8 sum
3329 : : ! write(std_out,*) 'ENTERING LUBKSB...'
3330 : : ! write(std_out,201) ((a(i,j),j=1,n),i=1,n)
3331 : : ! 201 FORMAT('A in lubksb ',/,3F16.8,/,3F16.8,/,3F16.8)
3332 : :
3333 : 6213 : ii=0
3334 [ + + ]: 24852 : do i=1,n
3335 : 18639 : ll=indx(i)
3336 : 18639 : sum=b(ll)
3337 : 18639 : b(ll)=b(i)
3338 [ + + ]: 18639 : if (ii.ne.0)then
3339 [ + + ]: 14497 : do j=ii,i-1
3340 : 14497 : sum=sum-a(i,j)*b(j)
3341 : : enddo
3342 [ + + ]: 12426 : else if (sum.ne.0.) then
3343 : 6213 : ii=i
3344 : : endif
3345 : 24852 : b(i)=sum
3346 : : enddo
3347 [ + + ]: 24852 : do i=n,1,-1
3348 : 18639 : sum=b(i)
3349 [ + + ]: 37278 : do j=i+1,n
3350 : 37278 : sum=sum-a(i,j)*b(j)
3351 : : enddo
3352 : 24852 : b(i)=sum/a(i,i)
3353 : : enddo
3354 : : ! write(std_out,*) 'LEAVING LUBKSB...'
3355 : 6213 : return
3356 : :
3357 : : END SUBROUTINE LUBKSB
3358 : : !!***
3359 : :
3360 : : !!****f* m_hide_lapack/dzegdi
3361 : : !! NAME
3362 : : !! dzgedi
3363 : : !!
3364 : : !! FUNCTION
3365 : : !! This routine is the clone of zgefa.F90 using real*8 a(2) instead of complex*16
3366 : : !! for the purpose of ABINIT
3367 : : !!
3368 : : !! NOTES
3369 : : !! This routine is deprecated, use lapack API
3370 : : !!
3371 : : !! SOURCE
3372 : :
3373 : 2394114 : subroutine dzgedi(a,lda,n,ipvt,det,work,job)
3374 : :
3375 : : integer :: lda,n,ipvt(n),job
3376 : : real*8 :: a(2,lda,n),det(2,2),work(2,n)
3377 : : !
3378 : : ! zgedi computes the determinant and inverse of a matrix
3379 : : ! using the factors computed by zgeco or zgefa.
3380 : : !
3381 : : ! on entry
3382 : : !
3383 : : ! a complex*16(lda, n)
3384 : : ! the output from zgeco or zgefa.
3385 : : !
3386 : : ! lda integer
3387 : : ! the leading dimension of the array a .
3388 : : !
3389 : : ! n integer
3390 : : ! the order of the matrix a .
3391 : : !
3392 : : ! ipvt integer(n)
3393 : : ! the pivot vector from zgeco or zgefa.
3394 : : !
3395 : : ! work complex*16(n)
3396 : : ! work vector. contents destroyed.
3397 : : !
3398 : : ! job integer
3399 : : ! = 11 both determinant and inverse.
3400 : : ! = 01 inverse only.
3401 : : ! = 10 determinant only.
3402 : : !
3403 : : ! on return
3404 : : !
3405 : : ! a inverse of original matrix if requested.
3406 : : ! otherwise unchanged.
3407 : : !
3408 : : ! det complex*16(2)
3409 : : ! determinant of original matrix if requested.
3410 : : ! otherwise not referenced.
3411 : : ! determinant = det(1) * 10.0**det(2)
3412 : : ! with 1.0 .le. cabs1(det(1)) .lt. 10.0
3413 : : ! or det(1) .eq. 0.0 .
3414 : : !
3415 : : ! error condition
3416 : : !
3417 : : ! a division by zero will occur if the input factor contains
3418 : : ! a zero on the diagonal and the inverse is requested.
3419 : : ! it will not occur if the subroutines are called correctly
3420 : : ! and if zgeco has set rcond .gt. 0.0 or zgefa has set
3421 : : ! info .eq. 0 .
3422 : : !
3423 : : ! linpack. this version dated 08/14/78 .
3424 : : ! cleve moler, university of new mexico, argonne national lab.
3425 : : !
3426 : : ! subroutines and functions
3427 : : !
3428 : : ! internal variables
3429 : : !
3430 : : double precision :: r(2),rk(2),rkj(2)
3431 : : double precision :: ten,rinv2,rabs
3432 : : integer :: i,j,k,kb,kp1,l,nm1
3433 : : !
3434 : : ! compute determinant
3435 : : !
3436 [ + + ]: 2394114 : if (job/10 .eq. 0) go to 70
3437 : 2363834 : det(1,1) = 1.0d0; det(2,1) = 0.0d0
3438 : 2363834 : det(1,2) = 0.0d0; det(2,2) = 0.0d0
3439 : 2363834 : ten = 10.0d0
3440 [ + + ]: 31739050 : do i = 1, n
3441 [ + + ]: 29375216 : if (ipvt(i) .ne. i) then
3442 : 12598326 : det(1,1) = -det(1,1)
3443 : 12598326 : det(2,1) = -det(2,1)
3444 : : end if
3445 : 29375216 : r(1)=det(1,1); r(2)=det(2,1)
3446 : 29375216 : det(1,1) = r(1)*a(1,i,i)-r(2)*a(2,i,i)
3447 : 29375216 : det(2,1) = r(2)*a(1,i,i)+r(1)*a(2,i,i)
3448 : : ! ...exit
3449 : 29375216 : rabs = abs(det(1,1))+abs(det(2,1))
3450 [ - + ]: 29375216 : if (rabs .eq. 0.0d0) go to 60
3451 : : 10 continue
3452 : 33686378 : rabs = abs(det(1,1))+abs(det(2,1))
3453 [ + + ]: 33686378 : if (rabs .ge. 1.0d0) go to 20
3454 : 4311162 : det(1,1) = ten*det(1,1); det(2,1) = ten*det(2,1)
3455 : 4311162 : det(1,2) = det(1,2) - 1.0d0
3456 : 33686378 : go to 10
3457 : : 20 continue
3458 : : 30 continue
3459 : 30755395 : rabs = abs(det(1,1))+abs(det(2,1))
3460 [ + + ]: 30755395 : if (rabs .lt. ten) go to 40
3461 : 1380179 : det(1,1) = det(1,1)/ten; det(2,1) = det(2,1)/ten
3462 : 1380179 : det(1,2) = det(1,2) + 1.0d0
3463 : 30755395 : go to 30
3464 : 2394114 : 40 continue
3465 : : end do
3466 : : 60 continue
3467 : : 70 continue
3468 : : !
3469 : : ! compute inverse(u)
3470 : : !
3471 [ + + ]: 2394114 : if (mod(job,10) .eq. 0) go to 150
3472 [ + + ]: 16029689 : do 100 k = 1, n
3473 : : !a(k,k) = (1.0d0,0.0d0)/a(k,k)
3474 : : !t = -a(k,k)
3475 : : !call zscal(k-1,t,a(1,k),1)
3476 : 14800065 : rinv2 = 1.d0/(a(1,k,k)**2+a(2,k,k)**2)
3477 : 14800065 : a(1,k,k) = rinv2*a(1,k,k)
3478 : 14800065 : a(2,k,k) = -rinv2*a(2,k,k)
3479 : 14800065 : rk(1) = -a(1,k,k); rk(2) = -a(2,k,k)
3480 [ + + ]: 201539483 : do i=1,k-1
3481 : 186739418 : r(1)=a(1,i,k)
3482 : 186739418 : r(2)=a(2,i,k)
3483 : 186739418 : a(1,i,k)=rk(1)*r(1)-rk(2)*r(2)
3484 : 201539483 : a(2,i,k)=rk(1)*r(2)+rk(2)*r(1)
3485 : : end do
3486 : 14800065 : kp1 = k + 1
3487 [ + + ]: 14800065 : if (n .lt. kp1) go to 90
3488 [ + + ]: 200309859 : do 80 j = kp1, n
3489 : : !t = a(k,j)
3490 : : !a(k,j) = (0.0d0,0.0d0)
3491 : : !call zaxpy(k,t,a(1,k),1,a(1,j),1)
3492 : 186739418 : rkj(1) = a(1,k,j); rkj(2) = a(2,k,j)
3493 : 186739418 : a(1,k,j) = 0.d0; a(2,k,j) = 0.d0
3494 [ + + ]: 2298677559 : do i=1,k
3495 : 2111938141 : a(1,i,j)=rkj(1)*a(1,i,k)-rkj(2)*a(2,i,k)+a(1,i,j)
3496 : 2298677559 : a(2,i,j)=rkj(2)*a(1,i,k)+rkj(1)*a(2,i,k)+a(2,i,j)
3497 : : end do
3498 : 14800065 : 80 continue
3499 : : 90 continue
3500 : 1229624 : 100 continue
3501 : : do i=1,n
3502 : : end do
3503 : : !
3504 : : ! form inverse(u)*inverse(l)
3505 : : !
3506 : 1229624 : nm1 = n - 1
3507 [ + + ]: 1229624 : if (nm1 .lt. 1) go to 140
3508 [ + + ]: 14799945 : do 130 kb = 1, nm1
3509 : 13570441 : k = n - kb
3510 : 13570441 : kp1 = k + 1
3511 [ + + ]: 200309859 : do 110 i = kp1, n
3512 : 186739418 : work(1,i) = a(1,i,k); work(2,i) = a(2,i,k)
3513 : 186739418 : a(1,i,k) = 0.0d0; a(2,i,k) = 0.d0
3514 : 13570441 : 110 continue
3515 [ + + ]: 200309859 : do 120 j = kp1, n
3516 : 186739418 : r(1) = work(1,j); r(2) = work(2,j)
3517 : : !call zaxpy(n,t,a(1,j),1,a(1,k),1)
3518 [ + + ]: 6335814423 : do i=1,n
3519 : 6149075005 : a(1,i,k)=r(1)*a(1,i,j)-r(2)*a(2,i,j)+a(1,i,k)
3520 : 6335814423 : a(2,i,k)=r(2)*a(1,i,j)+r(1)*a(2,i,j)+a(2,i,k)
3521 : : end do
3522 : 13570441 : 120 continue
3523 : 13570441 : l = ipvt(k)
3524 [ + + ]: 13570441 : if (l .ne. k) then
3525 : : !call zswap(n,a(1,k),1,a(1,l),1)
3526 [ + + ]: 184583184 : do i=1,n
3527 : 178896595 : r(1) = a(1,i,k); r(2) = a(2,i,k)
3528 : 178896595 : a(1,i,k) = a(1,i,l); a(2,i,k) = a(2,i,l)
3529 : 184583184 : a(1,i,l) = r(1); a(2,i,l) = r(2)
3530 : : end do
3531 : : end if
3532 : 2394114 : 130 continue
3533 : : 140 continue
3534 : : 150 continue
3535 : :
3536 : 2394114 : end subroutine dzgedi
3537 : : !!***
3538 : :
3539 : : !!****f* m_hide_lapack/dzgefa
3540 : : !! NAME
3541 : : !! dzgefa
3542 : : !!
3543 : : !! FUNCTION
3544 : : !! This routine is the clone of zgefa.F90 using real*8 a(2) instead of complex*16
3545 : : !! for the purpose of ABINIT (2008,TD)
3546 : : !!
3547 : : !! NOTES
3548 : : !! This routine is deprecated, use lapack API
3549 : : !!
3550 : : !! SOURCE
3551 : :
3552 : 2394114 : subroutine dzgefa(a,lda,n,ipvt,info)
3553 : :
3554 : : use m_linalg_interfaces
3555 : :
3556 : : !Arguments
3557 : : integer :: lda,n,ipvt(n),info
3558 : : real*8 :: a(2,lda,n)
3559 : : !
3560 : : ! zgefa factors a complex*16 matrix by gaussian elimination.
3561 : : !
3562 : : ! dzgefa is usually called by zgeco, but it can be called
3563 : : ! directly with a saving in time if rcond is not needed.
3564 : : ! (time for zgeco) = (1 + 9/n)*(time for zgefa) .
3565 : : !
3566 : : ! on entry
3567 : : !
3568 : : ! a complex*16(lda, n)
3569 : : ! the matrix to be factored.
3570 : : !
3571 : : ! lda integer
3572 : : ! the leading dimension of the array a .
3573 : : !
3574 : : ! n integer
3575 : : ! the order of the matrix a .
3576 : : !
3577 : : ! on return
3578 : : !
3579 : : ! a an upper triangular matrix and the multipliers
3580 : : ! which were used to obtain it.
3581 : : ! the factorization can be written a = l*u where
3582 : : ! l is a product of permutation and unit lower
3583 : : ! triangular matrices and u is upper triangular.
3584 : : !
3585 : : ! ipvt integer(n)
3586 : : ! an integer vector of pivot indices.
3587 : : !
3588 : : ! info integer
3589 : : ! = 0 normal value.
3590 : : ! = k if u(k,k) .eq. 0.0 . this is not an error
3591 : : ! condition for this subroutine, but it does
3592 : : ! indicate that zgesl or zgedi will divide by zero
3593 : : ! if called. use rcond in zgeco for a reliable
3594 : : ! indication of singularity.
3595 : : !
3596 : : ! linpack. this version dated 08/14/78 .
3597 : : ! cleve moler, university of new mexico, argonne national lab.
3598 : : !
3599 : : ! subroutines and functions
3600 : : !
3601 : : ! internal variables
3602 : : !
3603 : : !Local variables
3604 : : real*8 :: r(2),rk(2),rlj(2)
3605 : : real*8 :: rinv2,rmax,rabs
3606 : : integer :: i,j,k,kp1,l,nm1
3607 : :
3608 : : !
3609 : : ! gaussian elimination with partial pivoting
3610 : : !
3611 : 2394114 : info = 0
3612 : 2394114 : nm1 = n - 1
3613 [ + + ]: 2394114 : if (nm1 .lt. 1) go to 70
3614 [ + + ]: 29495797 : do 60 k = 1, nm1
3615 : 27101803 : kp1 = k + 1
3616 : : !
3617 : : ! find l = pivot index
3618 : : !
3619 : : !l = izamax(n-k+1,a(k,k),1) + k - 1
3620 : 27101803 : rmax=0.d0
3621 : 27101803 : l=0
3622 [ + + ]: 429945310 : do i=k,n
3623 : 402843507 : rabs=abs(a(1,i,k))+abs(a(2,i,k))
3624 [ + + ]: 429945310 : if(rmax<=rabs) then
3625 : 44409324 : rmax=rabs
3626 : 44409324 : l=i
3627 : : end if
3628 : : end do
3629 : 27101803 : ipvt(k) = l
3630 : : !
3631 : : ! zero pivot implies this column already triangularized
3632 : : !
3633 [ + - ]: 27101803 : if (abs(a(1,l,k))+abs(a(2,l,k)) .eq. 0.0d0) go to 40
3634 : : !
3635 : : ! interchange if necessary
3636 : : !
3637 [ + + ]: 27101803 : if (l .eq. k) go to 10
3638 : 12609440 : r(1) = a(1,l,k); r(2) = a(2,l,k)
3639 : 12609440 : a(1,l,k) = a(1,k,k); a(2,l,k) = a(2,k,k)
3640 : 27101803 : a(1,k,k) = r(1); a(2,k,k) = r(2)
3641 : : 10 continue
3642 : : !
3643 : : ! compute multipliers
3644 : : !
3645 : 27101803 : rinv2 = 1.d0/(a(1,k,k)**2+a(2,k,k)**2)
3646 : 27101803 : rk(1) = -rinv2*a(1,k,k)
3647 : 27101803 : rk(2) = rinv2*a(2,k,k)
3648 : : !call zscal(n-k,t,a(k+1,k),1)
3649 [ + + ]: 402843507 : do i=k+1,n
3650 : 375741704 : r(1)=a(1,i,k)
3651 : 375741704 : r(2)=a(2,i,k)
3652 : 375741704 : a(1,i,k)=rk(1)*r(1)-rk(2)*r(2)
3653 : 402843507 : a(2,i,k)=rk(1)*r(2)+rk(2)*r(1)
3654 : : end do
3655 : : !
3656 : : ! row elimination with column indexing
3657 : : !
3658 [ + + ]: 402843507 : do j = kp1, n
3659 : 375741704 : rlj(1) = a(1,l,j); rlj(2) = a(2,l,j)
3660 [ + + ]: 375741704 : if (l .eq. k) go to 20
3661 : 171271954 : a(1,l,j) = a(1,k,j); a(2,l,j) = a(2,k,j)
3662 : 375741704 : a(1,k,j) = rlj(1); a(2,k,j) = rlj(2)
3663 : : 20 continue
3664 : : !call zaxpy(n-k,t,a(1,k+1,k),1,a(1,k+1,j),1)
3665 [ + + ]: 8534416049 : do i=k+1,n
3666 : 8131572542 : a(1,i,j)=rlj(1)*a(1,i,k)-rlj(2)*a(2,i,k)+a(1,i,j)
3667 : 8507314246 : a(2,i,j)=rlj(2)*a(1,i,k)+rlj(1)*a(2,i,k)+a(2,i,j)
3668 : : end do
3669 : : end do
3670 : 0 : go to 50
3671 : : 40 continue
3672 : 0 : info = k
3673 : : 50 continue
3674 : 2394114 : 60 continue
3675 : : 70 continue
3676 : 2394114 : ipvt(n) = n
3677 [ - + ]: 2394114 : if (abs(a(1,n,n))+abs(a(2,n,n)) .eq. 0.0d0) info = n
3678 : :
3679 : 2394114 : end subroutine dzgefa
3680 : : !!***
3681 : :
3682 : : !!****f* m_hide_lapack/test_xginv
3683 : : !! NAME
3684 : : !! test_xginv
3685 : : !!
3686 : : !! FUNCTION
3687 : :
3688 : 0 : subroutine test_xginv(msize,skinds,do_check,Tres,comm)
3689 : :
3690 : : !Arguments -----------------------------------
3691 : : !scalars
3692 : : integer,intent(in) :: msize,comm
3693 : : logical,intent(in) :: do_check
3694 : : character(len=*),intent(in) :: skinds
3695 : : type(latime_t),intent(out) :: Tres
3696 : : !arrays
3697 : : !complex(sp),allocatable :: cmat_spc(:,:)
3698 : : !complex(sp),allocatable :: cmat_spc_check(:,:)
3699 : 0 : complex(dp),allocatable :: cmat_dpc(:,:)
3700 : 0 : complex(dp),allocatable :: cmat_dpc_check(:,:)
3701 : :
3702 : : !Local variables-------------------------------
3703 : : integer :: ii
3704 : : real(dp) :: max_abserr
3705 : : ! *************************************************************************
3706 : :
3707 : : if (.FALSE.) write(std_out,*)skinds
3708 : :
3709 [ # # ]: 0 : if (do_check) then
3710 [ # # # # : 0 : ABI_MALLOC(cmat_dpc_check,(msize,msize))
# # # # #
# # # # #
# # ]
3711 [ # # # # ]: 0 : cmat_dpc_check = czero
3712 [ # # ]: 0 : do ii=1,msize
3713 : 0 : cmat_dpc_check(ii,ii) = cone
3714 : : end do
3715 : : !call xginv(cmat_dpc_check,msize,comm=xmpi_comm_self)
3716 : : end if
3717 : :
3718 [ # # # # : 0 : ABI_MALLOC(cmat_dpc,(msize,msize))
# # # # #
# # # # #
# # ]
3719 [ # # ]: 0 : do ii=1,msize
3720 : 0 : cmat_dpc(ii,ii) = cone
3721 : : end do
3722 : :
3723 : 0 : call cwtime(Tres%ctime,Tres%wtime,Tres%gflops,"start")
3724 : :
3725 : 0 : call xginv(cmat_dpc,msize,comm)
3726 : :
3727 : 0 : call cwtime(Tres%ctime,Tres%wtime,Tres%gflops,"stop")
3728 : 0 : Tres%testname = 'test_xginv'
3729 : 0 : Tres%msize = msize
3730 : :
3731 : 0 : max_abserr = -one
3732 [ # # ]: 0 : if (do_check) then
3733 [ # # # # : 0 : max_abserr = MAXVAL( ABS(cmat_dpc - cmat_dpc_check) )
# # # # #
# # # ]
3734 : : end if
3735 : 0 : Tres%max_abserr = max_abserr
3736 : :
3737 [ # # ]: 0 : ABI_FREE(cmat_dpc)
3738 [ # # ]: 0 : ABI_SFREE(cmat_dpc_check)
3739 : :
3740 : 0 : end subroutine test_xginv
3741 : : !!***
3742 : :
3743 : : !!****f* m_cgtools/xhesv_cplex
3744 : : !! NAME
3745 : : !! xhesv_cplex
3746 : : !!
3747 : : !! FUNCTION
3748 : : !! ZHESV computes the solution to a complex (real) system of linear equations
3749 : : !! A * X = B,
3750 : : !!
3751 : : !! where A is an N-by-N Hermitian matrix and X and B are N-by-NRHS matrices.
3752 : : !! The value of cplex (1 or 2) defines whether we have a complex Hermitian or real symmetric matrix
3753 : : !!
3754 : : !! The diagonal pivoting method is used to factor A as
3755 : : !! A = U * D * U**H, if UPLO = 'U', or
3756 : : !! A = L * D * L**H, if UPLO = 'L',
3757 : : !!
3758 : : !! where U (or L) is a product of permutation and unit upper (lower)
3759 : : !! triangular matrices, and D is Hermitian and block diagonal with
3760 : : !! 1-by-1 and 2-by-2 diagonal blocks. The factored form of A is then
3761 : : !! used to solve the system of equations A * X = B.
3762 : : !!
3763 : : !! INPUTS
3764 : : !!
3765 : : !![in] UPLO
3766 : : !! UPLO is CHARACTER*1
3767 : : !! = 'U': Upper triangle of A is stored;
3768 : : !! = 'L': Lower triangle of A is stored.
3769 : : !![in] N
3770 : : !! N is INTEGER
3771 : : !! The number of linear equations, i.e., the order of the
3772 : : !! matrix A. N >= 0.
3773 : : !![in] NRHS
3774 : : !! NRHS is INTEGER
3775 : : !! The number of right hand sides, i.e., the number of columns
3776 : : !! of the matrix B. NRHS >= 0.
3777 : : !![in,out] A
3778 : : !! A is COMPLEX*16 array, dimension (LDA,N)
3779 : : !! On entry, the Hermitian matrix A. If UPLO = 'U', the leading
3780 : : !! N-by-N upper triangular part of A contains the upper
3781 : : !! triangular part of the matrix A, and the strictly lower
3782 : : !! triangular part of A is not referenced. If UPLO = 'L', the
3783 : : !! leading N-by-N lower triangular part of A contains the lower
3784 : : !! triangular part of the matrix A, and the strictly upper
3785 : : !! triangular part of A is not referenced.
3786 : : !!
3787 : : !! On exit, if INFO = 0, the block diagonal matrix D and the
3788 : : !! multipliers used to obtain the factor U or L from the
3789 : : !! factorization A = U*D*U**H or A = L*D*L**H as computed by
3790 : : !! ZHETRF.
3791 : : !![in,out] B
3792 : : !! B is COMPLEX*16 array, dimension (LDB,NRHS)
3793 : : !! On entry, the N-by-NRHS right hand side matrix B.
3794 : : !! On exit, if INFO = 0, the N-by-NRHS solution matrix X.
3795 : : !![out] INFO
3796 : : !! INFO is INTEGER
3797 : : !! = 0: successful exit
3798 : : !! < 0: if INFO = -i, the i-th argument had an illegal value
3799 : : !! > 0: if INFO = i, D(i,i) is exactly zero. The factorization
3800 : : !! has been completed, but the block diagonal matrix D is
3801 : : !! exactly singular, so the solution could not be computed.
3802 : : !! SOURCE
3803 : :
3804 : 9576 : subroutine xhesv_cplex(UPLO, cplex, N, NRHS, A, B, msg, info)
3805 : :
3806 : : !Arguments ------------------------------------
3807 : : character(len=1),intent(in) :: UPLO
3808 : : integer,intent(in) :: cplex, N, NRHS
3809 : : real(dp),intent(inout) :: A(cplex, N, N)
3810 : : real(dp),intent(inout) :: B(cplex, N, NRHS)
3811 : : character(len=*),intent(out) :: msg
3812 : : integer,intent(out) :: info
3813 : :
3814 : : !Local variables ------------------------------
3815 : : !scalars
3816 : : integer :: lwork, lda, ldb
3817 : 9576 : integer,allocatable :: ipiv(:)
3818 : 9576 : real(dp),allocatable :: work(:,:)
3819 : : !************************************************************************
3820 : :
3821 [ + - + + : 17148 : if (all(cplex /= [1, 2])) then
- + ]
3822 : 0 : write(msg,'(a,i0)')" Wrong value for cplex: ",cplex
3823 : 0 : info = 1; return
3824 : : end if
3825 : :
3826 : 9576 : lda = N; ldb = N
3827 : :
3828 [ - + + - : 28728 : ABI_MALLOC(ipiv, (N))
+ - + - -
+ - + ]
3829 [ - + + - : 28728 : ABI_MALLOC(work, (cplex, 1))
+ - + - -
+ - + ]
3830 : 9576 : lwork = -1
3831 : :
3832 [ + + ]: 9576 : if (cplex == 2) then
3833 : : ! Complex version
3834 : 7572 : call zhesv(uplo, N, NRHS, A, lda, ipiv, B, ldb, work, lwork, info)
3835 : 7572 : lwork = int(work(1, 1))
3836 [ - + + - : 30288 : ABI_REMALLOC(work, (cplex, lwork))
- + + - +
- + - - +
- + - + ]
3837 : :
3838 : 7572 : call zhesv(uplo, N, NRHS, A, lda, ipiv, B, ldb, work, lwork, info)
3839 : : else
3840 : : ! Read version
3841 : 2004 : call dsysv(uplo, N, NRHS, A, lda, ipiv, B, ldb, work, lwork, info)
3842 : 2004 : lwork = int(work(1, 1))
3843 [ - + + - : 8016 : ABI_REMALLOC(work, (cplex, lwork))
- + + - +
- + - - +
- + - + ]
3844 : :
3845 : 2004 : call dsysv(uplo, N, NRHS, A, lda, ipiv, B, ldb, work, lwork, info)
3846 : : end if
3847 : :
3848 [ - + ]: 9576 : ABI_FREE(ipiv)
3849 [ - + ]: 9576 : ABI_FREE(work)
3850 : :
3851 [ - + ]: 9576 : if (info < 0) then
3852 : 0 : write(msg,'(a,i0,a)')" The ",-info,"-th argument of ZPOTRI had an illegal value."
3853 [ - + ]: 9576 : else if (info > 0) then
3854 : : write(msg, "(a,i0,4a)") &
3855 : 0 : " D(i,i) is exactly zero for i= ", info, ch10, &
3856 : 0 : "The factorization has been completed, but the block diagonal matrix D is ", ch10, &
3857 : 0 : "exactly singular, so the solution could not be computed."
3858 : : end if
3859 : :
3860 : : end subroutine xhesv_cplex
3861 : : !!***
3862 : :
3863 : 0 : end module m_hide_lapack
3864 : : !!***
|
