DONT MERGE: temp. revert

Author: stephankramer

femtools/Multigrid.F90

@@ -707,11 +707,7 @@ function Prolongator(A, epsilon, omega, maxclustersize, cluster) result (P)
   PetscInt:: diagminloc
   PetscReal:: diagmin
   Vec:: sqrt_diag, inv_sqrt_diag, diag, one
-#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR<23)
   double precision, dimension(MAT_INFO_SIZE):: matrixinfo
-#else
-    MatInfo :: matrixinfo
-#endif
   integer, dimension(:), allocatable:: findN, N, R
   integer:: nrows, nentries, ncols
   integer:: jc, ccnt, base, end_of_range
@@ -720,11 +716,7 @@ function Prolongator(A, epsilon, omega, maxclustersize, cluster) result (P)
   call MatGetLocalSize(A, nrows, ncols, ierr)
   ! use Petsc_Tools's MatGetInfo because of bug in earlier patch levels of petsc 3.0
   call MatGetInfo(A, MAT_LOCAL, matrixinfo, ierr)
-#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR<23)
   nentries=matrixinfo(MAT_INFO_NZ_USED)
-#else
-  nentries=int(matrixinfo%nz_used)
-#endif
   call MatGetOwnerShipRange(A, base, end_of_range, ierr)
   ! we decrease by 1, so base+i gives 0-based petsc index if i is the local fortran index:
   base=base-1
@@ -866,7 +858,7 @@ subroutine create_prolongator(P, nrows, ncols, findN, N, R, A, base, omega)
     coarse_base=coarse_base-1
   else
     call MatCreateAIJ(MPI_COMM_SELF, nrows, ncols, nrows, ncols, &
-      0, dnnz, 0, PETSC_NULL_INTEGER_ARRAY, P, ierr)
+      0, dnnz, 0, PETSC_NULL_INTEGER, P, ierr)
     call MatSetOption(P, MAT_USE_INODES, PETSC_FALSE, ierr)
     ! subtract 1 from each cluster no to get petsc 0-based numbering
     coarse_base=-1
@@ -913,20 +905,13 @@ subroutine Prolongator_init(R, ccnt, findN, N, A, base, epsilon)
 PetscReal, intent(in):: epsilon
 
   PetscErrorCode:: ierr
-#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR<23)
   PetscReal, dimension(:), allocatable:: vals(:)
   integer, dimension(:), allocatable:: cols(:)
-#else
-  PetscReal, dimension(:), pointer:: vals(:)
-  integer, dimension(:), pointer:: cols(:)
-#endif
   PetscReal aij, eps_sqrt
   integer i, j, k, p, ncols
 
-#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR<23)
   ! workspace for MatGetRow
   allocate( vals(1:size(N)), cols(1:size(n)) )
-#endif
 
   eps_sqrt=sqrt(epsilon)
 

femtools/Petsc_Tools.F90

@@ -1113,7 +1113,7 @@ function CreatePrivateMatrixFromSparsity(sparsity) result (M)
     end do
 
     call MatCreateAIJ(MPI_COMM_SELF, nprows, npcols, nprows, npcols, &
-      0, nnz, 0, PETSC_NULL_INTEGER_ARRAY, M, ierr)
+      0, nnz, 0, PETSC_NULL_INTEGER, M, ierr)
     call MatSetup(M, ierr)
 
     call MatSetOption(M, MAT_USE_INODES, PETSC_FALSE, ierr)
@@ -1325,27 +1325,16 @@ function petsc2csr(matrix, column_numbering) result(A)
 
     PetscErrorCode ierr
     type(csr_sparsity) :: sparsity
-#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR<23)
+    double precision, dimension(MAT_INFO_SIZE):: matrixinfo
     PetscScalar, dimension(:), allocatable:: row_vals
