Actual source code: ex88.c

  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: };

 13: static PetscErrorCode MatMult_User(Mat A,Vec X,Vec Y)
 14: {
 15:   User user;

 19:   MatShellGetContext(A,&user);
 20:   MatMult(user->B,X,Y);
 21:   return(0);
 22: }

 26: static PetscErrorCode MatMultTranspose_User(Mat A,Vec X,Vec Y)
 27: {
 28:   User user;

 32:   MatShellGetContext(A,&user);
 33:   MatMultTranspose(user->B,X,Y);
 34:   return(0);
 35: }

 39: static PetscErrorCode MatGetDiagonal_User(Mat A,Vec X)
 40: {
 41:   User user;

 45:   MatShellGetContext(A,&user);
 46:   MatGetDiagonal(user->B,X);
 47:   return(0);
 48: }

 52: static PetscErrorCode TestMatrix(Mat A,Vec X,Vec Y,Vec Z)
 53: {
 55:   Vec            W1,W2;
 56:   Mat            E;
 57:   const char     *mattypename;
 58:   PetscViewer    viewer = PETSC_VIEWER_STDOUT_WORLD;

 61:   VecDuplicate(X,&W1);
 62:   VecDuplicate(X,&W2);
 63:   MatScale(A,31);
 64:   MatShift(A,37);
 65:   MatDiagonalScale(A,X,Y);
 66:   MatScale(A,41);
 67:   MatDiagonalScale(A,Y,Z);
 68:   MatComputeExplicitOperator(A,&E);

 70:   PetscObjectGetType((PetscObject)A,&mattypename);
 71:   PetscViewerASCIIPrintf(viewer,"Matrix of type: %s\n",mattypename);
 72:   MatView(E,viewer);
 73:   MatMult(A,Z,W1);
 74:   MatMultTranspose(A,W1,W2);
 75:   VecView(W2,viewer);
 76:   MatGetDiagonal(A,W2);
 77:   VecView(W2,viewer);
 78:   MatDestroy(&E);
 79:   VecDestroy(&W1);
 80:   VecDestroy(&W2);
 81:   return(0);
 82: }

 86: int main(int argc,char **args)
 87: {
 88:   const PetscScalar xvals[] = {11,13},yvals[] = {17,19},zvals[] = {23,29};
 89:   const PetscInt inds[] = {0,1};
 90:   PetscScalar avals[] = {2,3,5,7};
 91:   Mat            A,S,D[4],N;
 92:   Vec            X,Y,Z;
 93:   User           user;
 94:   PetscInt       i;

 97:   PetscInitialize(&argc,&args,(char *)0,help);
 98:   MatCreateSeqAIJ(PETSC_COMM_WORLD,2,2,2,PETSC_NULL,&A);
 99:   MatSetValues(A,2,inds,2,inds,avals,INSERT_VALUES);
100:   VecCreateSeq(PETSC_COMM_WORLD,2,&X);
101:   VecDuplicate(X,&Y);
102:   VecDuplicate(X,&Z);
103:   VecSetValues(X,2,inds,xvals,INSERT_VALUES);
104:   VecSetValues(Y,2,inds,yvals,INSERT_VALUES);
105:   VecSetValues(Z,2,inds,zvals,INSERT_VALUES);
106:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
107:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
108:   VecAssemblyBegin(X);
109:   VecAssemblyBegin(Y);
110:   VecAssemblyBegin(Z);
111:   VecAssemblyEnd(X);
112:   VecAssemblyEnd(Y);
113:   VecAssemblyEnd(Z);

115:   PetscNew(struct _n_User,&user);
116:   user->B = A;

118:   MatCreateShell(PETSC_COMM_WORLD,2,2,2,2,user,&S);
119:   MatShellSetOperation(S,MATOP_MULT,(void(*)(void))MatMult_User);
120:   MatShellSetOperation(S,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMultTranspose_User);
121:   MatShellSetOperation(S,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_User);

123:   for (i=0; i<4; i++) {
124:     MatCreateSeqDense(PETSC_COMM_WORLD,1,1,&avals[i],&D[i]);
125:   }
126:   MatCreateNest(PETSC_COMM_WORLD,2,PETSC_NULL,2,PETSC_NULL,D,&N);

128:   TestMatrix(S,X,Y,Z);
129:   TestMatrix(A,X,Y,Z);
130:   TestMatrix(N,X,Y,Z);

132:   for (i=0; i<4; i++) {MatDestroy(&D[i]);}
133:   MatDestroy(&A);
134:   MatDestroy(&S);
135:   MatDestroy(&N);
136:   VecDestroy(&X);
137:   VecDestroy(&Y);
138:   VecDestroy(&Z);
139:   PetscFree(user);
140:   PetscFinalize();
141:   return 0;
142: }