|
CP2K 2.4 (Revision 12889)
|
00001 !-----------------------------------------------------------------------------! 00002 ! CP2K: A general program to perform molecular dynamics simulations ! 00003 ! Copyright (C) 2000 - 2013 CP2K developers group ! 00004 !-----------------------------------------------------------------------------! 00005 00006 ! ***************************************************************************** 00012 MODULE cp_fm_basic_linalg 00013 USE cp_fm_struct, ONLY: cp_fm_struct_equivalent 00014 USE cp_fm_types, ONLY: & 00015 cp_fm_create, cp_fm_get_element, cp_fm_get_info, cp_fm_get_submatrix, & 00016 cp_fm_release, cp_fm_set_all, cp_fm_set_submatrix, cp_fm_to_fm, & 00017 cp_fm_type 00018 USE cp_para_types, ONLY: cp_blacs_env_type 00019 USE kahan_sum, ONLY: accurate_sum 00020 USE kinds, ONLY: dp,& 00021 sp 00022 USE mathlib, ONLY: invert_matrix 00023 USE message_passing, ONLY: mp_sum 00024 USE termination, ONLY: stop_program 00025 USE timings, ONLY: timeset,& 00026 timestop 00027 #include "cp_common_uses.h" 00028 00029 IMPLICIT NONE 00030 PRIVATE 00031 00032 LOGICAL, PRIVATE, PARAMETER :: debug_this_module=.TRUE. 00033 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cp_fm_basic_linalg' 00034 00035 PUBLIC :: cp_fm_scale, & ! scale a matrix 00036 cp_fm_scale_and_add, & ! scale and add two matrices 00037 cp_fm_column_scale, & ! scale colummns of a matrix 00038 cp_fm_trace, & ! trace of the transpose(A)*B 00039 cp_fm_schur_product, & ! schur product 00040 cp_fm_transpose, & ! transpose a matrix 00041 cp_fm_upper_to_full, & ! symmetrise an upper symmetric matrix 00042 cp_fm_syrk, & ! rank k update 00043 cp_fm_triangular_multiply, & ! triangular matrix multiply / solve 00044 cp_fm_symm, & ! multiply a symmetric with a non-symmetric matrix 00045 cp_fm_gemm, & ! multiply two matrices 00046 cp_fm_lu_decompose,& ! computes determinant (and lu decomp) 00047 cp_fm_invert,& ! computes the inverse and determinant 00048 cp_fm_frobenius_norm,& ! frobenius norm 00049 cp_fm_ger, & ! rank 1 operation 00050 cp_fm_triangular_invert,& ! compute the reciprocal of a tirangular matrix 00051 cp_fm_qr_factorization,& ! compute the QR factorization of a rectangular matrix 00052 cp_fm_solve ! solves the equation A*B=C A and C are input 00053 00054 CONTAINS 00055 00056 ! ***************************************************************************** 00061 SUBROUTINE cp_fm_scale_and_add(alpha,matrix_a,beta,matrix_b,error) 00062 00063 REAL(KIND=dp), INTENT(IN) :: alpha 00064 TYPE(cp_fm_type), POINTER :: matrix_a 00065 REAL(KIND=dp), INTENT(in), OPTIONAL :: beta 00066 TYPE(cp_fm_type), OPTIONAL, POINTER :: matrix_b 00067 TYPE(cp_error_type), INTENT(inout) :: error 00068 00069 CHARACTER(len=*), PARAMETER :: routineN = 'cp_fm_scale_and_add', 00070 routineP = moduleN//':'//routineN 00071 00072 INTEGER :: handle, size_a 00073 LOGICAL :: failure 00074 REAL(KIND=dp) :: my_beta 00075 REAL(KIND=dp), DIMENSION(:, :), POINTER :: a, b 00076 REAL(KIND=sp), DIMENSION(:, :), POINTER :: a_sp, b_sp 00077 00078 CALL timeset(routineN,handle) 00079 00080 failure=.FALSE. 00081 IF (PRESENT(matrix_b)) THEN 00082 my_beta=1.0_dp 00083 ELSE 00084 my_beta=0.0_dp 00085 ENDIF 00086 IF(PRESENT(beta)) my_beta=beta 00087 NULLIFY(a,b) 00088 00089 CPPrecondition(ASSOCIATED(matrix_a),cp_failure_level,routineP,error,failure) 00090 CPPrecondition(matrix_a%ref_count>0,cp_failure_level,routineP,error,failure) 00091 00092 IF (PRESENT(beta)) THEN 00093 CPPrecondition(PRESENT(matrix_b),cp_failure_level,routineP,error,failure) 00094 CPPrecondition(ASSOCIATED(matrix_b),cp_failure_level,routineP,error,failure) 00095 CPPrecondition(matrix_b%ref_count>0,cp_failure_level,routineP,error,failure) 00096 IF (matrix_a%id_nr==matrix_b%id_nr) THEN 00097 CALL cp_assert(matrix_a%id_nr/=matrix_b%id_nr, & 00098 cp_warning_level, cp_assertion_failed, & 00099 fromWhere=routineP, & 00100 message="Bad use of routine. Call cp_fm_scale instead: "// & 00101 CPSourceFileRef, & 00102 error=error, failure=failure) 00103 CALL cp_fm_scale(alpha+beta, matrix_a, error=error) 00104 CALL timestop(handle) 00105 RETURN 00106 END IF 00107 END IF 00108 00109 a => matrix_a%local_data 00110 a_sp => matrix_a%local_data_sp 00111 00112 IF(matrix_a%use_sp) THEN 00113 size_a = SIZE(a_sp,1)*SIZE(a_sp,2) 00114 ELSE 00115 size_a = SIZE(a,1)*SIZE(a,2) 00116 ENDIF 00117 00118 IF (alpha /= 1.0_dp) THEN 00119 IF(matrix_a%use_sp) THEN 00120 CALL sscal ( size_a, REAL(alpha,sp), a_sp, 1) 00121 ELSE 00122 CALL dscal ( size_a, alpha, a, 1) 00123 ENDIF 00124 ENDIF 00125 IF (my_beta.NE.0.0_dp) THEN 00126 CALL cp_assert(matrix_a%matrix_struct%context%group==& 00127 matrix_b%matrix_struct%context%group,cp_failure_level,& 00128 cp_assertion_failed,fromWhere=routineP,& 00129 message="matrixes must be in the same blacs context"//& 00130 CPSourceFileRef,& 00131 error=error,failure=failure) 00132 00133 IF (cp_fm_struct_equivalent(matrix_a%matrix_struct,& 00134 matrix_b%matrix_struct,error=error)) THEN 00135 00136 b => matrix_b%local_data 00137 b_sp => matrix_b%local_data_sp 00138 00139 IF(matrix_a%use_sp.AND.matrix_b%use_sp) THEN 00140 CALL saxpy ( size_a, REAL(my_beta,sp), b_sp, 1, a_sp, 1 ) 00141 ELSEIF(matrix_a%use_sp.AND..NOT.matrix_b%use_sp) THEN 00142 CALL saxpy ( size_a, REAL(my_beta,sp), REAL(b,sp), 1, a_sp, 1 ) 00143 ELSEIF(.NOT.matrix_a%use_sp.AND.matrix_b%use_sp) THEN 00144 CALL daxpy ( size_a, my_beta, REAL(b_sp,dp), 1, a, 1 ) 00145 ELSE 00146 CALL daxpy ( size_a, my_beta, b, 1, a, 1 ) 00147 ENDIF 00148 00149 ELSE 00150 #ifdef __SCALAPACK 00151 CALL cp_unimplemented_error(fromWhere=routineP, & 00152 message="to do (pdscal,pdcopy,pdaxpy)", error=error) 00153 #else 00154 CPAssert(.FALSE.,cp_failure_level,routineP,error,failure) 00155 #endif 00156 END IF 00157 00158 END IF 00159 00160 CALL timestop(handle) 00161 00162 END SUBROUTINE cp_fm_scale_and_add 00163 00164 ! ***************************************************************************** 00177 SUBROUTINE cp_fm_lu_decompose(matrix_a,almost_determinant,correct_sign) 00178 TYPE(cp_fm_type), POINTER :: matrix_a 00179 REAL(KIND=dp), INTENT(OUT) :: almost_determinant 00180 LOGICAL, INTENT(IN), OPTIONAL :: correct_sign 00181 00182 CHARACTER(len=*), PARAMETER :: routineN = 'cp_fm_lu_decompose', 00183 routineP = moduleN//':'//routineN 00184 00185 INTEGER :: handle, i, info, lda, n 00186 INTEGER, ALLOCATABLE, DIMENSION(:) :: ipivot 00187 INTEGER, DIMENSION(9) :: desca 00188 REAL(KIND=dp) :: determinant 00189 REAL(KIND=dp), DIMENSION(:), POINTER :: diag 00190 REAL(KIND=dp), DIMENSION(:, :), POINTER :: a 00191 00192 ! *** locals *** 00193 00194 CALL timeset(routineN,handle) 00195 00196 a => matrix_a%local_data 00197 n = matrix_a%matrix_struct%nrow_global 00198 ALLOCATE(ipivot(n+matrix_a%matrix_struct%nrow_block)) 00199 00200 #if defined(__SCALAPACK) 00201 desca(:) = matrix_a%matrix_struct%descriptor(:) 00202 CALL pdgetrf(n,n,a(1,1),1,1,desca,ipivot,info) 00203 00204 00205 ALLOCATE(diag(n)) 00206 diag(:)=0.0_dp 00207 DO i=1,n 00208 CALL cp_fm_get_element(matrix_a,i,i,diag(i)) ! not completely optimal in speed i would say 00209 ENDDO 00210 determinant=1.0_dp 00211 DO i=1,n 00212 determinant=determinant*diag(i) 00213 ENDDO 00214 DEALLOCATE(diag) 00215 #else 00216 lda=SIZE(a,1) 00217 CALL dgetrf(n,n,a(1,1),lda,ipivot,info) 00218 determinant=1.0_dp 00219 IF(correct_sign)THEN 00220 DO i=1,n 00221 IF(ipivot(i).NE.i)THEN 00222 determinant=-determinant*a(i,i) 00223 ELSE 00224 determinant=determinant*a(i,i) 00225 END IF 00226 END DO 00227 ELSE 00228 DO i=1,n 00229 determinant=determinant*a(i,i) 00230 ENDDO 00231 END IF 00232 #endif 00233 ! info is allowed to be zero 00234 ! this does just signal a zero diagonal element 00235 DEALLOCATE(ipivot) 00236 almost_determinant=determinant ! notice that the sign is random 00237 CALL timestop(handle) 00238 END SUBROUTINE 00239 00240 ! ***************************************************************************** 00252 SUBROUTINE cp_fm_gemm(transa,transb,m,n,k,alpha,matrix_a,matrix_b,beta,& 00253 matrix_c,error,b_first_col,a_first_row,b_first_row,c_first_col,c_first_row) 00254 00255 CHARACTER(LEN=1), INTENT(IN) :: transa, transb 00256 INTEGER, INTENT(IN) :: m, n, k 00257 REAL(KIND=dp), INTENT(IN) :: alpha 00258 TYPE(cp_fm_type), POINTER :: matrix_a, matrix_b 00259 REAL(KIND=dp), INTENT(IN) :: beta 00260 TYPE(cp_fm_type), POINTER :: matrix_c 00261 TYPE(cp_error_type), INTENT(inout) :: error 00262 INTEGER, INTENT(IN), OPTIONAL :: b_first_col, a_first_row, 00263 b_first_row, c_first_col, 00264 c_first_row 00265 00266 CHARACTER(len=*), PARAMETER :: routineN = 'cp_fm_gemm', 00267 routineP = moduleN//':'//routineN 00268 00269 INTEGER :: handle, i_a, i_b, i_c, j_b, 00270 j_c, lda, ldb, ldc 00271 INTEGER, DIMENSION(9) :: desca, descb, descc 00272 REAL(KIND=dp), DIMENSION(:, :), POINTER :: a, b, c 00273 REAL(KIND=sp), DIMENSION(:, :), POINTER :: a_sp, b_sp, c_sp 00274 00275 CALL timeset(routineN,handle) 00276 00277 a => matrix_a%local_data 00278 b => matrix_b%local_data 00279 c => matrix_c%local_data 00280 00281 a_sp => matrix_a%local_data_sp 00282 b_sp => matrix_b%local_data_sp 00283 c_sp => matrix_c%local_data_sp 00284 00285 IF (PRESENT(a_first_row)) THEN 00286 i_a = a_first_row 00287 ELSE 00288 i_a = 1 00289 END IF 00290 IF (PRESENT(b_first_row)) THEN 00291 i_b = b_first_row 00292 ELSE 00293 i_b = 1 00294 END IF 00295 IF (PRESENT(b_first_col)) THEN 00296 j_b = b_first_col 00297 ELSE 00298 j_b = 1 00299 END IF 00300 IF (PRESENT(c_first_row)) THEN 00301 i_c = c_first_row 00302 ELSE 00303 i_c = 1 00304 END IF 00305 00306 IF (PRESENT(c_first_col)) THEN 00307 j_c = c_first_col 00308 ELSE 00309 j_c = 1 00310 END IF 00311 00312 #if defined(__SCALAPACK) 00313 00314 desca(:) = matrix_a%matrix_struct%descriptor(:) 00315 descb(:) = matrix_b%matrix_struct%descriptor(:) 00316 descc(:) = matrix_c%matrix_struct%descriptor(:) 00317 00318 IF(matrix_a%use_sp.AND.matrix_b%use_sp.AND.matrix_c%use_sp) THEN 00319 00320 CALL psgemm(transa,transb,m,n,k,REAL(alpha,sp),a_sp(1,1),i_a,1,desca,b_sp(1,1),i_b,j_b, 00321 descb,REAL(beta,sp),c_sp(1,1),i_c,j_c,descc) 00322 00323 ELSEIF((.NOT.matrix_a%use_sp).AND.(.NOT.matrix_b%use_sp).AND.(.NOT.matrix_c%use_sp)) THEN 00324 00325 CALL pdgemm(transa,transb,m,n,k,alpha,a(1,1),i_a,1,desca,b(1,1),i_b,j_b,& 00326 descb,beta,c(1,1),i_c,j_c,descc) 00327 00328 ELSE 00329 CALL stop_program(routineN,moduleN,__LINE__,"Mixed precision gemm NYI") 00330 ENDIF 00331 #else 00332 00333 IF(matrix_a%use_sp.AND.matrix_b%use_sp.AND.matrix_c%use_sp) THEN 00334 00335 lda = SIZE(a_sp,1) 00336 ldb = SIZE(b_sp,1) 00337 ldc = SIZE(c_sp,1) 00338 00339 CALL sgemm(transa,transb,m,n,k,REAL(alpha,sp),a_sp(i_a,1),lda,b_sp(i_b,j_b),ldb, 00340 REAL(beta,sp),c_sp(i_c,j_c),ldc) 00341 00342 ELSEIF((.NOT.matrix_a%use_sp).AND.(.NOT.matrix_b%use_sp).AND.(.NOT.matrix_c%use_sp)) THEN 00343 00344 lda = SIZE(a,1) 00345 ldb = SIZE(b,1) 00346 ldc = SIZE(c,1) 00347 00348 CALL dgemm(transa,transb,m,n,k,alpha,a(i_a,1),lda,b(i_b,j_b),ldb,beta,c(i_c,j_c),ldc) 00349 00350 ELSE 00351 CALL stop_program(routineN,moduleN,__LINE__,"Mixed precision gemm NYI") 00352 ENDIF 00353 00354 #endif 00355 CALL timestop(handle) 00356 00357 END SUBROUTINE cp_fm_gemm 00358 00359 ! ***************************************************************************** 00375 SUBROUTINE cp_fm_symm(side,uplo,m,n,alpha,matrix_a,matrix_b,beta,matrix_c,& 00376 error) 00377 00378 CHARACTER(LEN=1), INTENT(IN) :: side, uplo 00379 INTEGER, INTENT(IN) :: m, n 00380 REAL(KIND=dp), INTENT(IN) :: alpha 00381 TYPE(cp_fm_type), POINTER :: matrix_a, matrix_b 00382 REAL(KIND=dp), INTENT(IN) :: beta 00383 TYPE(cp_fm_type), POINTER :: matrix_c 00384 TYPE(cp_error_type), INTENT(inout) :: error 00385 00386 CHARACTER(len=*), PARAMETER :: routineN = 'cp_fm_symm', 00387 routineP = moduleN//':'//routineN 00388 00389 INTEGER :: handle, lda, ldb, ldc 00390 INTEGER, DIMENSION(9) :: desca, descb, descc 00391 REAL(KIND=dp), DIMENSION(:, :), POINTER :: a, b, c 00392 00393 ! ------------------------------------------------------------------------- 00394 00395 CALL timeset(routineN,handle) 00396 00397 a => matrix_a%local_data 00398 b => matrix_b%local_data 00399 c => matrix_c%local_data 00400 00401 #if defined(__SCALAPACK) 00402 00403 desca(:) = matrix_a%matrix_struct%descriptor(:) 00404 descb(:) = matrix_b%matrix_struct%descriptor(:) 00405 descc(:) = matrix_c%matrix_struct%descriptor(:) 00406 00407 CALL pdsymm(side,uplo,m,n,alpha,a(1,1),1,1,desca,b(1,1),1,1,descb,beta,c(1,1),1,1,descc) 00408 00409 #else 00410 00411 lda = matrix_a%matrix_struct%local_leading_dimension 00412 ldb = matrix_b%matrix_struct%local_leading_dimension 00413 ldc = matrix_c%matrix_struct%local_leading_dimension 00414 00415 CALL dsymm(side,uplo,m,n,alpha,a(1,1),lda,b(1,1),ldb,beta,c(1,1),ldc) 00416 00417 #endif 00418 CALL timestop(handle) 00419 00420 END SUBROUTINE cp_fm_symm 00421 00422 ! ***************************************************************************** 00427 SUBROUTINE cp_fm_frobenius_norm(matrix_a,norm,error) 00428 TYPE(cp_fm_type), POINTER :: matrix_a 00429 REAL(KIND=dp), INTENT(inout) :: norm 00430 TYPE(cp_error_type), INTENT(inout) :: error 00431 00432 CHARACTER(len=*), PARAMETER :: routineN = 'cp_fm_frobenius_norm', 00433 routineP = moduleN//':'//routineN 00434 00435 INTEGER :: group, handle, size_a 00436 REAL(KIND=dp), DIMENSION(:, :), POINTER :: a 00437 REAL(KIND=dp), EXTERNAL :: DDOT 00438 00439 CALL timeset(routineN,handle) 00440 norm = 0.0_dp 00441 a => matrix_a%local_data 00442 size_a = SIZE(a,1)*SIZE(a,2) 00443 norm = DDOT( size_a, a(1,1), 1, a(1,1), 1 ) 00444 #if defined(__SCALAPACK) 00445 group = matrix_a%matrix_struct%para_env%group 00446 CALL mp_sum(norm,group) 00447 #endif 00448 norm = SQRT(norm) 00449 CALL timestop(handle) 00450 END SUBROUTINE cp_fm_frobenius_norm 00451 00452 ! ***************************************************************************** 00459 SUBROUTINE cp_fm_ger(alpha,vector_x,ix,jx,vector_y,iy,jy,matrix_a,error) 00460 00461 REAL(KIND=dp), INTENT(IN) :: alpha 00462 TYPE(cp_fm_type), POINTER :: vector_x 00463 INTEGER, INTENT(IN) :: ix, jx 00464 TYPE(cp_fm_type), POINTER :: vector_y 00465 INTEGER, INTENT(IN) :: iy, jy 00466 TYPE(cp_fm_type), POINTER :: matrix_a 00467 TYPE(cp_error_type), INTENT(inout) :: error 00468 00469 CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_fm_ger', 00470 routineP = moduleN//':'//routineN 00471 00472 INTEGER :: handle, ia, incx, incy, ja, 00473 lda, m, n 00474 INTEGER, DIMENSION(9) :: desca, descx, descy 00475 REAL(KIND=dp), DIMENSION(:, :), POINTER :: a, x, y 00476 00477 ! ------------------------------------------------------------------------- 00478 00479 CALL timeset(routineN,handle) 00480 00481 m = matrix_a%matrix_struct%nrow_global 00482 n = matrix_a%matrix_struct%ncol_global 00483 00484 ia = 1 00485 ja = 1 00486 00487 incx = 1 00488 incy = 1 00489 00490 a => matrix_a%local_data 00491 x => vector_x%local_data 00492 y => vector_y%local_data 00493 00494 #if defined(__SCALAPACK) 00495 00496 desca(:) = matrix_a%matrix_struct%descriptor(:) 00497 descx(:) = vector_x%matrix_struct%descriptor(:) 00498 descy(:) = vector_y%matrix_struct%descriptor(:) 00499 00500 CALL pdger(m,n,alpha,x(1,1),ix,jx,descx,incx,y(1,1),iy,jy,descy,incy,& 00501 a(1,1),ia,ja,desca) 00502 00503 #else 00504 00505 lda = SIZE(a,1) 00506 00507 CALL dger(m,n,alpha,x(ix,jx),incx,y(iy,jy),incy,a(ia,ja),lda) 00508 00509 #endif 00510 CALL timestop(handle) 00511 00512 END SUBROUTINE cp_fm_ger 00513 00514 ! ***************************************************************************** 00525 SUBROUTINE cp_fm_syrk(uplo,trans,k,alpha,matrix_a,ia,ja,beta,matrix_c,error) 00526 CHARACTER(LEN=1), INTENT(IN) :: uplo, trans 00527 INTEGER, INTENT(IN) :: k 00528 REAL(KIND=dp), INTENT(IN) :: alpha 00529 TYPE(cp_fm_type), POINTER :: matrix_a 00530 INTEGER, INTENT(IN) :: ia, ja 00531 REAL(KIND=dp), INTENT(IN) :: beta 00532 TYPE(cp_fm_type), POINTER :: matrix_c 00533 TYPE(cp_error_type), INTENT(inout) :: error 00534 00535 CHARACTER(len=*), PARAMETER :: routineN = 'cp_fm_syrk', 00536 routineP = moduleN//':'//routineN 00537 00538 INTEGER :: handle, lda, ldc, n 00539 INTEGER, DIMENSION(9) :: desca, descc 00540 REAL(KIND=dp), DIMENSION(:, :), POINTER :: a, c 00541 00542 CALL timeset(routineN,handle) 00543 00544 n = matrix_c%matrix_struct%nrow_global 00545 00546 a => matrix_a%local_data 00547 c => matrix_c%local_data 00548 00549 #if defined(__SCALAPACK) 00550 00551 desca(:) = matrix_a%matrix_struct%descriptor(:) 00552 descc(:) = matrix_c%matrix_struct%descriptor(:) 00553 00554 CALL pdsyrk(uplo,trans,n,k,alpha,a(1,1),ia,ja,desca,beta,c(1,1),1,1,descc) 00555 00556 #else 00557 00558 lda = SIZE(a,1) 00559 ldc = SIZE(c,1) 00560 00561 CALL dsyrk(uplo,trans,n,k,alpha,a(ia,ja),lda,beta,c(1,1),ldc) 00562 00563 #endif 00564 CALL timestop(handle) 00565 00566 END SUBROUTINE cp_fm_syrk 00567 00568 ! ***************************************************************************** 00573 SUBROUTINE cp_fm_schur_product(matrix_a,matrix_b,matrix_c,error) 00574 00575 TYPE(cp_fm_type), POINTER :: matrix_a, matrix_b, matrix_c 00576 TYPE(cp_error_type), INTENT(inout) :: error 00577 00578 CHARACTER(len=*), PARAMETER :: routineN = 'cp_fm_schur_product', 00579 routineP = moduleN//':'//routineN 00580 00581 INTEGER :: handle, icol_local, 00582 irow_local, mypcol, myprow, 00583 ncol_local, npcol, nprow, 00584 nrow_local 00585 REAL(KIND=dp), DIMENSION(:, :), POINTER :: a, b, c 00586 TYPE(cp_blacs_env_type), POINTER :: context 00587 00588 CALL timeset(routineN,handle) 00589 00590 context => matrix_a%matrix_struct%context 00591 myprow=context%mepos(1) 00592 mypcol=context%mepos(2) 00593 nprow=context%num_pe(1) 00594 npcol=context%num_pe(2) 00595 00596 a => matrix_a%local_data 00597 b => matrix_b%local_data 00598 c => matrix_c%local_data 00599 00600 nrow_local = matrix_a%matrix_struct%nrow_locals(myprow) 00601 ncol_local = matrix_a%matrix_struct%ncol_locals(mypcol) 00602 00603 DO icol_local=1,ncol_local 00604 DO irow_local=1,nrow_local 00605 c(irow_local,icol_local) = a(irow_local,icol_local)*b(irow_local,icol_local) 00606 END DO 00607 END DO 00608 00609 CALL timestop(handle) 00610 00611 END SUBROUTINE cp_fm_schur_product 00612 00613 ! ***************************************************************************** 00627 SUBROUTINE cp_fm_trace(matrix_a,matrix_b,trace,error) 00628 00629 TYPE(cp_fm_type), POINTER :: matrix_a, matrix_b 00630 REAL(KIND=dp), INTENT(OUT) :: trace 00631 TYPE(cp_error_type), INTENT(inout) :: error 00632 00633 CHARACTER(len=*), PARAMETER :: routineN = 'cp_fm_trace', 00634 routineP = moduleN//':'//routineN 00635 00636 INTEGER :: group, handle, mypcol, 00637 myprow, ncol_local, npcol, 00638 nprow, nrow_local 00639 REAL(KIND=dp), DIMENSION(:, :), POINTER :: a, b 00640 REAL(KIND=sp), DIMENSION(:, :), POINTER :: a_sp, b_sp 00641 TYPE(cp_blacs_env_type), POINTER :: context 00642 00643 CALL timeset(routineN,handle) 00644 00645 context => matrix_a%matrix_struct%context 00646 myprow=context%mepos(1) 00647 mypcol=context%mepos(2) 00648 nprow=context%num_pe(1) 00649 npcol=context%num_pe(2) 00650 00651 group = matrix_a%matrix_struct%para_env%group 00652 00653 a => matrix_a%local_data 00654 b => matrix_b%local_data 00655 00656 a_sp => matrix_a%local_data_sp 00657 b_sp => matrix_b%local_data_sp 00658 00659 nrow_local = MIN(matrix_a%matrix_struct%nrow_locals(myprow),matrix_b%matrix_struct%nrow_locals(myprow)) 00660 ncol_local = MIN(matrix_a%matrix_struct%ncol_locals(mypcol),matrix_b%matrix_struct%ncol_locals(mypcol)) 00661 00662 ! cries for an accurate_dot_product 00663 IF(matrix_a%use_sp.AND.matrix_b%use_sp) THEN 00664 trace = accurate_sum(REAL(a_sp(1:nrow_local,1:ncol_local) * & & b_sp(1:nrow_local,1:ncol_local), dp) ) 00665 ELSEIF(matrix_a%use_sp.AND..NOT.matrix_b%use_sp) THEN 00666 trace = accurate_sum(REAL(a_sp(1:nrow_local,1:ncol_local), dp) * 00667 b(1:nrow_local,1:ncol_local) ) 00668 ELSEIF(.NOT.matrix_a%use_sp.AND.matrix_b%use_sp) THEN 00669 trace = accurate_sum( a(1:nrow_local,1:ncol_local) * & 00670 & REAL(b_sp(1:nrow_local,1:ncol_local), dp) ) 00671 ELSE 00672 trace = accurate_sum(a(1:nrow_local,1:ncol_local) * & 00673 b(1:nrow_local,1:ncol_local) ) 00674 ENDIF 00675 00676 CALL mp_sum(trace,group) 00677 00678 CALL timestop(handle) 00679 00680 END SUBROUTINE cp_fm_trace 00681 00682 ! ***************************************************************************** 00716 SUBROUTINE cp_fm_triangular_multiply(triangular_matrix,matrix_b,side,& 00717 transpose_tr, invert_tr, uplo_tr,unit_diag_tr, n_rows, n_cols, & 00718 alpha,error) 00719 TYPE(cp_fm_type), POINTER :: triangular_matrix, matrix_b 00720 CHARACTER, INTENT(in), OPTIONAL :: side 00721 LOGICAL, INTENT(in), OPTIONAL :: transpose_tr, invert_tr 00722 CHARACTER, INTENT(in), OPTIONAL :: uplo_tr 00723 LOGICAL, INTENT(in), OPTIONAL :: unit_diag_tr 00724 INTEGER, INTENT(in), OPTIONAL :: n_rows, n_cols 00725 REAL(KIND=dp), INTENT(in), OPTIONAL :: alpha 00726 TYPE(cp_error_type), INTENT(inout) :: error 00727 00728 CHARACTER(len=*), PARAMETER :: routineN = 'cp_fm_triangular_multiply', 00729 routineP = moduleN//':'//routineN 00730 00731 CHARACTER :: side_char, transa, unit_diag, 00732 uplo 00733 INTEGER :: handle, m, n 00734 LOGICAL :: failure, invert 00735 REAL(KIND=dp) :: al 00736 00737 failure=.FALSE. 00738 00739 CALL timeset(routineN,handle) 00740 side_char='L' 00741 unit_diag='N' 00742 uplo='U' 00743 transa='N' 00744 invert=.FALSE. 00745 al=1.0_dp 00746 CALL cp_fm_get_info(matrix_b, nrow_global=m, ncol_global=n, error=error) 00747 IF (PRESENT(side)) side_char=side 00748 IF (PRESENT(invert_tr)) invert=invert_tr 00749 IF (PRESENT(uplo_tr)) uplo=uplo_tr 00750 IF (PRESENT(unit_diag_tr)) THEN 00751 IF (unit_diag_tr) THEN 00752 unit_diag='U' 00753 ELSE 00754 unit_diag='N' 00755 END IF 00756 END IF 00757 IF (PRESENT(transpose_tr)) THEN 00758 IF (transpose_tr) THEN 00759 transa='T' 00760 ELSE 00761 transa='N' 00762 END IF 00763 END IF 00764 IF (PRESENT(alpha)) al=alpha 00765 IF (PRESENT(n_rows)) m=n_rows 00766 IF (PRESENT(n_cols)) n=n_cols 00767 00768 IF (invert) THEN 00769 00770 #if defined(__SCALAPACK) 00771 CALL pdtrsm(side_char,uplo,transa,unit_diag,m,n,al,& 00772 triangular_matrix%local_data(1,1),1,1,& 00773 triangular_matrix%matrix_struct%descriptor,& 00774 matrix_b%local_data(1,1),1,1,& 00775 matrix_b%matrix_struct%descriptor(1)) 00776 #else 00777 CALL dtrsm(side_char,uplo,transa,unit_diag,m,n,al,& 00778 triangular_matrix%local_data(1,1),& 00779 SIZE(triangular_matrix%local_data,1),& 00780 matrix_b%local_data(1,1),SIZE(matrix_b%local_data,1)) 00781 #endif 00782 00783 ELSE 00784 00785 #if defined(__SCALAPACK) 00786 CALL pdtrmm(side_char,uplo,transa,unit_diag,m,n,al,& 00787 triangular_matrix%local_data(1,1),1,1,& 00788 triangular_matrix%matrix_struct%descriptor,& 00789 matrix_b%local_data(1,1),1,1,& 00790 matrix_b%matrix_struct%descriptor(1)) 00791 #else 00792 CALL dtrmm(side_char,uplo,transa,unit_diag,m,n,al,& 00793 triangular_matrix%local_data(1,1),& 00794 SIZE(triangular_matrix%local_data,1),& 00795 matrix_b%local_data(1,1),SIZE(matrix_b%local_data,1)) 00796 #endif 00797 00798 END IF 00799 00800 CALL timestop(handle) 00801 END SUBROUTINE cp_fm_triangular_multiply 00802 00803 ! ***************************************************************************** 00809 SUBROUTINE cp_fm_scale(alpha, matrix_a, error) 00810 REAL(KIND=dp), INTENT(IN) :: alpha 00811 TYPE(cp_fm_type), POINTER :: matrix_a 00812 TYPE(cp_error_type), INTENT(inout) :: error 00813 00814 CHARACTER(len=*), PARAMETER :: routineN = 'cp_fm_scale', 00815 routineP = moduleN//':'//routineN 00816 00817 INTEGER :: handle, size_a 00818 LOGICAL :: failure 00819 REAL(KIND=dp), DIMENSION(:, :), POINTER :: a 00820 00821 CALL timeset(routineN,handle) 00822 00823 failure=.FALSE. 00824 NULLIFY(a) 00825 00826 CPPrecondition(ASSOCIATED(matrix_a),cp_failure_level,routineP,error,failure) 00827 CPPrecondition(matrix_a%ref_count>0,cp_failure_level,routineP,error,failure) 00828 00829 a => matrix_a%local_data 00830 size_a = SIZE(a,1)*SIZE(a,2) 00831 00832 CALL DSCAL(size_a, alpha, a, 1) 00833 00834 CALL timestop(handle) 00835 00836 END SUBROUTINE cp_fm_scale 00837 00838 ! ***************************************************************************** 00845 SUBROUTINE cp_fm_transpose(matrix,matrixt,error) 00846 TYPE(cp_fm_type), POINTER :: matrix, matrixt 00847 TYPE(cp_error_type), INTENT(inout) :: error 00848 00849 CHARACTER(len=*), PARAMETER :: routineN = 'cp_fm_transpose', 00850 routineP = moduleN//':'//routineN 00851 00852 INTEGER :: handle, i, j, ncol_global, 00853 nrow_global 00854 INTEGER, DIMENSION(9) :: desca, descc 00855 LOGICAL :: failure 00856 REAL(KIND=dp), DIMENSION(:, :), POINTER :: a, c 00857 00858 failure = .FALSE. 00859 CPPrecondition(ASSOCIATED(matrix),cp_failure_level,routineP,error,failure) 00860 CPPrecondition(ASSOCIATED(matrixt),cp_failure_level,routineP,error,failure) 00861 IF (failure) RETURN 00862 nrow_global = matrix%matrix_struct%nrow_global 00863 ncol_global = matrix%matrix_struct%ncol_global 00864 CPPrecondition(nrow_global==ncol_global,cp_failure_level,routineP,error,failure) 00865 IF (failure) RETURN 00866 00867 CALL timeset(routineN,handle) 00868 00869 a => matrix%local_data 00870 c => matrixt%local_data 00871 00872 #if defined(__SCALAPACK) 00873 desca(:) = matrix%matrix_struct%descriptor(:) 00874 descc(:) = matrixt%matrix_struct%descriptor(:) 00875 CALL pdtran(nrow_global,ncol_global,1.0_dp,a(1,1),1,1,desca,0.0_dp,c(1,1),1,1,descc) 00876 #else 00877 DO j=1,ncol_global 00878 DO i=1,nrow_global 00879 c(j,i)=a(i,j) 00880 ENDDO 00881 ENDDO 00882 #endif 00883 CALL timestop(handle) 00884 00885 END SUBROUTINE cp_fm_transpose 00886 00887 ! ***************************************************************************** 00895 SUBROUTINE cp_fm_upper_to_full(matrix,work,error) 00896 00897 TYPE(cp_fm_type), POINTER :: matrix,work 00898 TYPE(cp_error_type), INTENT(inout) :: error 00899 00900 ! *** Local variables *** 00901 00902 INTEGER :: handle,icol_global,icol_local,ipcol,iprow,irow_global, 00903 irow_local,mypcol,myprow,ncol_block,ncol_global,ncol_local, 00904 npcol,nprow,nrow_block,nrow_global,nrow_local 00905 00906 INTEGER, DIMENSION(9) :: desca,descc 00907 00908 REAL(KIND = dp), DIMENSION(:,:), POINTER :: a,c 00909 REAL(KIND = sp), DIMENSION(:,:), POINTER :: a_sp,c_sp 00910 TYPE(cp_blacs_env_type), POINTER :: context 00911 CHARACTER(len=*), PARAMETER :: routineN='cp_fm_upper_to_full', 00912 routineP=moduleN//':'//routineN 00913 LOGICAL :: failure 00914 00915 #if defined(__SCALAPACK) 00916 INTEGER, EXTERNAL :: indxl2g 00917 00918 #endif 00919 00920 failure = .FALSE. 00921 CPPrecondition(ASSOCIATED(matrix),cp_failure_level,routineP,error,failure) 00922 CPPrecondition(ASSOCIATED(work),cp_failure_level,routineP,error,failure) 00923 IF (failure) RETURN 00924 nrow_global = matrix%matrix_struct%nrow_global 00925 ncol_global = matrix%matrix_struct%ncol_global 00926 CPPrecondition(nrow_global==ncol_global,cp_failure_level,routineP,error,failure) 00927 nrow_global = work%matrix_struct%nrow_global 00928 ncol_global = work%matrix_struct%ncol_global 00929 CPPrecondition(nrow_global==ncol_global,cp_failure_level,routineP,error,failure) 00930 CPPrecondition(matrix%use_sp.EQV.work%use_sp,cp_failure_level,routineP,error,failure) 00931 IF (failure) RETURN 00932 00933 CALL timeset(routineN,handle) 00934 00935 context => matrix%matrix_struct%context 00936 myprow=context%mepos(1) 00937 mypcol=context%mepos(2) 00938 nprow=context%num_pe(1) 00939 npcol=context%num_pe(2) 00940 00941 #if defined(__SCALAPACK) 00942 00943 nrow_block = matrix%matrix_struct%nrow_block 00944 ncol_block = matrix%matrix_struct%ncol_block 00945 00946 nrow_local = matrix%matrix_struct%nrow_locals(myprow) 00947 ncol_local = matrix%matrix_struct%ncol_locals(mypcol) 00948 00949 a => work%local_data 00950 a_sp => work%local_data_sp 00951 desca(:) = work%matrix_struct%descriptor(:) 00952 c => matrix%local_data 00953 c_sp => matrix%local_data_sp 00954 descc(:) = matrix%matrix_struct%descriptor(:) 00955 00956 DO icol_local=1,ncol_local 00957 icol_global = indxl2g(icol_local,ncol_block,mypcol,& 00958 matrix%matrix_struct%first_p_pos(2),npcol) 00959 DO irow_local=1,nrow_local 00960 irow_global = indxl2g(irow_local,nrow_block,myprow,& 00961 matrix%matrix_struct%first_p_pos(1),nprow) 00962 IF (irow_global > icol_global) THEN 00963 IF(matrix%use_sp) THEN 00964 c_sp(irow_local,icol_local) = 0.0_sp 00965 ELSE 00966 c(irow_local,icol_local) = 0.0_dp 00967 ENDIF 00968 ELSE IF (irow_global == icol_global) THEN 00969 IF(matrix%use_sp) THEN 00970 c_sp(irow_local,icol_local) = 0.5_sp*c_sp(irow_local,icol_local) 00971 ELSE 00972 c(irow_local,icol_local) = 0.5_dp*c(irow_local,icol_local) 00973 ENDIF 00974 END IF 00975 END DO 00976 END DO 00977 00978 DO icol_local=1,ncol_local 00979 DO irow_local=1,nrow_local 00980 IF(matrix%use_sp) THEN 00981 a_sp(irow_local,icol_local) = c_sp(irow_local,icol_local) 00982 ELSE 00983 a(irow_local,icol_local) = c(irow_local,icol_local) 00984 ENDIF 00985 ENDDO 00986 ENDDO 00987 00988 IF(matrix%use_sp) THEN 00989 CALL pstran(nrow_global,ncol_global,1.0_sp,a_sp(1,1),1,1,desca,1.0_sp,c_sp(1,1),1,1,descc) 00990 ELSE 00991 CALL pdtran(nrow_global,ncol_global,1.0_dp,a(1,1),1,1,desca,1.0_dp,c(1,1),1,1,descc) 00992 ENDIF 00993 00994 #else 00995 00996 a => matrix%local_data 00997 a_sp => matrix%local_data_sp 00998 DO irow_global=1,nrow_global 00999 DO icol_global=irow_global+1,ncol_global 01000 IF(matrix%use_sp) THEN 01001 a_sp(icol_global,irow_global)=a_sp(irow_global,icol_global) 01002 ELSE 01003 a(icol_global,irow_global)=a(irow_global,icol_global) 01004 ENDIF 01005 ENDDO 01006 ENDDO 01007 01008 #endif 01009 CALL timestop(handle) 01010 01011 END SUBROUTINE cp_fm_upper_to_full 01012 01013 ! ***************************************************************************** 01022 SUBROUTINE cp_fm_column_scale(matrixa,scaling) 01023 TYPE(cp_fm_type), POINTER :: matrixa 01024 REAL(KIND=dp), DIMENSION(:), INTENT(in) :: scaling 01025 01026 INTEGER :: i, icol_global, icol_local, ipcol, iprow, irow_local, k, 01027 mypcol, myprow, n, ncol_global, npcol, nprow 01028 REAL(KIND=dp), DIMENSION(:, :), POINTER :: a 01029 REAL(KIND=sp), DIMENSION(:, :), POINTER :: a_sp 01030 01031 myprow=matrixa%matrix_struct%context%mepos(1) 01032 mypcol=matrixa%matrix_struct%context%mepos(2) 01033 nprow=matrixa%matrix_struct%context%num_pe(1) 01034 npcol=matrixa%matrix_struct%context%num_pe(2) 01035 01036 ncol_global =matrixa%matrix_struct%ncol_global 01037 01038 a => matrixa%local_data 01039 a_sp => matrixa%local_data_sp 01040 IF(matrixa%use_sp) THEN 01041 n = SIZE(a_sp,1) 01042 ELSE 01043 n = SIZE(a,1) 01044 ENDIF 01045 k = MIN(SIZE(scaling),ncol_global) 01046 01047 #if defined(__SCALAPACK) 01048 01049 DO icol_global=1,k 01050 CALL infog2l(1,icol_global,matrixa%matrix_struct%descriptor,& 01051 nprow,npcol,myprow,mypcol,& 01052 irow_local,icol_local,iprow,ipcol) 01053 IF ((ipcol == mypcol)) THEN 01054 IF(matrixa%use_sp) THEN 01055 CALL SSCAL(n,REAL(scaling(icol_global),sp),a_sp(1,icol_local),1) 01056 ELSE 01057 CALL DSCAL(n,scaling(icol_global),a(1,icol_local),1) 01058 ENDIF 01059 END IF 01060 ENDDO 01061 #else 01062 DO i=1,k 01063 IF(matrixa%use_sp) THEN 01064 CALL SSCAL(n,REAL(scaling(i),sp),a_sp(1,i),1) 01065 ELSE 01066 CALL DSCAL(n,scaling(i),a(1,i),1) 01067 ENDIF 01068 ENDDO 01069 #endif 01070 END SUBROUTINE cp_fm_column_scale 01071 01072 ! ***************************************************************************** 01081 SUBROUTINE cp_fm_invert(matrix_a,matrix_inverse,det_a,error) 01082 01083 TYPE(cp_fm_type), POINTER :: matrix_a, matrix_inverse 01084 REAL(KIND=dp), INTENT(OUT), OPTIONAL :: det_a 01085 TYPE(cp_error_type), INTENT(inout) :: error 01086 01087 CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_fm_invert', 01088 routineP = moduleN//':'//routineN 01089 01090 INTEGER :: i, info, liwork, lwork, n 01091 INTEGER, ALLOCATABLE, DIMENSION(:) :: ipivot, iwork 01092 INTEGER, DIMENSION(9) :: desca 01093 LOGICAL :: failure, sign 01094 REAL(KIND=dp) :: alpha, beta, determinant, eps1 01095 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: berr, ferr, work 01096 REAL(KIND=dp), DIMENSION(:), POINTER :: diag 01097 REAL(KIND=dp), DIMENSION(:, :), POINTER :: a 01098 TYPE(cp_fm_type), POINTER :: matrix_B, matrix_lu 01099 01100 failure=.FALSE. 01101 01102 CALL cp_fm_create(matrix=matrix_lu,& 01103 matrix_struct=matrix_a%matrix_struct,& 01104 name="A_lu"//TRIM(ADJUSTL(cp_to_string(1)))//"MATRIX",& 01105 error=error) 01106 CALL cp_fm_to_fm(matrix_a,matrix_lu,error=error) 01107 01108 CALL cp_fm_create(matrix=matrix_B,& 01109 matrix_struct=matrix_a%matrix_struct,& 01110 name="B_mat"//TRIM(ADJUSTL(cp_to_string(1)))//"MATRIX",& 01111 error=error) 01112 a => matrix_lu%local_data 01113 n = matrix_lu%matrix_struct%nrow_global 01114 ALLOCATE(ipivot(n+matrix_a%matrix_struct%nrow_block)) 01115 #if defined(__SCALAPACK) 01116 ALLOCATE(ferr(n)) 01117 ALLOCATE(berr(n)) 01118 ALLOCATE(work(3*n)) 01119 ALLOCATE(iwork(3*N)) 01120 lwork=3*n 01121 liwork=3*n 01122 desca(:) = matrix_lu%matrix_struct%descriptor(:) 01123 CALL pdgetrf(n,n,a(1,1),1,1,desca,ipivot,info) 01124 ALLOCATE(diag(n)) 01125 diag(:)=0.0_dp 01126 DO i=1,n 01127 CALL cp_fm_get_element(matrix_lu,i,i,diag(i)) ! not completely optimal in speed i would say 01128 ENDDO 01129 determinant=1.0_dp 01130 DO i=1,n 01131 IF(ipivot(i)==i)THEN 01132 determinant=determinant*diag(i) 01133 ELSE 01134 determinant=determinant*diag(i)*(-1.0_dp) 01135 END IF 01136 ENDDO 01137 DEALLOCATE(diag) 01138 01139 alpha=0.0_dp 01140 beta=1.0_dp 01141 CALL cp_fm_set_all(matrix_inverse,alpha,beta,error) 01142 CALL pdgetrs('N',n,n,matrix_lu%local_data,1,1,desca,ipivot,matrix_inverse%local_data,1,1,desca,info) 01143 ! CALL cp_fm_set_all(matrix_B,alpha,beta,error) 01144 ! DO iter=1,10 01145 ! CALL pdgerfs('N',n,n,matrix_a%local_data,1,1,desca,matrix_lu%local_data,& 01146 ! 1,1,desca,ipivot,matrix_B%local_data,& 01147 ! 1,1,desca,matrix_inverse%local_data,1,1,& 01148 ! desca,ferr,berr,work,lwork,iwork,liwork,info) 01149 ! eps1=eps2 01150 ! eps2=MAXVAL(ferr) 01151 ! IF (ABS( eps2 - eps1) <= EPSILON(1.0_dp))THEN 01152 ! EXIT 01153 ! END IF 01154 ! END DO 01155 DEALLOCATE(ferr) 01156 DEALLOCATE(berr) 01157 DEALLOCATE(work) 01158 DEALLOCATE(iwork) 01159 01160 #else 01161 sign=.TRUE. 01162 CALL invert_matrix(matrix_a%local_data,matrix_inverse%local_data,& 01163 eval_error=eps1,error=error) 01164 CALL cp_fm_lu_decompose(matrix_lu,determinant,correct_sign=sign) 01165 #endif 01166 CALL cp_fm_release(matrix_lu,error=error) 01167 CALL cp_fm_release(matrix_B,error=error) 01168 DEALLOCATE(ipivot) 01169 IF(PRESENT(det_a)) det_a = determinant 01170 END SUBROUTINE cp_fm_invert 01171 01172 ! ***************************************************************************** 01176 SUBROUTINE cp_fm_triangular_invert(matrix_a,uplo_tr,error) 01177 01178 TYPE(cp_fm_type), POINTER :: matrix_a 01179 CHARACTER, INTENT(IN), OPTIONAL :: uplo_tr 01180 TYPE(cp_error_type), INTENT(inout) :: error 01181 01182 CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_fm_triangular_invert', 01183 routineP = moduleN//':'//routineN 01184 01185 CHARACTER :: unit_diag, uplo 01186 INTEGER :: handle, info, ncol_global 01187 INTEGER, DIMENSION(9) :: desca 01188 LOGICAL :: failure 01189 REAL(KIND=dp), DIMENSION(:, :), POINTER :: a 01190 01191 failure=.FALSE. 01192 01193 CALL timeset(routineN,handle) 01194 01195 unit_diag='N' 01196 uplo='U' 01197 IF(PRESENT(uplo_tr)) uplo=uplo_tr 01198 01199 ncol_global =matrix_a%matrix_struct%ncol_global 01200 01201 a => matrix_a%local_data 01202 01203 #if defined(__SCALAPACK) 01204 01205 desca(:) = matrix_a%matrix_struct%descriptor(:) 01206 01207 CALL pdtrtri( uplo, unit_diag, ncol_global, a(1,1), 1, 1, desca, info ) 01208 01209 #else 01210 CALL dtrtri( uplo, unit_diag, ncol_global, a(1,1), ncol_global, info ) 01211 #endif 01212 01213 01214 CALL timestop(handle) 01215 END SUBROUTINE cp_fm_triangular_invert 01216 01217 01218 ! ***************************************************************************** 01219 ! ***************************************************************************** 01225 SUBROUTINE cp_fm_qr_factorization(matrix_a, matrix_r, nrow_fact, ncol_fact, first_row, first_col, error) 01226 TYPE(cp_fm_type), POINTER :: matrix_a, matrix_r 01227 INTEGER, INTENT(IN), OPTIONAL :: nrow_fact, ncol_fact, 01228 first_row, first_col 01229 TYPE(cp_error_type), INTENT(inout) :: error 01230 01231 CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_fm_qr_factorization', 01232 routineP = moduleN//':'//routineN 01233 REAL(KIND=dp), PARAMETER :: rone = 1.0_dp, rzero = 0.0_dp 01234 01235 INTEGER :: handle, i, icol, info, irow, 01236 istat, j, lda, lwork, ncol, 01237 ndim, nrow 01238 INTEGER, DIMENSION(9) :: desca 01239 LOGICAL :: failure 01240 REAL(dp), ALLOCATABLE, DIMENSION(:) :: tau, work 01241 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: r_mat 01242 REAL(KIND=dp), DIMENSION(:, :), POINTER :: a 01243 01244 failure=.FALSE. 01245 01246 CALL timeset(routineN,handle) 01247 01248 ncol =matrix_a%matrix_struct%ncol_global 01249 nrow =matrix_a%matrix_struct%nrow_global 01250 lda = nrow 01251 01252 a => matrix_a%local_data 01253 01254 IF(PRESENT(nrow_fact)) nrow = nrow_fact 01255 IF(PRESENT(ncol_fact)) ncol = ncol_fact 01256 irow = 1 01257 IF(PRESENT(first_row)) irow = first_row 01258 icol = 1 01259 IF(PRESENT(first_col)) icol = first_col 01260 01261 01262 CPPrecondition(nrow>=ncol,cp_failure_level,routineP,error,failure) 01263 ndim = SIZE(a,2) 01264 ! ALLOCATE(ipiv(ndim),istat=STAT) 01265 ALLOCATE(tau(ndim),STAT=istat) 01266 01267 #if defined(__SCALAPACK) 01268 01269 desca(:) = matrix_a%matrix_struct%descriptor(:) 01270 01271 lwork = -1 01272 ALLOCATE(work(2*ndim),STAT=istat) 01273 CALL pdgeqrf( nrow, ncol, a, irow, icol, desca, tau, work, lwork, info ) 01274 lwork = work(1) 01275 DEALLOCATE(work,STAT=istat) 01276 ALLOCATE(work(lwork),STAT=istat) 01277 CALL pdgeqrf( nrow, ncol, a, irow, icol, desca, tau, work, lwork, info ) 01278 01279 #else 01280 lwork = -1 01281 ALLOCATE(work(2*ndim),STAT=istat) 01282 CALL dgeqrf( nrow, ncol, a, lda, tau, work, lwork, info ) 01283 lwork = work(1) 01284 DEALLOCATE(work,STAT=istat) 01285 ALLOCATE(work(lwork),STAT=istat) 01286 CALL dgeqrf( nrow, ncol, a, lda, tau, work, lwork, info ) 01287 01288 #endif 01289 01290 ALLOCATE(r_mat(ncol,ncol),STAT=istat) 01291 CALL cp_fm_get_submatrix(matrix_a,r_mat,1,1,ncol,ncol,error=error) 01292 DO i = 1,ncol 01293 DO j = i+1,ncol 01294 r_mat(j,i) = 0.0_dp 01295 END DO 01296 END DO 01297 CALL cp_fm_set_submatrix(matrix_r,r_mat,1,1,ncol,ncol,error=error) 01298 01299 01300 DEALLOCATE(tau, work, r_mat, STAT=istat) 01301 01302 CALL timestop(handle) 01303 01304 END SUBROUTINE cp_fm_qr_factorization 01305 01306 ! ***************************************************************************** 01311 SUBROUTINE cp_fm_solve(matrix_a,general_a,error) 01312 TYPE(cp_fm_type), POINTER :: matrix_a, general_a 01313 TYPE(cp_error_type), INTENT(inout) :: error 01314 01315 CHARACTER(len=*), PARAMETER :: routineN = 'cp_fm_solve', 01316 routineP = moduleN//':'//routineN 01317 REAL(KIND=dp), PARAMETER :: one = 1.0_dp, zero = 0.0_dp 01318 01319 INTEGER :: handle, info, lda, ldb, n 01320 INTEGER, ALLOCATABLE, DIMENSION(:) :: ipivot 01321 INTEGER, DIMENSION(9) :: desca, descb 01322 REAL(KIND=dp), DIMENSION(:, :), POINTER :: a, a_general 01323 01324 ! *** locals *** 01325 01326 CALL timeset(routineN,handle) 01327 01328 a => matrix_a%local_data 01329 a_general => general_a%local_data 01330 n = matrix_a%matrix_struct%nrow_global 01331 ALLOCATE(ipivot(n+matrix_a%matrix_struct%nrow_block)) 01332 01333 #if defined(__SCALAPACK) 01334 desca(:) = matrix_a%matrix_struct%descriptor(:) 01335 descb(:) = general_a%matrix_struct%descriptor(:) 01336 CALL pdgetrf(n,n,a(1,1),1,1,desca,ipivot,info) 01337 CALL pdgetrs("N" , n , n , a(1,1), 1, 1, desca ,ipivot, a_general(1,1) ,& 01338 1, 1, descb, info ) 01339 01340 #else 01341 lda=SIZE(a,1) 01342 ldb=SIZE(a_general,1) 01343 CALL dgetrf(n,n,a(1,1),lda,ipivot,info) 01344 CALL dgetrs("N",n,n,a(1,1),lda,ipivot,a_general,ldb,info) 01345 01346 #endif 01347 ! info is allowed to be zero 01348 ! this does just signal a zero diagonal element 01349 DEALLOCATE(ipivot) 01350 CALL timestop(handle) 01351 END SUBROUTINE 01352 01353 END MODULE cp_fm_basic_linalg 01354
1.7.3