-    integer, dimension(:), allocatable:: row_cols
-    MatInfo, dimension(MAT_INFO_SIZE):: matrixinfo
-#else
-    PetscScalar, dimension(:), pointer:: row_vals
-    integer, dimension(:), pointer :: row_cols
-    MatInfo :: matrixinfo
-#endif
-    integer, dimension(:), allocatable:: unn2gnn
+    integer, dimension(:), allocatable:: row_cols, unn2gnn
     integer private_columns
     integer i, j, k, ui, rows, columns, entries, ncols, offset, end_of_range
     logical parallel
 
     ! get the necessary info about the matrix:
     call MatGetInfo(matrix, MAT_LOCAL, matrixinfo, ierr)
-#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR<23)
-    entries = matrixinfo(MAT_INFO_NZ_USED)
-#else
-    entries = int(matrixinfo%nz_used)
-#endif
+    entries=matrixinfo(MAT_INFO_NZ_USED)
     ! note we're no longer using MAT_INFO for getting local n/o rows and cols
     ! as it's bugged in Petsc < 3.0 and obsoloted thereafter:
     call MatGetLocalSize(matrix, rows, columns, ierr)
@@ -1386,10 +1375,10 @@ function petsc2csr(matrix, column_numbering) result(A)
     end if
     call allocate(A, sparsity)
 
-#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR<23)
-    allocate(row_cols(1:entries), row_vals(1:entries))
-#endif
     if (parallel) then
+      ! in the parallel case we first copy in a temp. buffer
+      ! and only insert the entries A_ij with both i *and* j local
+      allocate(row_cols(1:entries), row_vals(1:entries))
       j=1
       do i=0, rows-1
         sparsity%findrm(i+1)=j
@@ -1408,30 +1397,41 @@ function petsc2csr(matrix, column_numbering) result(A)
             j=j+1
           end if
         end do
+        ! This is stupid, we were given copies in MatGetRow so it could
+        ! have restored its internal tmp arrays straight away, anyway:
         call MatRestoreRow(matrix, offset+i, ncols, row_cols, row_vals, ierr)
       end do
       A%sparsity%findrm(i+1)=j
 
+      deallocate(row_cols, row_vals)
     else
       ! Serial case:
       j=1
       do i=0, rows-1
         sparsity%findrm(i+1)=j
-        call MatGetRow(matrix, offset+i, ncols, row_cols, row_vals, ierr)
-        sparsity%colm(j:j+ncols-1) = row_cols(1:ncols)
-        A%val(j:j+ncols-1) = row_vals(1:ncols)
+#ifdef DOUBLEP
+        call MatGetRow(matrix, offset+i, ncols, sparsity%colm(j:), A%val(j:), ierr)
         j=j+ncols
-        call MatRestoreRow(matrix, offset+i, ncols, row_cols, row_vals, ierr)
+        ! This is stupid, we were given copies in MatGetRow so it could
+        ! have restored its internal tmp arrays straight away, anyway:
+        call MatRestoreRow(matrix, offset+i, ncols, sparsity%colm(j:), A%val(j:), ierr)
+#else
+        allocate(row_vals(size(A%val) - j + 1))
+        call MatGetRow(matrix, offset+i, ncols, sparsity%colm(j:), row_vals, ierr)
+        A%val(j:) = row_vals
+        j=j+ncols
+        ! This is stupid, we were given copies in MatGetRow so it could
+        ! have restored its internal tmp arrays straight away, anyway:
+        call MatRestoreRow(matrix, offset+i, ncols, sparsity%colm(j:), row_vals, ierr)
+        deallocate(row_vals)
+#endif
       end do
       A%sparsity%findrm(i+1)=j
 
       ! matrix is indexed from offset+0, colm should be indexed from 1:
       sparsity%colm=sparsity%colm-offset+1
 
     end if
-#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR<23)
-    deallocate(row_cols, row_vals)
-#endif
     call deallocate(sparsity)
 
     if (present(column_numbering)) then
@@ -1535,7 +1535,7 @@ function IsNullMatNullSpace(nullsp)
     ! MatNullSpace(0) is what is returned by MatGetNullspace if no nullspace is present
     ! (because a wrapper on the output is missing, and there isn't a PETSC_NULL_MATNULLSPACE
     ! in the first place)
