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 ierr = MPI_Comm_size(comm,&NP);CHKERRQ(ierr); 19 ierr = PetscOptionsGetInt(NULL,NULL,"-M",&M,NULL);CHKERRQ(ierr); 20 if (M < 6) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Matrix has to have more than 6 columns"); 21 /* Hypre matrix */ 22 ierr = MatCreate(comm,&B);CHKERRQ(ierr); 23 ierr = MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,M,M);CHKERRQ(ierr); 24 ierr = MatSetType(B,MATHYPRE);CHKERRQ(ierr); 25 ierr = MatHYPRESetPreallocation(B,9,NULL,9,NULL);CHKERRQ(ierr); 26 27 /* PETSc AIJ matrix */ 28 ierr = MatCreate(comm,&A);CHKERRQ(ierr); 29 ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,M,M);CHKERRQ(ierr); 30 ierr = MatSetType(A,MATAIJ);CHKERRQ(ierr); 31 ierr = MatSeqAIJSetPreallocation(A,9,NULL);CHKERRQ(ierr); 32 ierr = MatMPIAIJSetPreallocation(A,9,NULL,9,NULL);CHKERRQ(ierr); 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 ierr = MatSetValues(B,1,&i,6,cols,vals,ADD_VALUES);CHKERRQ(ierr); 43 ierr = MatSetValues(B,1,&i,1,&i,value,ADD_VALUES);CHKERRQ(ierr); 44 ierr = MatSetValues(A,1,&i,6,cols,vals,ADD_VALUES);CHKERRQ(ierr); 45 ierr = MatSetValues(A,1,&i,1,&i,value,ADD_VALUES);CHKERRQ(ierr); 46 } 47 48 /* MAT_FLUSH_ASSEMBLY currently not supported */ 49 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 50 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 51 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 52 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 53 54 /* Compare A and B */ 55 ierr = MatConvert(B,MATAIJ,MAT_INITIAL_MATRIX,&C);CHKERRQ(ierr); 56 ierr = MatAXPY(C,-1.,A,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 57 ierr = MatNorm(C,NORM_INFINITY,&err);CHKERRQ(ierr); 58 if (err > PETSC_SMALL) SETERRQ1(PetscObjectComm((PetscObject)B),PETSC_ERR_PLIB,"Error MatSetValues %g",err); 59 ierr = MatDestroy(&C);CHKERRQ(ierr); 60 61 /* MatZeroRows */ 62 ierr = PetscMalloc1(M, &rows);CHKERRQ(ierr); 63 for (i=0; i<M; i++) rows[i] = i; 64 ierr = MatZeroRows(B, M, rows, 10.0, NULL, NULL);CHKERRQ(ierr); 65 ierr = MatSetOption(A, MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE);CHKERRQ(ierr); 66 ierr = MatZeroRows(A, M, rows, 10.0,NULL, NULL);CHKERRQ(ierr); 67 ierr = MatConvert(B,MATAIJ,MAT_INITIAL_MATRIX,&C);CHKERRQ(ierr); 68 ierr = MatAXPY(C,-1.,A,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 69 ierr = MatNorm(C,NORM_INFINITY,&err);CHKERRQ(ierr); 70 if (err > PETSC_SMALL) SETERRQ1(PetscObjectComm((PetscObject)B),PETSC_ERR_PLIB,"Error MatZeroRows %g",err); 71 ierr = MatDestroy(&C);CHKERRQ(ierr); 72 ierr = PetscFree(rows);CHKERRQ(ierr); 73 74 /* Test MatZeroEntries */ 75 ierr = MatZeroEntries(B);CHKERRQ(ierr); 76 ierr = MatConvert(B,MATAIJ,MAT_INITIAL_MATRIX,&C);CHKERRQ(ierr); 77 ierr = MatNorm(C,NORM_INFINITY,&err);CHKERRQ(ierr); 78 if (err > PETSC_SMALL) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Error MatZeroEntries %g",err); 79 ierr = MatDestroy(&C);CHKERRQ(ierr); 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 ierr = MatSetValues(B,1,&i,6,cols,vals,INSERT_VALUES);CHKERRQ(ierr); 91 ierr = MatSetValues(B,1,&i,1,&i,value,INSERT_VALUES);CHKERRQ(ierr); 92 ierr = MatSetValues(A,1,&i,6,cols,vals,INSERT_VALUES);CHKERRQ(ierr); 93 ierr = MatSetValues(A,1,&i,1,&i,value,INSERT_VALUES);CHKERRQ(ierr); 94 } 95 96 /* MAT_FLUSH_ASSEMBLY currently not supported */ 97 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 98 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 99 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 100 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 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 ierr = MatGetOwnershipRange(A,&rstart,&rend);CHKERRQ(ierr); 114 for (i=rstart; i<rend; i++) { 115 ierr = MatGetRow(A,i,&nzA,&idxA,&vA);CHKERRQ(ierr); 116 ierr = MatGetRow(B,i,&nzB,&idxB,&vB);CHKERRQ(ierr); 117 if (nzA!=nzB) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error MatGetRow %D", nzA-nzB); 118 ierr = PetscSortIntWithScalarArray(nzB,(PetscInt*)idxB,(PetscScalar*)vB);CHKERRQ(ierr); 119 ierr = PetscArraycmp(idxA,idxB,nzA,&flg);CHKERRQ(ierr); 120 if (!flg) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error MatGetRow %D (indices)",i); 121 ierr = PetscArraycmp(vA,vB,nzA,&flg);CHKERRQ(ierr); 122 if (!flg) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error MatGetRow %D (values)",i); 123 ierr = MatRestoreRow(A,i,&nzA,&idxA,&vA);CHKERRQ(ierr); 124 ierr = MatRestoreRow(B,i,&nzB,&idxB,&vB);CHKERRQ(ierr); 125 } 126 127 ierr = MatGetOwnershipRange(A,&rstart,&rend);CHKERRQ(ierr); 128 ierr = PetscCalloc3((rend-rstart)*6,&valuesA,(rend-rstart)*6,&valuesB,rend-rstart,&rows);CHKERRQ(ierr); 129 for (i=rstart; i<rend; i++) rows[i-rstart] =i; 130 131 ierr = MatGetValues(A,rend-rstart,rows,6,cols,valuesA);CHKERRQ(ierr); 132 ierr = MatGetValues(B,rend-rstart,rows,6,cols,valuesB);CHKERRQ(ierr); 133 134 for (i=0; i<(rend-rstart); i++) { 135 ierr = PetscArraycmp(valuesA + 6*i,valuesB + 6*i,6,&flg);CHKERRQ(ierr); 136 if (!flg) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error MatGetValues %D",i + rstart); 137 } 138 ierr = PetscFree3(valuesA,valuesB,rows);CHKERRQ(ierr); 139 } 140 141 /* Compare A and B */ 142 ierr = MatConvert(B,MATAIJ,MAT_INITIAL_MATRIX,&C);CHKERRQ(ierr); 143 ierr = MatAXPY(C,-1.,A,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 144 ierr = MatNorm(C,NORM_INFINITY,&err);CHKERRQ(ierr); 145 if (err > PETSC_SMALL) SETERRQ1(PetscObjectComm((PetscObject)B),PETSC_ERR_PLIB,"Error MatSetValues with INSERT_VALUES %g",err); 146 147 ierr = MatDestroy(&A);CHKERRQ(ierr); 148 ierr = MatDestroy(&B);CHKERRQ(ierr); 149 ierr = MatDestroy(&C);CHKERRQ(ierr); 150 151 ierr = PetscFinalize(); 152 return ierr; 153 } 154 155 156 /*TEST 157 158 build: 159 requires: hypre 160 161 test: 162 suffix: 1 163 164 test: 165 output_file: output/ex225_1.out 166 nsize: 2 167 suffix: 2 168 169 TEST*/ 170