Line data Source code
1 : !!****m* ABINIT/m_gwls_QR_factorization
2 : !! NAME
3 : !! m_gwls_QR_factorization
4 : !!
5 : !! FUNCTION
6 : !! .
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 2009-2026 ABINIT group (JLJ, BR, MC)
10 : !! This file is distributed under the terms of the
11 : !! GNU General Public License, see ~abinit/COPYING
12 : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : !!
14 : !! SOURCE
15 :
16 : #if defined HAVE_CONFIG_H
17 : #include "config.h"
18 : #endif
19 :
20 : #include "abi_common.h"
21 :
22 :
23 :
24 : module m_gwls_QR_factorization
25 : !----------------------------------------------------------------------------------------------------
26 : ! This module implements the QR factorization using various algorithms, for the specific
27 : ! data distribution corresponding to FFT parallelism.
28 : !
29 : ! There are standard routines which do this (lapack, scalapack, etc), however it is complicated
30 : ! to get scalapack to run properly in parallel. Implementing ourselves is the shortest path to
31 : ! a working solution.
32 : !----------------------------------------------------------------------------------------------------
33 : !local modules
34 : use m_gwls_utility
35 : use m_gwls_TimingLog
36 : use m_gwls_wf
37 : use m_gwls_hamiltonian
38 :
39 : !abinit modules
40 : use defs_basis
41 : use defs_wvltypes
42 : use m_abicore
43 : use m_xmpi
44 : use m_errors
45 :
46 : use defs_abitypes, only : MPI_type
47 : use m_io_tools, only : get_unit
48 : use m_time, only : timab
49 :
50 :
51 : implicit none
52 : save
53 : private
54 : !!***
55 :
56 : logical, private :: debug = .false.
57 : !!***
58 :
59 : public :: extract_QR, extract_SVD
60 : !!***
61 :
62 : contains
63 :
64 : !!****f* m_hamiltonian/extract_QR
65 : !! NAME
66 : !! extract_QR
67 : !!
68 : !! FUNCTION
69 : !! .
70 : !!
71 : !! INPUTS
72 : !!
73 : !! OUTPUT
74 : !!
75 : !! SOURCE
76 :
77 2105 : subroutine extract_QR(mpi_communicator,Hsize,Xsize,Xmatrix,Rmatrix)
78 : !--------------------------------------------------------------------------
79 : ! This function computes the QR factorization:
80 : !
81 : ! X = Q . R
82 : !
83 : ! in order to extract the matrix of orthonormal vectors Q and
84 : ! the R matrix.
85 : !
86 : ! On output, the matrix X is replaced by Q.
87 : !
88 : ! If the code is running with only one processor, this routine
89 : ! simply invokes extract_QR_serial, which wraps standard Lapack routines.
90 : ! If we are running in MPI parallel, the serial Lapack routines no
91 : ! longer work (and understanding scalapack is too complicated right now).
92 : ! Thus, in that case, this routine implements some old school Gram-Schmidt
93 : ! algorithm.
94 : !--------------------------------------------------------------------------
95 :
96 : integer, intent(in) :: Hsize, Xsize, mpi_communicator
97 : complex(dp),intent(inout) :: Xmatrix(Hsize,Xsize)
98 :
99 : complex(dp), intent(out),optional :: Rmatrix(Xsize,Xsize)
100 :
101 : ! local variables
102 :
103 : real(dp) :: tsec(2)
104 : integer :: GWLS_TIMAB, OPTION_TIMAB
105 :
106 : ! *************************************************************************
107 :
108 :
109 :
110 2105 : GWLS_TIMAB = 1519
111 2105 : OPTION_TIMAB = 1
112 2105 : call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
113 :
114 :
115 : !--------------------------------------------------------------------------------
116 : ! Implement Gram-Schmidt.
117 : !--------------------------------------------------------------------------------
118 2105 : call extract_QR_Householder(mpi_communicator,Hsize,Xsize,Xmatrix,Rmatrix)
119 :
120 2105 : OPTION_TIMAB = 2
121 2105 : call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
122 :
123 2105 : end subroutine extract_QR
124 : !!***
125 :
126 :
127 : !!****f* m_hamiltonian/extract_SVD
128 : !! NAME
129 : !! extract_SVD
130 : !!
131 : !! FUNCTION
132 : !! .
133 : !!
134 : !! INPUTS
135 : !!
136 : !! OUTPUT
137 : !!
138 : !! SOURCE
139 :
140 28 : subroutine extract_SVD(mpi_communicator, Hsize,lsolutions_max,svd_matrix,svd_values)
141 : !--------------------------------------------------------------------------
142 : ! This function computes the singular value decomposition
143 : !
144 : ! X = U . SIGMA . V^dagger
145 : !
146 : ! More specifically, the matrix U of orthonormal vectors and SIGMA
147 : ! the eigenvalues are returned.
148 : !
149 : ! different algorithms are used, depending on parallelisation scheme.
150 : !--------------------------------------------------------------------------
151 :
152 : integer, intent(in) :: mpi_communicator
153 : integer, intent(in) :: Hsize, lsolutions_max
154 : complex(dp), intent(inout) :: svd_matrix(Hsize,lsolutions_max)
155 : real(dp), intent(out) :: svd_values(lsolutions_max)
156 :
157 :
158 28 : complex(dp), allocatable :: Rmatrix(:,:)
159 28 : complex(dp), allocatable :: svd_tmp(:,:)
160 :
161 : real(dp) :: tsec(2)
162 : integer :: GWLS_TIMAB, OPTION_TIMAB
163 :
164 : ! *************************************************************************
165 :
166 :
167 :
168 28 : GWLS_TIMAB = 1520
169 28 : OPTION_TIMAB = 1
170 28 : call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
171 :
172 :
173 : !if ( mpi_enreg%nproc_fft ==1 ) then
174 : if ( .false. ) then
175 :
176 : call extract_SVD_lapack(Hsize,lsolutions_max,svd_matrix,svd_values)
177 :
178 : else
179 :
180 112 : ABI_MALLOC(Rmatrix,(lsolutions_max,lsolutions_max))
181 :
182 : ! perform QR first
183 28 : call extract_QR(mpi_communicator, Hsize,lsolutions_max,svd_matrix,Rmatrix)
184 :
185 : ! perform SVD on the much smaller Rmatrix!
186 28 : call extract_SVD_lapack(lsolutions_max,lsolutions_max,Rmatrix,svd_values)
187 :
188 112 : ABI_MALLOC(svd_tmp,(Hsize,lsolutions_max))
189 :
190 : ! Rmatrix is overwritten with U matrix from SVD. Update the svd_matrix
191 : call ZGEMM( 'N', & ! Leave first array as is
192 : 'N', & ! Leave second array as is
193 : Hsize, & ! the number of rows of the matrix op( A )
194 : lsolutions_max, & ! the number of columns of the matrix op( B )
195 : lsolutions_max, & ! the number of columns of the matrix op( A ) == rows of matrix op( B )
196 : cmplx_1, & ! alpha constant
197 : svd_matrix, & ! matrix A
198 : Hsize, & ! LDA
199 : Rmatrix, & ! matrix B
200 : lsolutions_max, & ! LDB
201 : cmplx_0, & ! beta constant
202 : svd_tmp, & ! matrix C
203 28 : Hsize) ! LDC
204 :
205 66268 : svd_matrix(:,:) = svd_tmp(:,:)
206 :
207 28 : ABI_FREE(svd_tmp)
208 28 : ABI_FREE(Rmatrix)
209 : end if
210 28 : OPTION_TIMAB = 2
211 28 : call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
212 :
213 :
214 28 : end subroutine extract_SVD
215 : !!***
216 :
217 : !!****f* m_hamiltonian/extract_SVD_lapack
218 : !! NAME
219 : !! extract_SVD_lapack
220 : !!
221 : !! FUNCTION
222 : !! .
223 : !!
224 : !! INPUTS
225 : !!
226 : !! OUTPUT
227 : !!
228 : !! SOURCE
229 :
230 28 : subroutine extract_SVD_lapack(Hsize,lsolutions_max,svd_matrix,svd_values)
231 : !--------------------------------------------------------------------------
232 : ! This function computes the singular value decomposition
233 : ! using lapack routines. This is not appropriate in MPI parallel!
234 : !
235 : !
236 : !--------------------------------------------------------------------------
237 :
238 :
239 : integer, intent(in) :: Hsize, lsolutions_max
240 : complex(dp), intent(inout) :: svd_matrix(Hsize,lsolutions_max)
241 : real(dp), intent(out) :: svd_values(lsolutions_max)
242 :
243 :
244 :
245 : integer :: info_zgesvd
246 : integer :: lwork_svd
247 28 : complex(dp), allocatable :: work_svd(:)
248 28 : complex(dp), allocatable :: svd_U(:,:), svd_V(:,:)
249 28 : real (dp ), allocatable :: rwork_svd(:)
250 :
251 : integer :: debug_unit
252 : character(50) :: debug_filename
253 :
254 : ! *************************************************************************
255 :
256 :
257 :
258 :
259 : ! allocate arrays for the svd
260 28 : ABI_MALLOC(svd_U ,(1,1))
261 28 : ABI_MALLOC(svd_V ,(1,1))
262 :
263 :
264 : ! DIMENSION QUERRY for singluar decomposition problem
265 :
266 84 : ABI_MALLOC(rwork_svd ,(5*min(Hsize,lsolutions_max)))
267 28 : ABI_MALLOC(work_svd,(1))
268 28 : lwork_svd = -1
269 :
270 : call zgesvd('O', & ! The first min(m,n) columns of U (the left singular vectors) are overwritten on the array A;
271 : 'N', & ! no column vectors of V are computed
272 : Hsize, & ! number of rows of the matrix
273 : lsolutions_max, & ! number of columns of the matrix
274 : svd_matrix, & ! matrix to be decomposed
275 : Hsize, & ! LDA
276 : svd_values, & ! singular values
277 : svd_U, & ! dummy U; not referenced
278 : 1, & ! size of U
279 : svd_V, & ! dummy V; not referenced
280 : 1, & ! size of V
281 : work_svd, & ! work array
282 : lwork_svd, & ! size of work array
283 : rwork_svd, & ! work array
284 28 : info_zgesvd )
285 :
286 28 : if ( info_zgesvd /= 0) then
287 0 : debug_unit = get_unit()
288 0 : write(debug_filename,'(A,I4.4,A)') 'LAPACK_DEBUG_PROC=',mpi_enreg%me,'.log'
289 :
290 0 : open(debug_unit,file=trim(debug_filename),status='unknown')
291 :
292 0 : write(debug_unit,'(A)') '*********************************************************************************************'
293 0 : write(debug_unit,'(A,I4,A)') '* ERROR: info = ',info_zgesvd,' in ZGESVD(1), gwls_QR_factorization'
294 0 : write(debug_unit,'(A)') '*********************************************************************************************'
295 :
296 0 : close(debug_unit)
297 :
298 : end if
299 :
300 :
301 :
302 :
303 :
304 :
305 28 : lwork_svd = nint(dble(work_svd(1)))
306 :
307 28 : ABI_FREE(work_svd)
308 :
309 84 : ABI_MALLOC(work_svd,(lwork_svd))
310 :
311 : ! computation run
312 :
313 : call zgesvd('O', & ! The first min(m,n) columns of U (the left singular vectors) are overwritten on the array A;
314 : 'N', & ! no column vectors of V are computed
315 : Hsize, & ! number of rows of the matrix
316 : lsolutions_max, & ! number of columns of the matrix
317 : svd_matrix, & ! matrix to be decomposed
318 : Hsize, & ! LDA
319 : svd_values, & ! singular values
320 : svd_U, & ! dummy U; not referenced
321 : 1, & ! size of U
322 : svd_V, & ! dummy V; not referenced
323 : 1, & ! size of V
324 : work_svd, & ! work array
325 : lwork_svd, & ! size of work array
326 : rwork_svd, & ! work array
327 28 : info_zgesvd )
328 :
329 28 : if ( info_zgesvd /= 0) then
330 0 : debug_unit = get_unit()
331 0 : write(debug_filename,'(A,I4.4,A)') 'LAPACK_DEBUG_PROC=',mpi_enreg%me,'.log'
332 :
333 0 : open(debug_unit,file=trim(debug_filename),status='unknown')
334 :
335 0 : write(debug_unit,'(A)') '*********************************************************************************************'
336 0 : write(debug_unit,'(A,I4,A)') '* ERROR: info = ',info_zgesvd,' in ZGESVD(2), gwls_QR_factorization'
337 0 : write(debug_unit,'(A)') '*********************************************************************************************'
338 :
339 0 : close(debug_unit)
340 :
341 : end if
342 :
343 :
344 :
345 :
346 28 : ABI_FREE(work_svd)
347 28 : ABI_FREE(rwork_svd)
348 28 : ABI_FREE(svd_U)
349 28 : ABI_FREE(svd_V)
350 :
351 28 : end subroutine extract_SVD_lapack
352 : !!***
353 :
354 :
355 :
356 :
357 :
358 : !!****f* m_hamiltonian/extract_QR_Householder
359 : !! NAME
360 : !! extract_QR_Householder
361 : !!
362 : !! FUNCTION
363 : !! .
364 : !!
365 : !! INPUTS
366 : !!
367 : !! OUTPUT
368 : !!
369 : !! SOURCE
370 :
371 1932 : subroutine extract_QR_Householder(mpi_communicator,Hsize,Xsize,Xmatrix,Rmatrix)
372 : !--------------------------------------------------------------------------
373 : ! This function computes the QR factorization:
374 : !
375 : ! X = Q . R
376 : !
377 : ! in order to extract the matrix of orthonormal vectors Q and
378 : ! the R matrix.
379 : !
380 : ! On output, the matrix X is replaced by Q.
381 : !
382 : ! This routine uses Householder operations to generate Q and R.
383 : ! Special attention is given to the fact that the matrix may be
384 : ! distributed across processors and MPI communication is necessary.
385 : !
386 : !--------------------------------------------------------------------------
387 :
388 : integer, intent(in) :: Hsize, Xsize, mpi_communicator
389 : complex(dp),intent(inout) :: Xmatrix(Hsize,Xsize)
390 :
391 : complex(dp), intent(out),optional :: Rmatrix(Xsize,Xsize)
392 :
393 : ! local variables
394 : integer :: numbrer_of_plane_waves
395 :
396 : integer :: io_unit
397 : integer :: ierr
398 : character(50) :: filename
399 : logical :: file_exists
400 :
401 : integer,save :: counter = 0
402 :
403 : integer :: i, j, l_local
404 : integer :: l1, l2
405 :
406 2105 : integer, allocatable :: nproc_array(:)
407 :
408 2105 : complex(dp), allocatable :: Qinternal(:,:)
409 2105 : complex(dp), allocatable :: Rinternal(:,:)
410 2105 : complex(dp), allocatable :: vj(:)
411 2105 : complex(dp), allocatable :: A_matrix(:,:)
412 2105 : complex(dp), allocatable :: V_matrix(:,:)
413 2105 : complex(dp), allocatable :: list_beta(:)
414 : complex(dp) :: cmplx_value
415 : real (dp ) :: real_value
416 :
417 : complex(dp) :: norm_x
418 : complex(dp) :: phase
419 2105 : complex(dp), allocatable :: error(:,:)
420 2105 : complex(dp), allocatable :: coeff(:)
421 :
422 :
423 : integer :: mpi_rank
424 : integer :: mpi_nproc
425 : logical :: head_node
426 :
427 : ! *************************************************************************
428 :
429 :
430 : !--------------------------------------------------------------------------------
431 : ! Implement Householder algorithm, in parallel
432 : !--------------------------------------------------------------------------------
433 2105 : mpi_nproc = xmpi_comm_size(mpi_communicator)
434 :
435 : ! extract the rank of this processor in the communicator
436 2105 : mpi_rank = xmpi_comm_rank(mpi_communicator)
437 :
438 : ! only head node will write!
439 2105 : head_node = mpi_rank == 0
440 :
441 :
442 : !--------------------------------------------------------------------------------
443 : ! Open a log file for the output of extract_QR
444 : !--------------------------------------------------------------------------------
445 2105 : if (debug .and. head_node ) then
446 :
447 0 : io_unit = get_unit()
448 0 : write(filename,'(A,I0.4,A)') "extract_QR_",mpi_rank,".log"
449 0 : inquire(file=trim(filename),exist=file_exists)
450 :
451 0 : if (file_exists) then
452 0 : open(io_unit,file=trim(filename),position='append',status=files_status_old)
453 : else
454 0 : open(io_unit,file=trim(filename),status=files_status_new)
455 0 : write(io_unit,10) "#======================================================================================="
456 0 : write(io_unit,10) "# "
457 0 : write(io_unit,10) "# This file contains information regarding the QR factorization, from extract_QR "
458 0 : write(io_unit,25) "# The algorithm is running in MPI parallel with ",mpi_nproc," processors"
459 0 : write(io_unit,10) "# "
460 0 : write(io_unit,10) "#======================================================================================="
461 : end if
462 :
463 0 : counter = counter + 1
464 :
465 0 : write(io_unit,10) "# "
466 0 : write(io_unit,11) "# Call # ", counter
467 0 : write(io_unit,10) "# "
468 0 : write(io_unit,11) "# Hsize = ",Hsize
469 0 : write(io_unit,11) "# Xsize = ",Xsize
470 0 : write(io_unit,13) "# Rmatrix present? = ",present(Rmatrix)
471 :
472 :
473 : end if
474 :
475 : !--------------------------------------------------------------------------------
476 : ! Get the number of plane waves on every processor
477 : !--------------------------------------------------------------------------------
478 6315 : ABI_MALLOC(nproc_array,(mpi_nproc))
479 :
480 5426 : nproc_array = 0
481 :
482 2105 : numbrer_of_plane_waves = Hsize ! do this to avoid "intent" problems
483 2105 : call xmpi_allgather(numbrer_of_plane_waves, nproc_array, mpi_communicator, ierr)
484 :
485 :
486 : !--------------------------------------------------------------------------------
487 : ! Get the offset for each processor
488 : !
489 : ! The global index is then given by
490 : ! I_{global} = nproc_array(1+rank)+i_{local}
491 : !
492 : ! similarly, the local index is given by
493 : ! i_{local} = I_{global} - nproc_array(1+rank)
494 : ! which is only meaningful if 1 <= i_{local} <= Hsize
495 : !--------------------------------------------------------------------------------
496 :
497 5426 : do j = mpi_nproc, 1, -1
498 3321 : nproc_array(j) = 0
499 7026 : do i = 1, j-1
500 4921 : nproc_array(j) = nproc_array(j)+nproc_array(i)
501 : end do
502 : end do
503 :
504 :
505 : !--------------------------------------------------------------------------------
506 : ! Act on the A matrix, following the book by Golub (more or less ;) )
507 : !
508 : !--------------------------------------------------------------------------------
509 :
510 8420 : ABI_MALLOC(A_matrix, (Hsize,Xsize))
511 6315 : ABI_MALLOC(V_matrix, (Hsize,Xsize))
512 6315 : ABI_MALLOC(list_beta, (Xsize))
513 4210 : ABI_MALLOC(coeff , (Xsize))
514 :
515 486061 : A_matrix(:,:) = Xmatrix(:,:)
516 486061 : V_matrix(:,:) = cmplx_0
517 4700 : list_beta(:) = cmplx_0
518 :
519 :
520 6315 : ABI_MALLOC(vj, (Hsize))
521 :
522 4700 : do j = 1, Xsize
523 :
524 : ! Store xj in vj, for now
525 483956 : vj(:) = A_matrix(:,j)
526 :
527 :
528 2595 : if (j > 1) then
529 : !------------------------------------------
530 : ! set the array to zero all the way to j-1
531 : !------------------------------------------
532 490 : l_local = j-1-nproc_array(1+mpi_rank)
533 :
534 490 : if ( l_local > Hsize) then
535 0 : vj(:) = cmplx_0
536 490 : else if ( l_local <= Hsize .and. l_local >= 1) then
537 2260 : vj(1:l_local) = cmplx_0
538 : end if
539 :
540 : end if
541 :
542 : ! compute the norm of x
543 483956 : norm_x = sum(conjg(vj(:))*vj(:))
544 2595 : call xmpi_sum(norm_x,mpi_communicator,ierr) ! sum on all processors in communicator
545 2595 : norm_x = sqrt(norm_x)
546 :
547 :
548 4700 : if (abs(norm_x) > tol14) then
549 : ! if |x| ~ 0, there is nothing to do! the column in A is full of zeros.
550 :
551 : ! find the j^th component of x
552 2595 : l_local = j-nproc_array(1+mpi_rank)
553 :
554 : ! update vj, on the right processor!
555 2595 : if ( l_local <= Hsize .and. l_local >= 1) then
556 :
557 1873 : phase = vj(l_local)
558 :
559 1873 : if (abs(phase) < tol14) then
560 : phase = cmplx_1
561 : else
562 1829 : phase = phase/abs(phase)
563 : end if
564 :
565 1873 : vj(l_local) = vj(l_local) + phase*norm_x
566 :
567 : end if
568 :
569 : !compute beta
570 483956 : cmplx_value = sum(conjg(vj(:))*vj(:))
571 2595 : call xmpi_sum(cmplx_value,mpi_communicator,ierr) ! sum on all processors
572 :
573 2595 : list_beta(j) = 2.0_dp/cmplx_value
574 :
575 : ! store v for later use; this is less efficient than storing it in the null part of A,
576 : ! but it is less of a pain to implement. Feel free to clean this up.
577 483956 : V_matrix(:,j) = vj(:)
578 :
579 :
580 : ! Update the A matrix
581 :
582 : ! Compute v^dagger . A
583 : !call ZGEMM( 'C', & ! Hermitian conjugate the first array
584 : ! 'N', & ! Leave second array as is
585 : ! 1, & ! the number of rows of the matrix op( A )
586 : ! Xsize-j+1, & ! the number of columns of the matrix op( B )
587 : ! Hsize, & ! the number of columns of the matrix op( A ) == rows of matrix op( B )
588 : ! cmplx_1, & ! alpha constant
589 : ! vj, & ! matrix A
590 : ! Hsize, & ! LDA
591 : ! A_matrix(:,j:Xsize), & ! matrix B
592 : ! Hsize, & ! LDB
593 : ! cmplx_0, & ! beta constant
594 : ! coeff(:,j:Xsize), & ! matrix C
595 : ! 1) ! LDC
596 8655 : do i = j, Xsize
597 998876 : coeff(i) = sum(conjg(vj)*A_matrix(:,i))
598 : end do
599 2595 : call xmpi_sum(coeff,mpi_communicator,ierr) ! sum on all processors in the communicator
600 :
601 : ! update A
602 8655 : do i = j, Xsize
603 998876 : A_matrix(:,i) = A_matrix(:,i) - list_beta(j)*coeff(i)*vj(:)
604 : end do
605 :
606 : end if
607 :
608 : end do
609 :
610 : !--------------------------------------------------------------------------------
611 : ! Extract the R matrix
612 : !
613 : !--------------------------------------------------------------------------------
614 :
615 8420 : ABI_MALLOC(Rinternal,(Xsize,Xsize))
616 :
617 14225 : Rinternal = cmplx_0
618 :
619 :
620 4700 : do j = 1, Xsize
621 10760 : do i = 1, j
622 6060 : l_local = i-nproc_array(1+mpi_rank)
623 :
624 8655 : if ( l_local <= Hsize .and. l_local >= 1) then
625 3853 : Rinternal(i,j) = A_matrix(l_local,j)
626 : end if
627 :
628 : end do ! i
629 : end do ! j
630 :
631 2105 : call xmpi_sum(Rinternal,mpi_communicator,ierr) ! sum on all processors
632 :
633 :
634 :
635 : !--------------------------------------------------------------------------------
636 : ! Extract the Q matrix
637 : !
638 : !--------------------------------------------------------------------------------
639 :
640 6315 : ABI_MALLOC( Qinternal, (Hsize,Xsize))
641 :
642 : ! initialize Q to the identity in the top corner
643 486061 : Qinternal = cmplx_0
644 :
645 4700 : do j = 1, Xsize
646 2595 : l_local = j-nproc_array(1+mpi_rank)
647 :
648 4700 : if ( l_local <= Hsize .and. l_local >= 1 ) then
649 1873 : Qinternal(l_local,j) = cmplx_1
650 : end if
651 :
652 : end do ! j
653 :
654 :
655 : ! Build Q interatively
656 4700 : do j = Xsize,1, -1
657 :
658 483956 : vj(:) = V_matrix(:,j)
659 :
660 :
661 : ! Update the A matrix
662 :
663 : ! Compute v^dagger . A
664 : !call ZGEMM( 'C', & ! Hermitian conjugate the first array
665 : ! 'N', & ! Leave second array as is
666 : ! 1, & ! the number of rows of the matrix op( A )
667 : ! Xsize, & ! the number of columns of the matrix op( B )
668 : ! Hsize, & ! the number of columns of the matrix op( A ) == rows of matrix op( B )
669 : ! cmplx_1, & ! alpha constant
670 : ! vj, & ! matrix A
671 : ! Hsize, & ! LDA
672 : ! Qinternal, & ! matrix B
673 : ! Hsize, & ! LDB
674 : ! cmplx_0, & ! beta constant
675 : ! coeff, & ! matrix C
676 : ! 1) ! LDC
677 :
678 12120 : do i = 1, Xsize
679 1511201 : coeff(i) = sum(conjg(vj)*Qinternal(:,i))
680 : end do
681 2595 : call xmpi_sum(coeff,mpi_communicator,ierr) ! sum on all processors in communicator
682 :
683 :
684 : ! update Q
685 14225 : do i = 1, Xsize
686 1511201 : Qinternal(:,i) = Qinternal(:,i) - list_beta(j)*coeff(i)*vj(:)
687 : end do
688 :
689 : end do ! j
690 :
691 :
692 :
693 : ! clean up
694 2105 : ABI_FREE(V_matrix)
695 2105 : ABI_FREE(coeff)
696 2105 : ABI_FREE(vj)
697 :
698 : !--------------------------------------------------------------------------------
699 : ! Do some debug, if requested
700 : !
701 : !--------------------------------------------------------------------------------
702 :
703 2105 : if (debug ) then
704 :
705 0 : if ( head_node ) then
706 :
707 0 : write(io_unit,20) "# nproc_array = ",nproc_array
708 0 : flush(io_unit)
709 :
710 0 : write(io_unit,40) "# list_beta = ",real(list_beta)
711 0 : flush(io_unit)
712 : end if
713 :
714 :
715 0 : ABI_MALLOC(error,(Xsize,Xsize))
716 :
717 0 : error = cmplx_0
718 :
719 0 : do l2=1,Xsize
720 0 : error(l2,l2) = error(l2,l2) - cmplx_1
721 0 : do l1=1,Xsize
722 :
723 0 : cmplx_value = complex_vector_product(Qinternal(:,l1),Qinternal(:,l2),Hsize)
724 0 : call xmpi_sum(cmplx_value,mpi_communicator,ierr) ! sum on all processors working on FFT!
725 :
726 0 : error(l1,l2) = error(l1,l2)+cmplx_value
727 :
728 : end do
729 : end do
730 :
731 0 : if ( head_node ) then
732 0 : write(io_unit,12) "# || Q^t.Q - I || = ",sqrt(sum(abs(error(:,:))**2))
733 0 : flush(io_unit)
734 : end if
735 :
736 :
737 0 : ABI_FREE(error)
738 :
739 0 : ABI_MALLOC(error,(Hsize,Xsize))
740 :
741 0 : error = Xmatrix
742 :
743 0 : do l2=1,Xsize
744 0 : do l1=1,Xsize
745 0 : error(:,l2) = error(:,l2) - Qinternal(:,l1)*Rinternal(l1,l2)
746 : end do
747 : end do
748 :
749 0 : real_value = zero
750 0 : do l2=1,Xsize
751 0 : do l1=1,Xsize
752 0 : cmplx_value = complex_vector_product(error(:,l1),error(:,l2),Hsize)
753 :
754 0 : real_value = real_value + abs(cmplx_value)**2
755 : end do
756 : end do
757 :
758 0 : call xmpi_sum(real_value,mpi_communicator,ierr) ! sum on all processors
759 :
760 :
761 0 : real_value = sqrt(real_value)
762 :
763 0 : if ( head_node) then
764 0 : write(io_unit,12) "# || Xin - Q.R || = ",real_value
765 :
766 0 : if ( real_value > 1.0D-10 ) write(io_unit,10) "# ERROR! "
767 :
768 0 : write(io_unit,10) '#'
769 0 : write(io_unit,10) '# R matrix'
770 0 : write(io_unit,10) '#'
771 0 : do l1=1, Xsize
772 0 : write(io_unit,30) Rinternal(l1,:)
773 : end do
774 0 : write(io_unit,10) ''
775 0 : write(io_unit,10) ''
776 :
777 :
778 0 : write(io_unit,10) '#'
779 0 : write(io_unit,10) '# top of A matrix'
780 0 : write(io_unit,10) '#'
781 0 : do l1=1, 2*Xsize
782 0 : write(io_unit,30) A_matrix(l1,:)
783 : end do
784 0 : write(io_unit,10) ''
785 0 : write(io_unit,10) ''
786 :
787 0 : close(io_unit)
788 :
789 : end if
790 :
791 0 : ABI_FREE(error)
792 :
793 : end if
794 :
795 : !--------------------------------------------------------------------------------
796 : ! Final assignment
797 : !
798 : !--------------------------------------------------------------------------------
799 :
800 486061 : Xmatrix = Qinternal
801 :
802 2105 : if (present(Rmatrix)) then
803 13636 : Rmatrix = Rinternal
804 : end if
805 :
806 2105 : ABI_FREE(Qinternal)
807 2105 : ABI_FREE(Rinternal)
808 2105 : ABI_FREE(nproc_array)
809 2105 : ABI_FREE(A_matrix)
810 2105 : ABI_FREE(list_beta)
811 :
812 :
813 : 10 format(A)
814 : 11 format(A,I8)
815 : 12 format(A,E24.16)
816 : 13 format(A,2X,L10)
817 : 20 format(A,1000I10)
818 : 25 format(A,I5,A)
819 : 30 format(1000(F22.16,2X,F22.16,5X))
820 : 40 format(A,1000(F22.16,5X))
821 :
822 6315 : end subroutine extract_QR_Householder
823 : !!***
824 :
825 :
826 :
827 :
828 : end module m_gwls_QR_factorization
829 : !!***
|
