23 #include "../base/base_uses.f90"
28 LOGICAL,
PRIVATE,
PARAMETER :: debug_this_module = .true.
29 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'cp_fm_cholesky'
50 TYPE(cp_fm_type),
INTENT(IN) :: matrix
51 INTEGER,
INTENT(in),
OPTIONAL :: n
52 INTEGER,
INTENT(out),
OPTIONAL :: info_out
54 CHARACTER(len=*),
PARAMETER :: routinen =
'cp_fm_cholesky_decompose'
56 INTEGER :: handle, info, my_n
57 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: a
58 REAL(kind=
sp),
DIMENSION(:, :),
POINTER :: a_sp
59 #if defined(__SCALAPACK)
60 INTEGER,
DIMENSION(9) :: desca
63 CALL timeset(routinen, handle)
65 my_n = min(matrix%matrix_struct%nrow_global, &
66 matrix%matrix_struct%ncol_global)
72 a => matrix%local_data
73 a_sp => matrix%local_data_sp
75 #if defined(__SCALAPACK)
76 desca(:) = matrix%matrix_struct%descriptor(:)
78 IF (matrix%use_sp)
THEN
84 IF (matrix%use_sp)
THEN
85 CALL pspotrf(
'U', my_n, a_sp(1, 1), 1, 1, desca, info)
87 CALL pdpotrf(
'U', my_n, a(1, 1), 1, 1, desca, info)
92 IF (matrix%use_sp)
THEN
93 CALL spotrf(
'U', my_n, a_sp(1, 1),
SIZE(a_sp, 1), info)
95 CALL dpotrf(
'U', my_n, a(1, 1),
SIZE(a, 1), info)
100 IF (
PRESENT(info_out))
THEN
104 CALL cp_abort(__location__, &
105 "Cholesky decompose failed: the matrix is not positive definite or ill-conditioned.")
108 CALL timestop(handle)
122 TYPE(cp_fm_type),
INTENT(IN) :: matrix
123 INTEGER,
INTENT(in),
OPTIONAL :: n
124 INTEGER,
INTENT(OUT),
OPTIONAL :: info_out
126 CHARACTER(len=*),
PARAMETER :: routinen =
'cp_fm_cholesky_invert'
127 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: a
128 REAL(kind=
sp),
DIMENSION(:, :),
POINTER :: a_sp
129 INTEGER :: info, handle
131 #if defined(__SCALAPACK)
132 INTEGER,
DIMENSION(9) :: desca
135 CALL timeset(routinen, handle)
137 my_n = min(matrix%matrix_struct%nrow_global, &
138 matrix%matrix_struct%ncol_global)
144 a => matrix%local_data
145 a_sp => matrix%local_data_sp
147 #if defined(__SCALAPACK)
149 desca(:) = matrix%matrix_struct%descriptor(:)
151 IF (matrix%use_sp)
THEN
152 CALL pspotri(
'U', my_n, a_sp(1, 1), 1, 1, desca, info)
154 CALL pdpotri(
'U', my_n, a(1, 1), 1, 1, desca, info)
159 IF (matrix%use_sp)
THEN
160 CALL spotri(
'U', my_n, a_sp(1, 1),
SIZE(a_sp, 1), info)
162 CALL dpotri(
'U', my_n, a(1, 1),
SIZE(a, 1), info)
167 IF (
PRESENT(info_out))
THEN
171 cpabort(
"Cholesky invert failed: the matrix is not positive definite or ill-conditioned.")
174 CALL timestop(handle)
192 TYPE(cp_fm_type),
INTENT(IN) :: matrix, matrixb
193 INTEGER,
OPTIONAL :: itype
195 CHARACTER(len=*),
PARAMETER :: routinen =
'cp_fm_cholesky_reduce'
196 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: a, b
197 INTEGER :: info, handle
198 INTEGER :: n, my_itype
199 #if defined(__SCALAPACK)
200 REAL(kind=
dp) :: scale
201 INTEGER,
DIMENSION(9) :: desca, descb
204 CALL timeset(routinen, handle)
206 n = matrix%matrix_struct%nrow_global
209 IF (
PRESENT(itype)) my_itype = itype
211 a => matrix%local_data
212 b => matrixb%local_data
214 #if defined(__SCALAPACK)
216 desca(:) = matrix%matrix_struct%descriptor(:)
217 descb(:) = matrixb%matrix_struct%descriptor(:)
219 CALL pdsygst(my_itype,
'U', n, a(1, 1), 1, 1, desca, b(1, 1), 1, 1, descb, scale, info)
223 IF (scale /= 1.0_dp) &
224 cpabort(
"scale not equal 1 (scale="//cp_to_string(scale)//
")")
227 CALL dsygst(my_itype,
'U', n, a(1, 1), n, b(1, 1), n, info)
233 CALL timestop(handle)
254 TYPE(cp_fm_type),
INTENT(IN) :: matrix, matrixb, matrixout
255 INTEGER,
INTENT(IN) :: neig
256 CHARACTER(LEN=*),
INTENT(IN) :: op
257 CHARACTER(LEN=*),
INTENT(IN),
OPTIONAL :: pos
258 CHARACTER(LEN=*),
INTENT(IN),
OPTIONAL :: transa
260 CHARACTER(len=*),
PARAMETER :: routinen =
'cp_fm_cholesky_restore'
261 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: a, b, out
262 REAL(kind=
sp),
DIMENSION(:, :),
POINTER :: a_sp, b_sp, out_sp
263 INTEGER :: itype, handle
265 REAL(kind=
dp) :: alpha
266 INTEGER :: myprow, mypcol
267 TYPE(cp_blacs_env_type),
POINTER :: context
268 CHARACTER :: chol_pos, chol_transa
269 #if defined(__SCALAPACK)
271 INTEGER,
DIMENSION(9) :: desca, descb, descout
274 CALL timeset(routinen, handle)
276 context => matrix%matrix_struct%context
277 myprow = context%mepos(1)
278 mypcol = context%mepos(2)
279 n = matrix%matrix_struct%nrow_global
281 IF (op /=
"SOLVE" .AND. op /=
"MULTIPLY") &
282 cpabort(
"wrong argument op")
284 IF (
PRESENT(pos))
THEN
291 cpabort(
"wrong argument pos")
298 IF (
PRESENT(transa)) chol_transa = transa
300 IF ((matrix%use_sp .NEQV. matrixb%use_sp) .OR. (matrix%use_sp .NEQV. matrixout%use_sp)) &
301 cpabort(
"not the same precision")
304 a => matrix%local_data
305 b => matrixb%local_data
306 out => matrixout%local_data
307 a_sp => matrix%local_data_sp
308 b_sp => matrixb%local_data_sp
309 out_sp => matrixout%local_data_sp
311 #if defined(__SCALAPACK)
313 desca(:) = matrix%matrix_struct%descriptor(:)
314 descb(:) = matrixb%matrix_struct%descriptor(:)
315 descout(:) = matrixout%matrix_struct%descriptor(:)
318 IF (matrix%use_sp)
THEN
319 CALL pscopy(n, a_sp(1, 1), 1, i, desca, 1, out_sp(1, 1), 1, i, descout, 1)
321 CALL pdcopy(n, a(1, 1), 1, i, desca, 1, out(1, 1), 1, i, descout, 1)
324 IF (op .EQ.
"SOLVE")
THEN
325 IF (matrix%use_sp)
THEN
326 CALL pstrsm(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)
329 CALL pdtrsm(chol_pos,
'U', chol_transa,
'N', n, neig, alpha, b(1, 1), 1, 1, descb, out(1, 1), 1, 1, descout)
332 IF (matrix%use_sp)
THEN
333 CALL pstrmm(chol_pos,
'U', chol_transa,
'N', n, neig, real(alpha,
sp), b_sp(1, 1), 1, 1, descb, &
334 out_sp(1, 1), 1, 1, descout)
336 CALL pdtrmm(chol_pos,
'U', chol_transa,
'N', n, neig, alpha, b(1, 1), 1, 1, descb, out(1, 1), 1, 1, descout)
342 IF (matrix%use_sp)
THEN
343 CALL scopy(neig*n, a_sp(1, 1), 1, out_sp(1, 1), 1)
345 CALL dcopy(neig*n, a(1, 1), 1, out(1, 1), 1)
347 IF (op .EQ.
"SOLVE")
THEN
348 IF (matrix%use_sp)
THEN
349 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)
351 CALL dtrsm(chol_pos,
'U', chol_transa,
'N', n, neig, alpha, b(1, 1),
SIZE(b, 1), out(1, 1), n)
354 IF (matrix%use_sp)
THEN
355 CALL strmm(chol_pos,
'U', chol_transa,
'N', n, neig, real(alpha,
sp), b_sp(1, 1), n, out_sp(1, 1), n)
357 CALL dtrmm(chol_pos,
'U', chol_transa,
'N', n, neig, alpha, b(1, 1), n, out(1, 1), n)
363 CALL timestop(handle)
methods related to the blacs parallel environment
various cholesky decomposition related routines
subroutine, public cp_fm_cholesky_invert(matrix, n, info_out)
used to replace the cholesky decomposition by the inverse
subroutine, public cp_fm_cholesky_decompose(matrix, n, info_out)
used to replace a symmetric positive def. matrix M with its cholesky decomposition U: M = U^T * U,...
subroutine, public cp_fm_cholesky_reduce(matrix, matrixb, itype)
reduce a matrix pencil A,B to normal form B has to be cholesky decomposed with cp_fm_cholesky_decompo...
subroutine, public cp_fm_cholesky_restore(matrix, neig, matrixb, matrixout, op, pos, transa)
...
subroutine, public cp_pdpotrf_dlaf(uplo, n, a, ia, ja, desca, info)
Cholesky factorization using DLA-Future.
subroutine, public cp_pspotrf_dlaf(uplo, n, a, ia, ja, desca, info)
Cholesky factorization using DLA-Future.
represent a full matrix distributed on many processors
various routines to log and control the output. The idea is that decisions about where to log should ...
Defines the basic variable types.
integer, parameter, public dp
integer, parameter, public sp