1 2 static char help[] = "Test Hypre matrix APIs\n"; 3 4 #include <petscmathypre.h> 5 6 int main(int argc,char **args) 7 { 8 Mat A, B, C; 9 PetscReal err; 10 PetscInt i,j,M = 20; 11 PetscErrorCode ierr; 12 PetscMPIInt NP; 13 MPI_Comm comm; 14 PetscInt *rows; 15 16 ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr; 17 comm = PETSC_COMM_WORLD; 18 CHKERRMPI(MPI_Comm_size(comm,&NP)); 19 CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-M",&M,NULL)); 20 PetscCheckFalse(M < 6,PETSC_COMM_WORLD,PETSC_ERR_SUP,"Matrix has to have more than 6 columns"); 21 /* Hypre matrix */ 22 CHKERRQ(MatCreate(comm,&B)); 23 CHKERRQ(MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,M,M)); 24 CHKERRQ(MatSetType(B,MATHYPRE)); 25 CHKERRQ(MatHYPRESetPreallocation(B,9,NULL,9,NULL)); 26 27 /* PETSc AIJ matrix */ 28 CHKERRQ(MatCreate(comm,&A)); 29 CHKERRQ(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,M,M)); 30 CHKERRQ(MatSetType(A,MATAIJ)); 31 CHKERRQ(MatSeqAIJSetPreallocation(A,9,NULL)); 32 CHKERRQ(MatMPIAIJSetPreallocation(A,9,NULL,9,NULL)); 33 34 /*Set Values */ 35 for (i=0; i<M; i++) { 36 PetscInt cols[] = {0,1,2,3,4,5}; 37 PetscScalar vals[6] = {0}; 38 PetscScalar value[] = {100}; 39 for (j=0; j<6; j++) 40 vals[j] = ((PetscReal)j)/NP; 41 42 CHKERRQ(MatSetValues(B,1,&i,6,cols,vals,ADD_VALUES)); 43 CHKERRQ(MatSetValues(B,1,&i,1,&i,value,ADD_VALUES)); 44 CHKERRQ(MatSetValues(A,1,&i,6,cols,vals,ADD_VALUES)); 45 CHKERRQ(MatSetValues(A,1,&i,1,&i,value,ADD_VALUES)); 46 } 47 48 /* MAT_FLUSH_ASSEMBLY currently not supported */ 49 CHKERRQ(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 50 CHKERRQ(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 51 CHKERRQ(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 52 CHKERRQ(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 53 54 /* Compare A and B */ 55 CHKERRQ(MatConvert(B,MATAIJ,MAT_INITIAL_MATRIX,&C)); 56 CHKERRQ(MatAXPY(C,-1.,A,SAME_NONZERO_PATTERN)); 57 CHKERRQ(MatNorm(C,NORM_INFINITY,&err)); 58 PetscCheckFalse(err > PETSC_SMALL,PetscObjectComm((PetscObject)B),PETSC_ERR_PLIB,"Error MatSetValues %g",err); 59 CHKERRQ(MatDestroy(&C)); 60 61 /* MatZeroRows */ 62 CHKERRQ(PetscMalloc1(M, &rows)); 63 for (i=0; i<M; i++) rows[i] = i; 64 CHKERRQ(MatZeroRows(B, M, rows, 10.0, NULL, NULL)); 65 CHKERRQ(MatSetOption(A, MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE)); 66 CHKERRQ(MatZeroRows(A, M, rows, 10.0,NULL, NULL)); 67 CHKERRQ(MatConvert(B,MATAIJ,MAT_INITIAL_MATRIX,&C)); 68 CHKERRQ(MatAXPY(C,-1.,A,SAME_NONZERO_PATTERN)); 69 CHKERRQ(MatNorm(C,NORM_INFINITY,&err)); 70 PetscCheckFalse(err > PETSC_SMALL,PetscObjectComm((PetscObject)B),PETSC_ERR_PLIB,"Error MatZeroRows %g",err); 71 CHKERRQ(MatDestroy(&C)); 72 CHKERRQ(PetscFree(rows)); 73 74 /* Test MatZeroEntries */ 75 CHKERRQ(MatZeroEntries(B)); 76 CHKERRQ(MatConvert(B,MATAIJ,MAT_INITIAL_MATRIX,&C)); 77 CHKERRQ(MatNorm(C,NORM_INFINITY,&err)); 78 PetscCheckFalse(err > PETSC_SMALL,PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Error MatZeroEntries %g",err); 79 CHKERRQ(MatDestroy(&C)); 80 81 /* Insert Values */ 82 for (i=0; i<M; i++) { 83 PetscInt cols[] = {0,1,2,3,4,5}; 84 PetscScalar vals[6] = {0}; 85 PetscScalar value[] = {100}; 86 87 for (j=0; j<6; j++) 88 vals[j] = ((PetscReal)j)/NP; 89 90 CHKERRQ(MatSetValues(B,1,&i,6,cols,vals,INSERT_VALUES)); 91 CHKERRQ(MatSetValues(B,1,&i,1,&i,value,INSERT_VALUES)); 92 CHKERRQ(MatSetValues(A,1,&i,6,cols,vals,INSERT_VALUES)); 93 CHKERRQ(MatSetValues(A,1,&i,1,&i,value,INSERT_VALUES)); 94 } 95 96 /* MAT_FLUSH_ASSEMBLY currently not supported */ 97 CHKERRQ(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 98 CHKERRQ(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 99 CHKERRQ(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 100 CHKERRQ(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 101 102 /* Rows are not sorted with HYPRE so we need an intermediate sort 103 They use a temporary buffer, so we can sort inplace the const memory */ 104 { 105 const PetscInt *idxA,*idxB; 106 const PetscScalar *vA, *vB; 107 PetscInt rstart, rend, nzA, nzB; 108 PetscInt cols[] = {0,1,2,3,4,-5}; 109 PetscInt *rows; 110 PetscScalar *valuesA, *valuesB; 111 PetscBool flg; 112 113 CHKERRQ(MatGetOwnershipRange(A,&rstart,&rend)); 114 for (i=rstart; i<rend; i++) { 115 CHKERRQ(MatGetRow(A,i,&nzA,&idxA,&vA)); 116 CHKERRQ(MatGetRow(B,i,&nzB,&idxB,&vB)); 117 PetscCheckFalse(nzA!=nzB,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error MatGetRow %" PetscInt_FMT, nzA-nzB); 118 CHKERRQ(PetscSortIntWithScalarArray(nzB,(PetscInt*)idxB,(PetscScalar*)vB)); 119 CHKERRQ(PetscArraycmp(idxA,idxB,nzA,&flg)); 120 PetscCheckFalse(!flg,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error MatGetRow %" PetscInt_FMT " (indices)",i); 121 CHKERRQ(PetscArraycmp(vA,vB,nzA,&flg)); 122 PetscCheckFalse(!flg,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error MatGetRow %" PetscInt_FMT " (values)",i); 123 CHKERRQ(MatRestoreRow(A,i,&nzA,&idxA,&vA)); 124 CHKERRQ(MatRestoreRow(B,i,&nzB,&idxB,&vB)); 125 } 126 127 CHKERRQ(MatGetOwnershipRange(A,&rstart,&rend)); 128 CHKERRQ(PetscCalloc3((rend-rstart)*6,&valuesA,(rend-rstart)*6,&valuesB,rend-rstart,&rows)); 129 for (i=rstart; i<rend; i++) rows[i-rstart] =i; 130 131 CHKERRQ(MatGetValues(A,rend-rstart,rows,6,cols,valuesA)); 132 CHKERRQ(MatGetValues(B,rend-rstart,rows,6,cols,valuesB)); 133 134 for (i=0; i<(rend-rstart); i++) { 135 CHKERRQ(PetscArraycmp(valuesA + 6*i,valuesB + 6*i,6,&flg)); 136 PetscCheckFalse(!flg,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error MatGetValues %" PetscInt_FMT,i + rstart); 137 } 138 CHKERRQ(PetscFree3(valuesA,valuesB,rows)); 139 } 140 141 /* Compare A and B */ 142 CHKERRQ(MatConvert(B,MATAIJ,MAT_INITIAL_MATRIX,&C)); 143 CHKERRQ(MatAXPY(C,-1.,A,SAME_NONZERO_PATTERN)); 144 CHKERRQ(MatNorm(C,NORM_INFINITY,&err)); 145 PetscCheckFalse(err > PETSC_SMALL,PetscObjectComm((PetscObject)B),PETSC_ERR_PLIB,"Error MatSetValues with INSERT_VALUES %g",err); 146 147 CHKERRQ(MatDestroy(&A)); 148 CHKERRQ(MatDestroy(&B)); 149 CHKERRQ(MatDestroy(&C)); 150 151 ierr = PetscFinalize(); 152 return ierr; 153 } 154 155 /*TEST 156 157 build: 158 requires: hypre 159 160 test: 161 suffix: 1 162 requires: !defined(PETSC_HAVE_HYPRE_DEVICE) 163 164 test: 165 suffix: 2 166 requires: !defined(PETSC_HAVE_HYPRE_DEVICE) 167 output_file: output/ex225_1.out 168 nsize: 2 169 170 TEST*/ 171