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 52387 : 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 52387 : real(dp),allocatable :: rwork(:)
332 52387 : 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 52387 : use_scalapack=.FALSE.
342 52387 : 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 52387 : lwork = MAX(1,2*n-1)
353 157161 : ABI_MALLOC(work, (lwork))
354 157161 : ABI_MALLOC(rwork, (MAX(1,3*n-2)))
355 :
356 52387 : call ZHEEV(jobz,uplo,n,a,n,w,work,lwork,rwork,info)
357 :
358 52387 : 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 52387 : 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 52387 : ABI_FREE(rwork)
371 52387 : 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 8463483 : 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 8463483 : integer,allocatable :: ipiv(:)
2708 8463483 : 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 8463483 : use_scalapack=.FALSE.
2717 8463483 : if (PRESENT(comm)) then
2718 5977 : 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 25390449 : ABI_MALLOC(ipiv, (n))
2727 8463483 : call ZGETRF(n,n,a,n,ipiv,info) ! P* L* U Factorization.
2728 :
2729 8463483 : 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 8463483 : 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 8463483 : lwork=MAX(1,n)
2744 25390449 : ABI_MALLOC(work,(lwork))
2745 :
2746 8463483 : call ZGETRI(n,a,n,ipiv,work,lwork,info) ! Invert U and then compute inv(A)
2747 :
2748 8463483 : 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 8463483 : 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 8463483 : ABI_FREE(ipiv)
2761 8463483 : 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 : !!***
|
