2: static char help[] = "Tests incorrect use of MatDiagonalSet() for SHELL matrices\n\n";

  4:  #include <petscmat.h>

  6: typedef struct _n_User *User;
  7: struct _n_User {
  8:   Mat B;
  9: };

 11: static PetscErrorCode MatGetDiagonal_User(Mat A,Vec X)
 12: {
 13:   User           user;

 17:   MatShellGetContext(A,&user);
 18:   MatGetDiagonal(user->B,X);
 19:   return(0);
 20: }

 22: static PetscErrorCode MatDiagonalSet_User(Mat A,Vec D,InsertMode is)
 23: {
 24:   User           user;

 28:   MatShellGetContext(A,&user);
 29:   MatDiagonalSet(user->B,D,is);
 30:   return(0);
 31: }

 33: PetscErrorCode ErrorHandler(MPI_Comm comm,int line,const char *func,const char *file,PetscErrorCode n,PetscErrorType p,const char *mess,void *ctx)
 34: {
 36:   PetscPrintf(PETSC_COMM_SELF,"[ERROR] %s: %s\n",func,mess);
 37:   return(0);
 38: }

 40: int main(int argc,char **args)
 41: {
 42:   const PetscScalar xvals[] = {11,13};
 43:   const PetscInt    inds[]  = {0,1};
 44:   PetscScalar       avals[] = {2,3,5,7};
 45:   Mat               A,S;
 46:   Vec               X,Y;
 47:   User              user;
 48:   PetscBool         flag;
 49:   PetscErrorCode    ierr;

 51:   PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
 52:   MatCreateSeqAIJ(PETSC_COMM_WORLD,2,2,2,NULL,&A);
 53:   MatSetUp(A);
 54:   MatSetValues(A,2,inds,2,inds,avals,INSERT_VALUES);
 55:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
 56:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
 57:   VecCreateSeq(PETSC_COMM_WORLD,2,&X);
 58:   VecSetValues(X,2,inds,xvals,INSERT_VALUES);
 59:   VecAssemblyBegin(X);
 60:   VecAssemblyEnd(X);
 61:   VecDuplicate(X,&Y);

 63:   PetscNew(&user);
 64:   user->B = A;

 66:   MatCreateShell(PETSC_COMM_WORLD,2,2,2,2,user,&S);
 67:   MatSetUp(S);
 68:   MatShellSetOperation(S,MATOP_GET_DIAGONAL,(void (*)(void))MatGetDiagonal_User);
 69:   MatShellSetOperation(S,MATOP_DIAGONAL_SET,(void (*)(void))MatDiagonalSet_User);

 71:   MatShift(S,42);
 72:   MatGetDiagonal(S,Y);
 73:   VecView(Y,PETSC_VIEWER_STDOUT_WORLD);
 74:   MatDiagonalSet(S,X,INSERT_VALUES);
 75:   MatGetDiagonal(S,Y);
 76:   VecView(Y,PETSC_VIEWER_STDOUT_WORLD);
 77:   VecEqual(X,Y,&flag);
 78:   if (!flag) {
 79:     PetscPrintf(PETSC_COMM_WORLD,"Diagonal not set correctly!\n");
 80:   }

 82:   /* NOTE: This test case will fail */
 83:   PetscPushErrorHandler(ErrorHandler, NULL);
 84:   MatScale(S,42);
 85:   MatGetDiagonal(S,Y);
 86:   VecView(Y,PETSC_VIEWER_STDOUT_WORLD);
 87:   MatDiagonalSet(S,X,INSERT_VALUES);
 88:   MatGetDiagonal(S,Y);
 89:   VecView(Y,PETSC_VIEWER_STDOUT_WORLD);
 90:   PetscPopErrorHandler();

 92:   MatDestroy(&A);
 93:   MatDestroy(&S);
 94:   VecDestroy(&X);
 95:   VecDestroy(&Y);
 96:   PetscFree(user);
 97:   PetscFinalize();
 98:   return ierr;
 99: }