Branch data 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 : : !!***
|
