Actual source code: test1.c

slepc-3.18.2 2023-01-26
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain

  6:    This file is part of SLEPc.
  7:    SLEPc is distributed under a 2-clause BSD license (see LICENSE).
  8:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  9: */

 11: static char help[] = "Test ST with shell matrices.\n\n";

 13: #include <slepcst.h>

 15: static PetscErrorCode MatGetDiagonal_Shell(Mat S,Vec diag);
 16: static PetscErrorCode MatMultTranspose_Shell(Mat S,Vec x,Vec y);
 17: static PetscErrorCode MatMult_Shell(Mat S,Vec x,Vec y);
 18: static PetscErrorCode MatDuplicate_Shell(Mat S,MatDuplicateOption op,Mat *M);

 20: static PetscErrorCode MyShellMatCreate(Mat *A,Mat *M)
 21: {
 22:   MPI_Comm       comm;
 23:   PetscInt       n;

 26:   MatGetSize(*A,&n,NULL);
 27:   PetscObjectGetComm((PetscObject)*A,&comm);
 28:   MatCreateShell(comm,PETSC_DECIDE,PETSC_DECIDE,n,n,A,M);
 29:   MatShellSetOperation(*M,MATOP_MULT,(void(*)(void))MatMult_Shell);
 30:   MatShellSetOperation(*M,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMultTranspose_Shell);
 31:   MatShellSetOperation(*M,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_Shell);
 32:   MatShellSetOperation(*M,MATOP_DUPLICATE,(void(*)(void))MatDuplicate_Shell);
 33:   return 0;
 34: }

 36: int main(int argc,char **argv)
 37: {
 38:   Mat            A,S,mat[1];
 39:   ST             st;
 40:   Vec            v,w;
 41:   STType         type;
 42:   KSP            ksp;
 43:   PC             pc;
 44:   PetscScalar    sigma;
 45:   PetscInt       n=10,i,Istart,Iend;

 48:   SlepcInitialize(&argc,&argv,(char*)0,help);
 49:   PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
 50:   PetscPrintf(PETSC_COMM_WORLD,"\n1-D Laplacian with shell matrices, n=%" PetscInt_FMT "\n\n",n);

 52:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 53:      Compute the operator matrix for the 1-D Laplacian
 54:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 56:   MatCreate(PETSC_COMM_WORLD,&A);
 57:   MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n);
 58:   MatSetFromOptions(A);
 59:   MatSetUp(A);

 61:   MatGetOwnershipRange(A,&Istart,&Iend);
 62:   for (i=Istart;i<Iend;i++) {
 63:     if (i>0) MatSetValue(A,i,i-1,-1.0,INSERT_VALUES);
 64:     if (i<n-1) MatSetValue(A,i,i+1,-1.0,INSERT_VALUES);
 65:     MatSetValue(A,i,i,2.0,INSERT_VALUES);
 66:   }
 67:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
 68:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

 70:   /* create the shell version of A */
 71:   MyShellMatCreate(&A,&S);

 73:   /* work vectors */
 74:   MatCreateVecs(A,&v,&w);
 75:   VecSet(v,1.0);

 77:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 78:                 Create the spectral transformation object
 79:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 81:   STCreate(PETSC_COMM_WORLD,&st);
 82:   mat[0] = S;
 83:   STSetMatrices(st,1,mat);
 84:   STSetTransform(st,PETSC_TRUE);
 85:   STSetFromOptions(st);

 87:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 88:               Apply the transformed operator for several ST's
 89:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 91:   /* shift, sigma=0.0 */
 92:   STSetUp(st);
 93:   STGetType(st,&type);
 94:   PetscPrintf(PETSC_COMM_WORLD,"ST type %s\n",type);
 95:   STApply(st,v,w);
 96:   VecView(w,NULL);
 97:   STApplyTranspose(st,v,w);
 98:   VecView(w,NULL);

100:   /* shift, sigma=0.1 */
101:   sigma = 0.1;
102:   STSetShift(st,sigma);
103:   STGetShift(st,&sigma);
104:   PetscPrintf(PETSC_COMM_WORLD,"With shift=%g\n",(double)PetscRealPart(sigma));
105:   STApply(st,v,w);
106:   VecView(w,NULL);

108:   /* sinvert, sigma=0.1 */
109:   STPostSolve(st);   /* undo changes if inplace */
110:   STSetType(st,STSINVERT);
111:   STGetKSP(st,&ksp);
112:   KSPSetType(ksp,KSPGMRES);
113:   KSPGetPC(ksp,&pc);
114:   PCSetType(pc,PCJACOBI);
115:   STGetType(st,&type);
116:   PetscPrintf(PETSC_COMM_WORLD,"ST type %s\n",type);
117:   STGetShift(st,&sigma);
118:   PetscPrintf(PETSC_COMM_WORLD,"With shift=%g\n",(double)PetscRealPart(sigma));
119:   STApply(st,v,w);
120:   VecView(w,NULL);

122:   /* sinvert, sigma=-0.5 */
123:   sigma = -0.5;
124:   STSetShift(st,sigma);
125:   STGetShift(st,&sigma);
126:   PetscPrintf(PETSC_COMM_WORLD,"With shift=%g\n",(double)PetscRealPart(sigma));
127:   STApply(st,v,w);
128:   VecView(w,NULL);

130:   STDestroy(&st);
131:   MatDestroy(&A);
132:   MatDestroy(&S);
133:   VecDestroy(&v);
134:   VecDestroy(&w);
135:   SlepcFinalize();
136:   return 0;
137: }

139: static PetscErrorCode MatMult_Shell(Mat S,Vec x,Vec y)
140: {
141:   Mat               *A;

144:   MatShellGetContext(S,&A);
145:   MatMult(*A,x,y);
146:   return 0;
147: }

149: static PetscErrorCode MatMultTranspose_Shell(Mat S,Vec x,Vec y)
150: {
151:   Mat               *A;

154:   MatShellGetContext(S,&A);
155:   MatMultTranspose(*A,x,y);
156:   return 0;
157: }

159: static PetscErrorCode MatGetDiagonal_Shell(Mat S,Vec diag)
160: {
161:   Mat               *A;

164:   MatShellGetContext(S,&A);
165:   MatGetDiagonal(*A,diag);
166:   return 0;
167: }

169: static PetscErrorCode MatDuplicate_Shell(Mat S,MatDuplicateOption op,Mat *M)
170: {
171:   Mat            *A;

174:   MatShellGetContext(S,&A);
175:   MyShellMatCreate(A,M);
176:   return 0;
177: }

179: /*TEST

181:    test:
182:       suffix: 1
183:       args: -st_matmode {{inplace shell}}
184:       output_file: output/test1_1.out
185:       requires: !single

187: TEST*/