-    IsNullMatNullSpace = nullsp%v<= 0
+    IsNullMatNullSpace = nullsp%v==-1 .or. nullsp%v==0
 
 end function IsNullMatNullSpace
 

femtools/Solvers.F90

@@ -1366,16 +1366,15 @@ subroutine petsc_solve_destroy_petsc_csr(y, b, solver_option_path)
 
 end subroutine petsc_solve_destroy_petsc_csr
 
-subroutine ConvergenceCheck(converged_reason, iterations, name, solver_option_path, &
+subroutine ConvergenceCheck(reason, iterations, name, solver_option_path, &
   startfromzero, A, b, petsc_numbering, x0, vector_x0, checkconvergence, nomatrixdump)
   !!< Checks reason of convergence. If negative (not converged)
   !!< writes out a scary warning and dumps matrix (if first time),
   !!< and if reason<0 but reason/=-3
   !!< (i.e. not converged due to other reasons than reaching max_its)
   !!< it sets sig_int to .true. causing the run to halt and dump
   !!< at the end of the time step.
-  KSPConvergedReason, intent(in):: converged_reason
-  integer, intent(in):: iterations
+  integer, intent(in):: reason, iterations
   !! name of the thing we're solving for, used in log output:
   character(len=*), intent(in):: name
   !! for new options path to solver options
@@ -1399,7 +1398,6 @@ subroutine ConvergenceCheck(converged_reason, iterations, name, solver_option_pa
   PetscErrorCode ierr
   character(len=30) reasons(10)
   real spin_up_time, current_time
-  integer reason
 
   reasons(1)  = "Undefined"
   reasons(2)  = "KSP_DIVERGED_NULL"
@@ -1412,11 +1410,6 @@ subroutine ConvergenceCheck(converged_reason, iterations, name, solver_option_pa
   reasons(9)  = "KSP_DIVERGED_NAN"
   reasons(10) = "KSP_DIVERGED_INDEFINITE_MAT"
 
-#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR<23)
-  reason = converged_reason
-#else
-  reason = converged_reason%v  ! recover the (old style) numerical version (probably a bad idea)
-#endif
   if (reason<=0) then
      if(present_and_true(nomatrixdump)) matrixdumped = .true.
      if (present(checkconvergence)) then
@@ -1822,11 +1815,6 @@ recursive subroutine setup_pc_from_options(pc, pmat, option_path, &
   logical, optional, intent(in) :: is_subpc
 
     KSP:: subksp
-#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR<23)
-    KSP:: subksps
-#else
-    KSP, dimension(:), pointer :: subksps
-#endif
     PC:: subpc
     MatNullSpace:: nullsp
     PCType:: pctype, hypretype
@@ -1894,19 +1882,11 @@ recursive subroutine setup_pc_from_options(pc, pmat, option_path, &
       ! need to call this before the subpc can be retrieved:
       call PCSetup(pc, ierr)
 
-#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>=23)
-      nullify(subksps)
-#endif
       if (pctype==PCBJACOBI) then
-        call PCBJACOBIGetSubKSP(pc, n_local, first_local, subksps, ierr)
+        call PCBJACOBIGetSubKSP(pc, n_local, first_local, subksp, ierr)
       else
-        call PCASMGetSubKSP(pc, n_local, first_local, subksps, ierr)
+        call PCASMGetSubKSP(pc, n_local, first_local, subksp, ierr)
       end if
-#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>=23)
-      subksp = subksps(1)
-#else
-      subksp = subksps
-#endif
 
       call KSPGetPC(subksp, subpc, ierr)
       ! recursively call to setup the subpc
@@ -1928,15 +1908,7 @@ recursive subroutine setup_pc_from_options(pc, pmat, option_path, &
        call PCSetType(pc, PCBJACOBI, ierr)
        ! need to call this before the subpc can be retrieved:
        call PCSetup(pc, ierr)
-#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>=23)
-       nullify(subksps)
-#endif
-       call PCBJACOBIGetSubKSP(pc, n_local, first_local, subksps, ierr)
-#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>=23)
-       subksp = subksps(1)
-#else
-       subksp = subksps
-#endif
+       call PCBJACOBIGetSubKSP(pc, n_local, first_local, subksp, ierr)
        call KSPGetPC(subksp, subpc, ierr)
        call PCSetType(subpc, pctype, ierr)
 
@@ -2012,11 +1984,7 @@ recursive subroutine setup_fieldsplit_preconditioner(pc, option_path, &
   type(petsc_numbering_type), intent(in):: petsc_numbering
 
     character(len=128):: fieldsplit_type
-#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR<23)
     KSP, dimension(size(petsc_numbering%gnn2unn,2)):: subksps
-#else
-    KSP, dimension(:), pointer:: subksps => null()
-#endif
     Mat :: mat, pmat
     MatNullSpace :: null_space
     IS:: index_set

femtools/Sparse_Tools_Petsc.F90

@@ -411,7 +411,7 @@ subroutine allocate_petsc_csr_matrix_from_nnz(matrix, rows, columns, &
       ! Create serial block matrix:
       call MatCreateBAIJ(MPI_COMM_SELF, element_size, &
          urows, ucols, urows, ucols, &
-         0, dnnz, 0, PETSC_NULL_INTEGER_ARRAY, matrix%M, ierr)
+         0, dnnz, 0, PETSC_NULL_INTEGER, matrix%M, ierr)
 
     elseif (use_element_blocks) then
 
@@ -429,7 +429,7 @@ subroutine allocate_petsc_csr_matrix_from_nnz(matrix, rows, columns, &
 
       ! Create serial matrix:
       call MatCreateAIJ(MPI_COMM_SELF, urows, ucols, urows, ucols, &
-         0, dnnz, 0, PETSC_NULL_INTEGER_ARRAY, matrix%M, ierr)
+         0, dnnz, 0, PETSC_NULL_INTEGER, matrix%M, ierr)
       call MatSetBlockSizes(matrix%M, lgroup_size(1), lgroup_size(2), ierr)
 
     else
@@ -735,20 +735,12 @@ function petsc_csr_entries(matrix) result (entries)
     integer :: entries
     type(petsc_csr_matrix), intent(in) :: matrix
 
-#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR<23)
     double precision, dimension(MAT_INFO_SIZE):: matrixinfo
-#else
-    MatInfo :: matrixinfo
-#endif
     PetscErrorCode:: ierr
 
     ! get the necessary info about the matrix:
     call MatGetInfo(matrix%M, MAT_LOCAL, matrixinfo, ierr)
-#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR<23)
     entries=matrixinfo(MAT_INFO_NZ_USED)
-#else
-    entries=int(matrixinfo%nz_used)
-#endif
 
   end function petsc_csr_entries
 
@@ -780,26 +772,6 @@ subroutine petsc_csr_addto(matrix, blocki, blockj, i, j, val)
 
   end subroutine petsc_csr_addto
 
-  subroutine fix_column_numbers(idxm, idxn)
-    ! for some rows that we are not assembling (indicated with negative index)
-    ! the column indices are not reliable, as they are outside any halo and thus
-    ! have not been filled in in the gnn2unn numbering (we don't know their universal petsc number)
-    ! this is fine because we're telling petsc to skip these entries (with a negative row index)
-    ! but the random column index may still trip some internal petsc checks
-    ! therefore here we set these to negative as well
-
-    PetscInt, dimension(:), intent(in) :: idxm  ! row numbers
-    PetscInt, dimension(:), intent(inout) :: idxn  ! column numbers
-    integer i
-
-    do i=1, size(idxm)
-      if (idxm(i)<0) then
-        idxn(i) = -1
-      end if
-    end do
-
-  end subroutine fix_column_numbers
-
   subroutine petsc_csr_vaddto(matrix, blocki, blockj, i, j, val)
     !!< Add multiple values to matrix(blocki, blockj, i,j)
     type(petsc_csr_matrix), intent(inout) :: matrix
@@ -813,10 +785,8 @@ subroutine petsc_csr_vaddto(matrix, blocki, blockj, i, j, val)
 
     idxm=matrix%row_numbering%gnn2unn(i,blocki)
     idxn=matrix%column_numbering%gnn2unn(j,blockj)
-    call fix_column_numbers(idxm, idxn)
 
-    call MatSetValues(matrix%M, size(i), idxm, size(j), idxn, &
-        real(reshape(val, (/ size(val) /)), kind=PetscScalar_kind), &
+    call MatSetValues(matrix%M, size(i), idxm, size(j), idxn, real(val, kind=PetscScalar_kind), &
         ADD_VALUES, ierr)
 
     matrix%is_assembled=.false.
@@ -837,10 +807,9 @@ subroutine petsc_csr_block_addto(matrix, i, j, val)
 
     idxm=matrix%row_numbering%gnn2unn(i,:)
     idxn=matrix%column_numbering%gnn2unn(j,:)
-    call fix_column_numbers(idxm, idxn)
 
     call MatSetValues(matrix%M, size(idxm), idxm, size(idxn), idxn, &
-                  real(reshape(val, (/ size(val) /)), kind=PetscScalar_kind), ADD_VALUES, ierr)
+                  real(val, kind=PetscScalar_kind), ADD_VALUES, ierr)
 
     matrix%is_assembled=.false.
 
@@ -855,7 +824,7 @@ subroutine petsc_csr_blocks_addto(matrix, i, j, val)
     integer, dimension(:), intent(in) :: j
     real, dimension(:,:,:,:), intent(in) :: val
 
-    PetscScalar, dimension(size(i)*size(j)):: value
+    PetscScalar, dimension(size(i), size(j)):: value
     PetscInt, dimension(size(i)):: idxm
     PetscInt, dimension(size(j)):: idxn
     PetscErrorCode:: ierr
@@ -865,9 +834,8 @@ subroutine petsc_csr_blocks_addto(matrix, i, j, val)
       idxm=matrix%row_numbering%gnn2unn(i,blocki)
       do blockj=1, size(matrix%column_numbering%gnn2unn,2)
         idxn=matrix%column_numbering%gnn2unn(j,blockj)
-        call fix_column_numbers(idxm, idxn)
         ! unfortunately we need a copy here to pass contiguous memory
-        value=reshape(val(blocki, blockj, :, :), (/ size(value) /))
+        value=val(blocki, blockj, :, :)
         call MatSetValues(matrix%M, size(i), idxm, size(j), idxn, &
               value, ADD_VALUES, ierr)
       end do
@@ -887,7 +855,7 @@ subroutine petsc_csr_blocks_addto_withmask(matrix, i, j, val, block_mask)
     real, dimension(:,:,:,:), intent(in) :: val
     logical, dimension(:,:), intent(in) :: block_mask
 
-    PetscScalar, dimension(size(i)*size(j)):: value
+    PetscScalar, dimension(size(i), size(j)):: value
     PetscInt, dimension(size(i)):: idxm
     PetscInt, dimension(size(j)):: idxn
     PetscErrorCode:: ierr
@@ -898,9 +866,8 @@ subroutine petsc_csr_blocks_addto_withmask(matrix, i, j, val, block_mask)
       do blockj=1, size(matrix%column_numbering%gnn2unn,2)
         if (block_mask(blocki,blockj)) then
           idxn=matrix%column_numbering%gnn2unn(j,blockj)
-          call fix_column_numbers(idxm, idxn)
           ! unfortunately we need a copy here to pass contiguous memory
-          value=reshape(val(blocki, blockj, :, :), (/ size(value) /))
+          value=val(blocki, blockj, :, :)
           call MatSetValues(matrix%M, size(i), idxm, size(j), idxn, &
                 value, ADD_VALUES, ierr)
         end if