[petsc-dev] petsc's hypre interface

[Alexander Grayver]

Dear petsc devs,

I have came across strange behavior trying to make use of PETSc routines 
(from mhyp.c) to convert /rectangular/ matrices into hypre format in 
parallel (i.e. MATMPIAIJ).
Whenever I use MatHYPRE_IJMatrixCreate -> MatHYPRE_IJMatrixCopy -> 
hypre_ParCSRMatrixTranspose I get following error:
hypre_assert failed: AT_i[A_j[j]] < num_nonzerosAT
hypre_assert failed: AT_i[A_j[j]] < num_nonzerosAT

it comes from csr_matop.c : hypre_CSRMatrixTranspose.

The error message never appears if construct hypre matrix directly and 
transpose it.

Attached is a minimal example that reproduces the problem (try to 
(un)comment line 58).

I am wondering where is an origin of this error: in my code, in petsc or 
in hypre?

Thanks,
Alexander
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20130812/7a508560/attachment.html>
-------------- next part --------------
#include <petscmat.h>
#include "../src/dm/impls/da/hypre/mhyp.h"

#undef __FUNCT__
#define __FUNCT__ "main"
int main(int argc,char **args)
{
  Mat A;  
  PetscInt ierr;
  PetscInt M = 1000,N = 500,ilower,iupper,jlower,jupper,i,j[2];
  PetscScalar v[2] = {1, 1};
  
  HYPRE_IJMatrix B, C;
  HYPRE_ParCSRMatrix par_B, par_BT, par_C, par_CT;
  HYPRE_Int herr, data = 0;
  HYPRE_Int nnz = 2;
  
  PetscInitialize(&argc,&args,(char*)0,PETSC_NULL);
  
  ierr = MatCreateAIJ(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,M,N,2,PETSC_NULL,2,PETSC_NULL,&A);CHKERRQ(ierr);
  
  MatGetOwnershipRange(A, &ilower, &iupper);
  MatGetOwnershipRangeColumn(A, &jlower, &jupper);
  
  herr = HYPRE_IJMatrixCreate(PETSC_COMM_WORLD, ilower, iupper - 1, jlower, jupper - 1, &C);
  if(herr) printf("Error in HYPRE_IJMatrixCreate!");
  herr = HYPRE_IJMatrixSetObjectType(C, HYPRE_PARCSR);
  if(herr) printf("Error in HYPRE_IJMatrixSetObjectType!");
  herr = HYPRE_IJMatrixInitialize(C);
  if(herr) printf("Error in HYPRE_IJMatrixInitialize!");
    
  for(i = ilower; i < iupper; ++i)
  {
    j[0] = i % N;
    j[1] = (i + 100) % N;
    ierr = MatSetValues(A,1,&i,2,&j[0],&v[0],INSERT_VALUES);CHKERRQ(ierr);
	
	herr = HYPRE_IJMatrixSetValues(C, 1, &nnz, &i, (const int*)&j[0], &v[0]);
	if(herr) printf("Error in HYPRE_IJMatrixSetValues!");
  }
  
  ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  
  herr = HYPRE_IJMatrixAssemble(C);
  if(herr) printf("Error in HYPRE_IJMatrixAssemble!");
  
  herr = HYPRE_IJMatrixGetObject(C, (void **) &par_C);
  if(herr) printf("Error in HYPRE_IJMatrixGetObject!");
  herr = hypre_ParCSRMatrixTranspose(par_C, &par_CT, data); 
  if(herr) printf("Error in hypre_ParCSRMatrixTranspose!");
  
  ierr = MatHYPRE_IJMatrixCreate(A,&B);CHKERRQ(ierr);
  ierr = MatHYPRE_IJMatrixCopy(A,B);CHKERRQ(ierr);
  
  herr = HYPRE_IJMatrixGetObject(B, (void **) &par_B);
  if(herr) printf("Error in HYPRE_IJMatrixGetObject!");
  herr = hypre_ParCSRMatrixTranspose(par_B, &par_BT, data); 
  if(herr) printf("Error in hypre_ParCSRMatrixTranspose!");
      
  ierr = MatDestroy(&A);CHKERRQ(ierr);
  HYPRE_IJMatrixDestroy(B);
  HYPRE_IJMatrixDestroy(C);
  ierr = PetscFinalize();
  return 0;
}