1!--------------------------------------------------------------------------------------------------!
2!   CP2K: A general program to perform molecular dynamics simulations                              !
3!   Copyright (C) 2000 - 2020  CP2K developers group                                               !
4!--------------------------------------------------------------------------------------------------!
5
6! **************************************************************************************************
7!> \brief various cholesky decomposition related routines
8!> \par History
9!>      09.2002 created [fawzi]
10!> \author Fawzi Mohamed
11! **************************************************************************************************
12MODULE cp_fm_cholesky
13   USE cp_blacs_env,                    ONLY: cp_blacs_env_type
14   USE cp_fm_types,                     ONLY: cp_fm_type
15   USE cp_log_handling,                 ONLY: cp_to_string
16   USE kinds,                           ONLY: dp,&
17                                              sp
18#include "../base/base_uses.f90"
19
20   IMPLICIT NONE
21   PRIVATE
22
23   LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
24   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cp_fm_cholesky'
25
26   PUBLIC :: cp_fm_cholesky_decompose, cp_fm_cholesky_invert, &
27             cp_fm_cholesky_reduce, cp_fm_cholesky_restore
28
29!***
30CONTAINS
31
32! **************************************************************************************************
33!> \brief used to replace a symmetric positive def. matrix M with its cholesky
34!>      decomposition U: M = U^T * U, with U upper triangular
35!> \param matrix the matrix to replace with its cholesky decomposition
36!> \param n the number of row (and columns) of the matrix &
37!>        (defaults to the min(size(matrix)))
38!> \param info_out ...
39!> \par History
40!>      05.2002 created [JVdV]
41!>      12.2002 updated, added n optional parm [fawzi]
42!> \author Joost
43! **************************************************************************************************
44   SUBROUTINE cp_fm_cholesky_decompose(matrix, n, info_out)
45      TYPE(cp_fm_type), POINTER                :: matrix
46      INTEGER, INTENT(in), OPTIONAL            :: n
47      INTEGER, INTENT(out), OPTIONAL           :: info_out
48
49      CHARACTER(len=*), PARAMETER :: routineN = 'cp_fm_cholesky_decompose', &
50                                     routineP = moduleN//':'//routineN
51
52      INTEGER                                  :: handle, info, my_n
53      REAL(KIND=dp), DIMENSION(:, :), POINTER  :: a
54      REAL(KIND=sp), DIMENSION(:, :), POINTER  :: a_sp
55#if defined(__SCALAPACK)
56      INTEGER, DIMENSION(9)                    :: desca
57#endif
58
59      CALL timeset(routineN, handle)
60
61      my_n = MIN(matrix%matrix_struct%nrow_global, &
62                 matrix%matrix_struct%ncol_global)
63      IF (PRESENT(n)) THEN
64         CPASSERT(n <= my_n)
65         my_n = n
66      END IF
67
68      a => matrix%local_data
69      a_sp => matrix%local_data_sp
70
71#if defined(__SCALAPACK)
72      desca(:) = matrix%matrix_struct%descriptor(:)
73
74      IF (matrix%use_sp) THEN
75         CALL pspotrf('U', my_n, a_sp(1, 1), 1, 1, desca, info)
76      ELSE
77         CALL pdpotrf('U', my_n, a(1, 1), 1, 1, desca, info)
78      ENDIF
79
80#else
81
82      IF (matrix%use_sp) THEN
83         CALL spotrf('U', my_n, a_sp(1, 1), SIZE(a_sp, 1), info)
84      ELSE
85         CALL dpotrf('U', my_n, a(1, 1), SIZE(a, 1), info)
86      ENDIF
87
88#endif
89
90      IF (PRESENT(info_out)) THEN
91         info_out = info
92      ELSE
93         IF (info /= 0) &
94            CALL cp_abort(__LOCATION__, &
95                          "Cholesky decompose failed: the matrix is not positive definite or ill-conditioned.")
96      END IF
97
98      CALL timestop(handle)
99
100   END SUBROUTINE cp_fm_cholesky_decompose
101
102! **************************************************************************************************
103!> \brief used to replace the cholesky decomposition by the inverse
104!> \param matrix the matrix to invert (must be an upper triangular matrix)
105!> \param n size of the matrix to invert (defaults to the min(size(matrix)))
106!> \param info_out ...
107!> \par History
108!>      05.2002 created [JVdV]
109!> \author Joost VandeVondele
110! **************************************************************************************************
111   SUBROUTINE cp_fm_cholesky_invert(matrix, n, info_out)
112      TYPE(cp_fm_type), POINTER                 :: matrix
113      INTEGER, INTENT(in), OPTIONAL             :: n
114      INTEGER, INTENT(OUT), OPTIONAL            :: info_out
115
116      CHARACTER(len=*), PARAMETER :: routineN = 'cp_fm_cholesky_invert', &
117                                     routineP = moduleN//':'//routineN
118      REAL(KIND=dp), DIMENSION(:, :), POINTER  :: a
119      REAL(KIND=sp), DIMENSION(:, :), POINTER  :: a_sp
120      INTEGER                                   :: info, handle
121      INTEGER                                   :: my_n
122#if defined(__SCALAPACK)
123      INTEGER, DIMENSION(9)                     :: desca
124#endif
125
126      CALL timeset(routineN, handle)
127
128      my_n = MIN(matrix%matrix_struct%nrow_global, &
129                 matrix%matrix_struct%ncol_global)
130      IF (PRESENT(n)) THEN
131         CPASSERT(n <= my_n)
132         my_n = n
133      END IF
134
135      a => matrix%local_data
136      a_sp => matrix%local_data_sp
137
138#if defined(__SCALAPACK)
139
140      desca(:) = matrix%matrix_struct%descriptor(:)
141
142      IF (matrix%use_sp) THEN
143         CALL pspotri('U', my_n, a_sp(1, 1), 1, 1, desca, info)
144      ELSE
145         CALL pdpotri('U', my_n, a(1, 1), 1, 1, desca, info)
146      ENDIF
147
148#else
149
150      IF (matrix%use_sp) THEN
151         CALL spotri('U', my_n, a_sp(1, 1), SIZE(a_sp, 1), info)
152      ELSE
153         CALL dpotri('U', my_n, a(1, 1), SIZE(a, 1), info)
154      ENDIF
155
156#endif
157
158      IF (PRESENT(info_out)) THEN
159         info_out = info
160      ELSE
161         IF (info /= 0) &
162            CPABORT("Cholesky invert failed: the matrix is not positive definite or ill-conditioned.")
163      END IF
164
165      CALL timestop(handle)
166
167   END SUBROUTINE cp_fm_cholesky_invert
168
169! **************************************************************************************************
170!> \brief reduce a matrix pencil A,B to normal form
171!>      B has to be cholesky decomposed with  cp_fm_cholesky_decompose
172!>      before calling this routine
173!>      A,B -> inv(U^T)*A*inv(U),1
174!>      (AX=BX -> inv(U^T)*A*inv(U)*U*X=U*X hence evecs U*X)
175!> \param matrix the symmetric matrix A
176!> \param matrixb the cholesky decomposition of matrix B
177!> \param itype ...
178!> \par History
179!>      05.2002 created [JVdV]
180!> \author Joost VandeVondele
181! **************************************************************************************************
182   SUBROUTINE cp_fm_cholesky_reduce(matrix, matrixb, itype)
183      TYPE(cp_fm_type), POINTER           :: matrix, matrixb
184      INTEGER, OPTIONAL                   :: itype
185
186      CHARACTER(len=*), PARAMETER :: routineN = 'cp_fm_cholesky_reduce', &
187                                     routineP = moduleN//':'//routineN
188      REAL(KIND=dp), DIMENSION(:, :), POINTER  :: a, b
189      INTEGER                                   :: info, handle
190      INTEGER                                   :: n, my_itype
191#if defined(__SCALAPACK)
192      REAL(KIND=dp)                             :: scale
193      INTEGER, DIMENSION(9)                     :: desca, descb
194#endif
195
196      CALL timeset(routineN, handle)
197
198      n = matrix%matrix_struct%nrow_global
199
200      my_itype = 1
201      IF (PRESENT(itype)) my_itype = itype
202
203      a => matrix%local_data
204      b => matrixb%local_data
205
206#if defined(__SCALAPACK)
207
208      desca(:) = matrix%matrix_struct%descriptor(:)
209      descb(:) = matrixb%matrix_struct%descriptor(:)
210
211      CALL pdsygst(my_itype, 'U', n, a(1, 1), 1, 1, desca, b(1, 1), 1, 1, descb, scale, info)
212
213      ! this is supposed to be one in current version of lapack
214      ! if not, eigenvalues have to be scaled by this number
215      IF (scale /= 1.0_dp) &
216         CPABORT("scale not equal 1 (scale="//cp_to_string(scale)//")")
217#else
218
219      CALL dsygst(my_itype, 'U', n, a(1, 1), n, b(1, 1), n, info)
220
221#endif
222
223      CPASSERT(info == 0)
224
225      CALL timestop(handle)
226
227   END SUBROUTINE cp_fm_cholesky_reduce
228
229!
230! op can be "SOLVE" (out = U^-1 * in ) or "MULTIPLY"   (out = U * in )
231! pos can be "LEFT" or "RIGHT" (U at the left or at the right)
232!
233! DEPRECATED, see cp_fm_basic_linalg:cp_fm_triangular_multiply
234!
235! **************************************************************************************************
236!> \brief ...
237!> \param matrix ...
238!> \param neig ...
239!> \param matrixb ...
240!> \param matrixout ...
241!> \param op ...
242!> \param pos ...
243!> \param transa ...
244! **************************************************************************************************
245   SUBROUTINE cp_fm_cholesky_restore(matrix, neig, matrixb, matrixout, op, pos, transa)
246      TYPE(cp_fm_type), POINTER          :: matrix, matrixb, matrixout
247      INTEGER, INTENT(IN)                         :: neig
248      CHARACTER(LEN=*), INTENT(IN)        :: op
249      CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: pos
250      CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: transa
251
252      CHARACTER(len=*), PARAMETER :: routineN = 'cp_fm_cholesky_restore', &
253                                     routineP = moduleN//':'//routineN
254      REAL(KIND=dp), DIMENSION(:, :), POINTER         :: a, b, out
255      REAL(KIND=sp), DIMENSION(:, :), POINTER         :: a_sp, b_sp, out_sp
256      INTEGER                                   :: itype, handle
257      INTEGER                                   :: n
258      REAL(KIND=dp)                           :: alpha
259      INTEGER                                   :: myprow, mypcol
260      TYPE(cp_blacs_env_type), POINTER          :: context
261      CHARACTER                                 :: chol_pos, chol_transa
262#if defined(__SCALAPACK)
263      INTEGER                                   :: i
264      INTEGER, DIMENSION(9)                     :: desca, descb, descout
265#endif
266
267      CALL timeset(routineN, handle)
268
269      context => matrix%matrix_struct%context
270      myprow = context%mepos(1)
271      mypcol = context%mepos(2)
272      n = matrix%matrix_struct%nrow_global
273      itype = 1
274      IF (op /= "SOLVE" .AND. op /= "MULTIPLY") &
275         CPABORT("wrong argument op")
276
277      IF (PRESENT(pos)) THEN
278         SELECT CASE (pos)
279         CASE ("LEFT")
280            chol_pos = 'L'
281         CASE ("RIGHT")
282            chol_pos = 'R'
283         CASE DEFAULT
284            CPABORT("wrong argument pos")
285         END SELECT
286      ELSE
287         chol_pos = 'L'
288      ENDIF
289
290      chol_transa = 'N'
291      IF (PRESENT(transa)) chol_transa = transa
292
293      IF ((matrix%use_sp .NEQV. matrixb%use_sp) .OR. (matrix%use_sp .NEQV. matrixout%use_sp)) &
294         CPABORT("not the same precision")
295
296      ! notice b is the cholesky guy
297      a => matrix%local_data
298      b => matrixb%local_data
299      out => matrixout%local_data
300      a_sp => matrix%local_data_sp
301      b_sp => matrixb%local_data_sp
302      out_sp => matrixout%local_data_sp
303
304#if defined(__SCALAPACK)
305
306      desca(:) = matrix%matrix_struct%descriptor(:)
307      descb(:) = matrixb%matrix_struct%descriptor(:)
308      descout(:) = matrixout%matrix_struct%descriptor(:)
309      alpha = 1.0_dp
310      DO i = 1, neig
311         IF (matrix%use_sp) THEN
312            CALL pscopy(n, a_sp(1, 1), 1, i, desca, 1, out_sp(1, 1), 1, i, descout, 1)
313         ELSE
314            CALL pdcopy(n, a(1, 1), 1, i, desca, 1, out(1, 1), 1, i, descout, 1)
315         ENDIF
316      ENDDO
317      IF (op .EQ. "SOLVE") THEN
318         IF (matrix%use_sp) THEN
319            CALL pstrsm(chol_pos, 'U', chol_transa, 'N', n, neig, REAL(alpha, sp), b_sp(1, 1), 1, 1, descb, &
320                        out_sp(1, 1), 1, 1, descout)
321         ELSE
322            CALL pdtrsm(chol_pos, 'U', chol_transa, 'N', n, neig, alpha, b(1, 1), 1, 1, descb, out(1, 1), 1, 1, descout)
323         ENDIF
324      ELSE
325         IF (matrix%use_sp) THEN
326            CALL pstrmm(chol_pos, 'U', chol_transa, 'N', n, neig, REAL(alpha, sp), b_sp(1, 1), 1, 1, descb, &
327                        out_sp(1, 1), 1, 1, descout)
328         ELSE
329            CALL pdtrmm(chol_pos, 'U', chol_transa, 'N', n, neig, alpha, b(1, 1), 1, 1, descb, out(1, 1), 1, 1, descout)
330         ENDIF
331      ENDIF
332#else
333
334      alpha = 1.0_dp
335      IF (matrix%use_sp) THEN
336         CALL scopy(neig*n, a_sp(1, 1), 1, out_sp(1, 1), 1)
337      ELSE
338         CALL dcopy(neig*n, a(1, 1), 1, out(1, 1), 1)
339      ENDIF
340      IF (op .EQ. "SOLVE") THEN
341         IF (matrix%use_sp) THEN
342            CALL strsm(chol_pos, 'U', chol_transa, 'N', n, neig, REAL(alpha, sp), b_sp(1, 1), SIZE(b_sp, 1), out_sp(1, 1), n)
343         ELSE
344            CALL dtrsm(chol_pos, 'U', chol_transa, 'N', n, neig, alpha, b(1, 1), SIZE(b, 1), out(1, 1), n)
345         ENDIF
346      ELSE
347         IF (matrix%use_sp) THEN
348            CALL strmm(chol_pos, 'U', chol_transa, 'N', n, neig, REAL(alpha, sp), b_sp(1, 1), n, out_sp(1, 1), n)
349         ELSE
350            CALL dtrmm(chol_pos, 'U', chol_transa, 'N', n, neig, alpha, b(1, 1), n, out(1, 1), n)
351         ENDIF
352      ENDIF
353
354#endif
355
356      CALL timestop(handle)
357
358   END SUBROUTINE cp_fm_cholesky_restore
359
360END MODULE cp_fm_cholesky
361