Actual source code: ex225.c
petsc-3.14.6 2021-03-30
2: static char help[] = "Test Hypre matrix APIs\n";
4: #include <petscmathypre.h>
6: int main(int argc,char **args)
7: {
8: Mat A, B, C;
9: PetscReal err;
10: PetscInt i,j,M = 20;
12: PetscMPIInt NP;
13: MPI_Comm comm;
14: PetscInt *rows;
16: PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
17: comm = PETSC_COMM_WORLD;
18: MPI_Comm_size(comm,&NP);
19: PetscOptionsGetInt(NULL,NULL,"-M",&M,NULL);
20: if (M < 6) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Matrix has to have more than 6 columns");
21: /* Hypre matrix */
22: MatCreate(comm,&B);
23: MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,M,M);
24: MatSetType(B,MATHYPRE);
25: MatHYPRESetPreallocation(B,9,NULL,9,NULL);
27: /* PETSc AIJ matrix */
28: MatCreate(comm,&A);
29: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,M,M);
30: MatSetType(A,MATAIJ);
31: MatSeqAIJSetPreallocation(A,9,NULL);
32: MatMPIAIJSetPreallocation(A,9,NULL,9,NULL);
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;
42: MatSetValues(B,1,&i,6,cols,vals,ADD_VALUES);
43: MatSetValues(B,1,&i,1,&i,value,ADD_VALUES);
44: MatSetValues(A,1,&i,6,cols,vals,ADD_VALUES);
45: MatSetValues(A,1,&i,1,&i,value,ADD_VALUES);
46: }
48: /* MAT_FLUSH_ASSEMBLY currently not supported */
49: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
50: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
51: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
52: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
54: /* Compare A and B */
55: MatConvert(B,MATAIJ,MAT_INITIAL_MATRIX,&C);
56: MatAXPY(C,-1.,A,SAME_NONZERO_PATTERN);
57: MatNorm(C,NORM_INFINITY,&err);
58: if (err > PETSC_SMALL) SETERRQ1(PetscObjectComm((PetscObject)B),PETSC_ERR_PLIB,"Error MatSetValues %g",err);
59: MatDestroy(&C);
61: /* MatZeroRows */
62: PetscMalloc1(M, &rows);
63: for (i=0; i<M; i++) rows[i] = i;
64: MatZeroRows(B, M, rows, 10.0, NULL, NULL);
65: MatSetOption(A, MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE);
66: MatZeroRows(A, M, rows, 10.0,NULL, NULL);
67: MatConvert(B,MATAIJ,MAT_INITIAL_MATRIX,&C);
68: MatAXPY(C,-1.,A,SAME_NONZERO_PATTERN);
69: MatNorm(C,NORM_INFINITY,&err);
70: if (err > PETSC_SMALL) SETERRQ1(PetscObjectComm((PetscObject)B),PETSC_ERR_PLIB,"Error MatZeroRows %g",err);
71: MatDestroy(&C);
72: PetscFree(rows);
74: /* Test MatZeroEntries */
75: MatZeroEntries(B);
76: MatConvert(B,MATAIJ,MAT_INITIAL_MATRIX,&C);
77: MatNorm(C,NORM_INFINITY,&err);
78: if (err > PETSC_SMALL) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Error MatZeroEntries %g",err);
79: MatDestroy(&C);
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};
87: for (j=0; j<6; j++)
88: vals[j] = ((PetscReal)j)/NP;
90: MatSetValues(B,1,&i,6,cols,vals,INSERT_VALUES);
91: MatSetValues(B,1,&i,1,&i,value,INSERT_VALUES);
92: MatSetValues(A,1,&i,6,cols,vals,INSERT_VALUES);
93: MatSetValues(A,1,&i,1,&i,value,INSERT_VALUES);
94: }
96: /* MAT_FLUSH_ASSEMBLY currently not supported */
97: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
98: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
99: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
100: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
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;
113: MatGetOwnershipRange(A,&rstart,&rend);
114: for (i=rstart; i<rend; i++) {
115: MatGetRow(A,i,&nzA,&idxA,&vA);
116: MatGetRow(B,i,&nzB,&idxB,&vB);
117: if (nzA!=nzB) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error MatGetRow %D", nzA-nzB);
118: PetscSortIntWithScalarArray(nzB,(PetscInt*)idxB,(PetscScalar*)vB);
119: PetscArraycmp(idxA,idxB,nzA,&flg);
120: if (!flg) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error MatGetRow %D (indices)",i);
121: PetscArraycmp(vA,vB,nzA,&flg);
122: if (!flg) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error MatGetRow %D (values)",i);
123: MatRestoreRow(A,i,&nzA,&idxA,&vA);
124: MatRestoreRow(B,i,&nzB,&idxB,&vB);
125: }
127: MatGetOwnershipRange(A,&rstart,&rend);
128: PetscCalloc3((rend-rstart)*6,&valuesA,(rend-rstart)*6,&valuesB,rend-rstart,&rows);
129: for (i=rstart; i<rend; i++) rows[i-rstart] =i;
131: MatGetValues(A,rend-rstart,rows,6,cols,valuesA);
132: MatGetValues(B,rend-rstart,rows,6,cols,valuesB);
134: for (i=0; i<(rend-rstart); i++) {
135: PetscArraycmp(valuesA + 6*i,valuesB + 6*i,6,&flg);
136: if (!flg) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error MatGetValues %D",i + rstart);
137: }
138: PetscFree3(valuesA,valuesB,rows);
139: }
141: /* Compare A and B */
142: MatConvert(B,MATAIJ,MAT_INITIAL_MATRIX,&C);
143: MatAXPY(C,-1.,A,SAME_NONZERO_PATTERN);
144: MatNorm(C,NORM_INFINITY,&err);
145: if (err > PETSC_SMALL) SETERRQ1(PetscObjectComm((PetscObject)B),PETSC_ERR_PLIB,"Error MatSetValues with INSERT_VALUES %g",err);
147: MatDestroy(&A);
148: MatDestroy(&B);
149: MatDestroy(&C);
151: PetscFinalize();
152: return ierr;
153: }
156: /*TEST
158: build:
159: requires: hypre
161: test:
162: suffix: 1
164: test:
165: output_file: output/ex225_1.out
166: nsize: 2
167: suffix: 2
169: TEST*/