2: static char help[] = "Tests MatShift(), MatScale(), and MatDiagonalScale() for SHELL and NEST matrices\n\n";

  4:  #include <petscmat.h>

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

 11: static PetscErrorCode MatView_User(Mat A,PetscViewer viewer)
 12: {
 13:   User           user;

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

 22: static PetscErrorCode MatMult_User(Mat A,Vec X,Vec Y)
 23: {
 24:   User           user;

 28:   MatShellGetContext(A,&user);
 29:   MatMult(user->B,X,Y);
 30:   return(0);
 31: }

 33: static PetscErrorCode MatMultTranspose_User(Mat A,Vec X,Vec Y)
 34: {
 35:   User           user;

 39:   MatShellGetContext(A,&user);
 40:   MatMultTranspose(user->B,X,Y);
 41:   return(0);
 42: }

 44: static PetscErrorCode MatGetDiagonal_User(Mat A,Vec X)
 45: {
 46:   User           user;

 50:   MatShellGetContext(A,&user);
 51:   MatGetDiagonal(user->B,X);
 52:   return(0);
 53: }

 55: static PetscErrorCode MatDiagonalSet_User(Mat A,Vec D,InsertMode is)
 56: {
 57:   User           user;

 61:   MatShellGetContext(A,&user);
 62:   MatDiagonalSet(user->B,D,is);
 63:   return(0);
 64: }


 67: static PetscErrorCode TestMatrix(Mat A,Vec X,Vec Y,Vec Z)
 68: {
 70:   Vec            W1,W2;
 71:   Mat            E;
 72:   const char     *mattypename;
 73:   PetscViewer    viewer = PETSC_VIEWER_STDOUT_WORLD;

 76:   VecDuplicate(X,&W1);
 77:   VecDuplicate(X,&W2);
 78:   MatScale(A,31);
 79:   MatShift(A,37);
 80:   MatDiagonalScale(A,X,Y);
 81:   MatScale(A,41);
 82:   MatDiagonalScale(A,Y,Z);
 83:   MatComputeExplicitOperator(A,&E);

 85:   PetscObjectGetType((PetscObject)A,&mattypename);
 86:   PetscViewerASCIIPrintf(viewer,"Matrix of type: %s\n",mattypename);
 87:   MatView(E,viewer);
 88:   MatMult(A,Z,W1);
 89:   MatMultTranspose(A,W1,W2);
 90:   VecView(W2,viewer);
 91:   MatGetDiagonal(A,W2);
 92:   VecView(W2,viewer);
 93:   /* MATSHELL does not support MatDiagonalSet after MatScale */
 94:   if (strncmp(mattypename, "shell", 5)) {
 95:     MatDiagonalSet(A,X,INSERT_VALUES);
 96:     MatGetDiagonal(A,W1);
 97:     VecView(W1,viewer);
 98:   } else {
 99:     PetscViewerASCIIPrintf(viewer,"MatDiagonalSet not tested on MATSHELL\n");
100:   }
101:   MatDestroy(&E);
102:   VecDestroy(&W1);
103:   VecDestroy(&W2);
104:   return(0);
105: }

107: int main(int argc,char **args)
108: {
109:   const PetscScalar xvals[] = {11,13},yvals[] = {17,19},zvals[] = {23,29};
110:   const PetscInt    inds[]  = {0,1};
111:   PetscScalar       avals[] = {2,3,5,7};
112:   Mat               A,S,D[4],N;
113:   Vec               X,Y,Z;
114:   User              user;
115:   PetscInt          i;
116:   PetscErrorCode    ierr;

118:   PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
119:   MatCreateSeqAIJ(PETSC_COMM_WORLD,2,2,2,NULL,&A);
120:   MatSetUp(A);
121:   MatSetValues(A,2,inds,2,inds,avals,INSERT_VALUES);
122:   VecCreateSeq(PETSC_COMM_WORLD,2,&X);
123:   VecDuplicate(X,&Y);
124:   VecDuplicate(X,&Z);
125:   VecSetValues(X,2,inds,xvals,INSERT_VALUES);
126:   VecSetValues(Y,2,inds,yvals,INSERT_VALUES);
127:   VecSetValues(Z,2,inds,zvals,INSERT_VALUES);
128:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
129:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
130:   VecAssemblyBegin(X);
131:   VecAssemblyBegin(Y);
132:   VecAssemblyBegin(Z);
133:   VecAssemblyEnd(X);
134:   VecAssemblyEnd(Y);
135:   VecAssemblyEnd(Z);

137:   PetscNew(&user);
138:   user->B = A;

140:   MatCreateShell(PETSC_COMM_WORLD,2,2,2,2,user,&S);
141:   MatSetUp(S);
142:   MatShellSetOperation(S,MATOP_VIEW,(void (*)(void))MatView_User);
143:   MatShellSetOperation(S,MATOP_MULT,(void (*)(void))MatMult_User);
144:   MatShellSetOperation(S,MATOP_MULT_TRANSPOSE,(void (*)(void))MatMultTranspose_User);
145:   MatShellSetOperation(S,MATOP_GET_DIAGONAL,(void (*)(void))MatGetDiagonal_User);
146:   MatShellSetOperation(S,MATOP_DIAGONAL_SET,(void (*)(void))MatDiagonalSet_User);

148:   for (i=0; i<4; i++) {
149:     MatCreateSeqDense(PETSC_COMM_WORLD,1,1,&avals[i],&D[i]);
150:   }
151:   MatCreateNest(PETSC_COMM_WORLD,2,NULL,2,NULL,D,&N);
152:   MatSetUp(N);

154:   TestMatrix(S,X,Y,Z);
155:   TestMatrix(A,X,Y,Z);
156:   TestMatrix(N,X,Y,Z);

158:   for (i=0; i<4; i++) {MatDestroy(&D[i]);}
159:   MatDestroy(&A);
160:   MatDestroy(&S);
161:   MatDestroy(&N);
162:   VecDestroy(&X);
163:   VecDestroy(&Y);
164:   VecDestroy(&Z);
165:   PetscFree(user);
166:   PetscFinalize();
167:   return ierr;
168: }