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;
 11:   PetscErrorCode     ierr;
 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:     for (j=0; j<6; j++)
 39:       vals[j] = ((PetscReal)j)/NP;
 40:     PetscScalar value[] = {100};

 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:   PetscCalloc1(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:     for (j=0; j<6; j++)
 86:       vals[j] = ((PetscReal)j)/NP;

 88:     PetscScalar value[] = {100};

 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:   {
103:     const PetscInt *idxA,*idxB;
104:     const PetscScalar *vA, *vB;
105:     PetscInt rstart, rend, nzA, nzB, j;
106:     err = 0.0;
107:     MatGetOwnershipRange(A,&rstart,&rend);
108:     for (i=rstart; i<rend; i++) {
109:       MatGetRow(A,i,&nzA, &idxA,&vA);
110:       MatGetRow(B,i,&nzB, &idxB,&vB);
111:       if (nzA!=nzB) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error MatGetRow %d", nzA-nzB);
112:       for (j=0;j<nzA;j++) {
113:         err +=idxA[j]-idxB[j];
114:         err +=vA[j]-vA[j];
115:       }
116:       MatRestoreRow(A,i,&nzA, &idxA,&vA);
117:       MatRestoreRow(B,i,&nzB, &idxB,&vB);
118:     }
119:     if (err > PETSC_SMALL) SETERRQ1(PetscObjectComm((PetscObject)B),PETSC_ERR_PLIB,"Error MatGetRow %g",err);

121:     PetscInt    cols[] = {0,1,2,3,4,-5};
122:     PetscInt    *rows;
123:     PetscScalar *valuesA, *valuesB;
124:     MatGetOwnershipRange(A,&rstart,&rend);
125:     PetscCalloc3((rend-rstart)*6,&valuesA,(rend-rstart)*6,&valuesB,rend-rstart,&rows);
126:     for (i=rstart; i<rend; i++)
127:       rows[i-rstart] =i;

129:     MatGetValues(A,rend-rstart,rows,6,cols,valuesA);
130:     MatGetValues(B,rend-rstart,rows,6,cols,valuesB);

132:     err = 0;
133:     for (i=0; i<(rend-rstart); i++)
134:       for (j=0; j<6; j++)
135:         err += valuesA[i*6+j] - valuesB[i*6+j];

137:     PetscFree3(valuesA,valuesB,rows);

139:     if (err > PETSC_SMALL) SETERRQ1(PetscObjectComm((PetscObject)B),PETSC_ERR_PLIB,"Error MatGetValueS %g",err);
140:   }

142:   /* Compare A and B */
143:   MatConvert(B,MATAIJ,MAT_INITIAL_MATRIX,&C);
144:   MatAXPY(C,-1.,A,SAME_NONZERO_PATTERN);
145:   MatNorm(C,NORM_INFINITY,&err);
146:   if (err > PETSC_SMALL) SETERRQ1(PetscObjectComm((PetscObject)B),PETSC_ERR_PLIB,"Error MatSetValues with INSERT_VALUES %g",err);

148:   MatDestroy(&A);
149:   MatDestroy(&B);
150:   MatDestroy(&C);

152:   PetscFinalize();
153:   return ierr;
154: }


157: /*TEST

159:    build:
160:       requires: hypre

162:    test:
163:       suffix: 1
164:       requires: hypre

166: TEST*/