1C 2C This file is part of MUMPS 5.1.2, released 3C on Mon Oct 2 07:37:01 UTC 2017 4C 5C 6C Copyright 1991-2017 CERFACS, CNRS, ENS Lyon, INP Toulouse, Inria, 7C University of Bordeaux. 8C 9C This version of MUMPS is provided to you free of charge. It is 10C released under the CeCILL-C license: 11C http://www.cecill.info/licences/Licence_CeCILL-C_V1-en.html 12C 13 SUBROUTINE CMUMPS_FACTO_ROOT( MYID, MASTER_OF_ROOT, 14 & root, N, IROOT, 15 & COMM, IW, LIW, IFREE, 16 & A, LA, PTRAST, PTLUST_S, PTRFAC, 17 & STEP, INFO, LDLT, QR, 18 & WK, LWK, KEEP,KEEP8,DKEEP,OPELIW) 19 USE CMUMPS_LR_STATS, ONLY: UPDATE_FLOPS_STATS_ROOT 20 IMPLICIT NONE 21 INCLUDE 'cmumps_root.h' 22 INCLUDE 'mpif.h' 23 TYPE ( CMUMPS_ROOT_STRUC ) :: root 24 INTEGER N, IROOT, COMM, LIW, MYID, IFREE, MASTER_OF_ROOT 25 INTEGER(8) :: LA 26 INTEGER(8) :: LWK 27 COMPLEX WK( LWK ) 28 INTEGER KEEP(500) 29 REAL DKEEP(230) 30 INTEGER(8) KEEP8(150) 31 INTEGER(8) :: PTRAST(KEEP(28)) 32 INTEGER(8) :: PTRFAC(KEEP(28)) 33 INTEGER PTLUST_S(KEEP(28)), STEP(N), IW( LIW ) 34 INTEGER INFO( 2 ), LDLT, QR 35 COMPLEX A( LA ) 36 DOUBLE PRECISION, intent(inout) :: OPELIW 37 INTEGER IOLDPS 38 INTEGER(8) :: IAPOS 39 INTEGER(8) :: ENTRIES_ROOT 40 INTEGER LOCAL_M, LOCAL_N, LPIV, IERR, allocok 41 INTEGER FWD_LOCAL_N_RHS, FWD_MTYPE 42 INCLUDE 'mumps_headers.h' 43 EXTERNAL numroc 44 INTEGER numroc 45 IF ( .NOT. root%yes ) RETURN 46 IF ( KEEP(60) .NE. 0 ) THEN 47 IF ((LDLT == 1 .OR. LDLT == 2) .AND. KEEP(60) == 3 ) THEN 48 CALL CMUMPS_SYMMETRIZE( WK, root%MBLOCK, 49 & root%MYROW, root%MYCOL, root%NPROW, root%NPCOL, 50 & root%SCHUR_POINTER(1), 51 & root%SCHUR_LLD, root%SCHUR_NLOC, 52 & root%TOT_ROOT_SIZE, MYID, COMM ) 53 ENDIF 54 RETURN 55 ENDIF 56 IOLDPS = PTLUST_S(STEP(IROOT))+KEEP(IXSZ) 57 IAPOS = PTRAST(STEP(IROOT)) 58 LOCAL_M = IW( IOLDPS + 2 ) 59 LOCAL_N = IW( IOLDPS + 1 ) 60 IAPOS = PTRFAC(IW ( IOLDPS + 4 )) 61 IF ( LDLT.EQ.0 .OR. LDLT.EQ.2 .OR. QR.ne.0 ) THEN 62 LPIV = LOCAL_M + root%MBLOCK 63 ELSE 64 LPIV = 1 65 END IF 66 IF (associated( root%IPIV )) DEALLOCATE(root%IPIV) 67 root%LPIV = LPIV 68 ALLOCATE( root%IPIV( LPIV ), stat = allocok ) 69 IF ( allocok .GT. 0 ) THEN 70 INFO(1) = -13 71 INFO(2) = LPIV 72 WRITE(*,*) MYID,': problem allocating IPIV(',LPIV,') in root' 73 CALL MUMPS_ABORT() 74 END IF 75 CALL DESCINIT( root%DESCRIPTOR(1), root%TOT_ROOT_SIZE, 76 & root%TOT_ROOT_SIZE, root%MBLOCK, root%NBLOCK, 77 & 0, 0, root%CNTXT_BLACS, LOCAL_M, IERR ) 78 IF ( LDLT.EQ.2 ) THEN 79 IF(root%MBLOCK.NE.root%NBLOCK) THEN 80 WRITE(*,*) ' Error: symmetrization only works for' 81 WRITE(*,*) ' square block sizes, MBLOCK/NBLOCK=', 82 & root%MBLOCK, root%NBLOCK 83 CALL MUMPS_ABORT() 84 END IF 85 IF ( LWK .LT. min( 86 & int(root%MBLOCK,8) * int(root%NBLOCK,8), 87 & int(root%TOT_ROOT_SIZE,8)* int(root%TOT_ROOT_SIZE,8 ) 88 & )) THEN 89 WRITE(*,*) 'Not enough workspace for symmetrization.' 90 CALL MUMPS_ABORT() 91 END IF 92 CALL CMUMPS_SYMMETRIZE( WK, root%MBLOCK, 93 & root%MYROW, root%MYCOL, root%NPROW, root%NPCOL, 94 & A( IAPOS ), LOCAL_M, LOCAL_N, 95 & root%TOT_ROOT_SIZE, MYID, COMM ) 96 END IF 97 IF (LDLT.EQ.0.OR.LDLT.EQ.2) THEN 98 CALL pcgetrf( root%TOT_ROOT_SIZE, root%TOT_ROOT_SIZE, 99 & A( IAPOS ), 100 & 1, 1, root%DESCRIPTOR(1), root%IPIV(1), IERR ) 101 IF ( IERR .GT. 0 ) THEN 102 INFO(1)=-10 103 INFO(2)=IERR-1 104 END IF 105 ELSE 106 CALL pcpotrf('L',root%TOT_ROOT_SIZE,A(IAPOS), 107 & 1,1,root%DESCRIPTOR(1),IERR) 108 IF ( IERR .GT. 0 ) THEN 109 INFO(1)=-40 110 INFO(2)=IERR-1 111 END IF 112 END IF 113 IF (IERR .GT. 0) THEN 114 CALL MUMPS_UPDATE_FLOPS_ROOT( OPELIW, LDLT, 115 & root%TOT_ROOT_SIZE, INFO(2), 116 & root%NPROW, root%NPCOL, MYID ) 117 IF (KEEP(486) .GT. 0) THEN 118 CALL UPDATE_FLOPS_STATS_ROOT( LDLT, 119 & root%TOT_ROOT_SIZE, INFO(2), 120 & root%NPROW, root%NPCOL, MYID ) 121 ENDIF 122 ELSE 123 CALL MUMPS_UPDATE_FLOPS_ROOT( OPELIW, LDLT, 124 & root%TOT_ROOT_SIZE, root%TOT_ROOT_SIZE, 125 & root%NPROW, root%NPCOL, MYID ) 126 IF (KEEP(486) .GT. 0) THEN 127 CALL UPDATE_FLOPS_STATS_ROOT( LDLT, 128 & root%TOT_ROOT_SIZE, root%TOT_ROOT_SIZE, 129 & root%NPROW, root%NPCOL, MYID ) 130 ENDIF 131 ENDIF 132 IF ( LDLT .EQ. 0 ) THEN 133 ENTRIES_ROOT = int(root%TOT_ROOT_SIZE,8) 134 & * int(root%TOT_ROOT_SIZE,8) 135 ELSE 136 ENTRIES_ROOT = int(root%TOT_ROOT_SIZE,8) 137 & * int(root%TOT_ROOT_SIZE,8) 138 ENDIF 139 KEEP8(10)=KEEP8(10) + ENTRIES_ROOT / 140 & int(root%NPROW * root%NPCOL,8) 141 IF (MYID .eq. MASTER_OF_ROOT) THEN 142 KEEP8(10)=KEEP8(10) + 143 & mod(ENTRIES_ROOT, int(root%NPROW*root%NPCOL,8)) 144 ENDIF 145 IF (KEEP(258).NE.0) THEN 146 IF (root%MBLOCK.NE.root%NBLOCK) THEN 147 write(*,*) "Internal error in CMUMPS_FACTO_ROOT:", 148 & "Block size different for rows and columns", 149 & root%MBLOCK, root%NBLOCK 150 CALL MUMPS_ABORT() 151 ENDIF 152 CALL CMUMPS_GETDETER2D(root%MBLOCK, root%IPIV(1),root%MYROW, 153 & root%MYCOL, root%NPROW, root%NPCOL, A(IAPOS), LOCAL_M, 154 & LOCAL_N, root%TOT_ROOT_SIZE, MYID, DKEEP(6), KEEP(259), 155 & LDLT) 156 ENDIF 157 IF (KEEP(252) .NE. 0) THEN 158 FWD_LOCAL_N_RHS = numroc(KEEP(253), root%NBLOCK, 159 & root%MYCOL, 0, root%NPCOL) 160 FWD_LOCAL_N_RHS = max(1,FWD_LOCAL_N_RHS) 161 FWD_MTYPE = 1 162 CALL CMUMPS_SOLVE_2D_BCYCLIC( 163 & root%TOT_ROOT_SIZE, 164 & KEEP(253), 165 & FWD_MTYPE, 166 & A(IAPOS), 167 & root%DESCRIPTOR(1), 168 & LOCAL_M, LOCAL_N, FWD_LOCAL_N_RHS, 169 & root%IPIV(1), LPIV, 170 & root%RHS_ROOT(1,1), LDLT, 171 & root%MBLOCK, root%NBLOCK, 172 & root%CNTXT_BLACS, IERR) 173 ENDIF 174 RETURN 175 END SUBROUTINE CMUMPS_FACTO_ROOT 176