Actual source code: test1.c

      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 MatMult_Shell(Mat S,Vec x,Vec y);
     16: static PetscErrorCode MatMultTranspose_Shell(Mat S,Vec x,Vec y);
     17: #if defined(PETSC_USE_COMPLEX)
     18: static PetscErrorCode MatMultHermitianTranspose_Shell(Mat S,Vec x,Vec y);
     19: #endif
     20: static PetscErrorCode MatGetDiagonal_Shell(Mat S,Vec diag);
     21: static PetscErrorCode MatDuplicate_Shell(Mat S,MatDuplicateOption op,Mat *M);

     23: static PetscErrorCode MyShellMatCreate(Mat *A,Mat *M)
     24: {
     25:   MPI_Comm       comm;
     26:   PetscInt       n;

     28:   PetscFunctionBeginUser;
     29:   PetscCall(MatGetSize(*A,&n,NULL));
     30:   PetscCall(PetscObjectGetComm((PetscObject)*A,&comm));
     31:   PetscCall(MatCreateShell(comm,PETSC_DECIDE,PETSC_DECIDE,n,n,A,M));
     32:   PetscCall(MatShellSetOperation(*M,MATOP_MULT,(PetscErrorCodeFn*)MatMult_Shell));
     33:   PetscCall(MatShellSetOperation(*M,MATOP_MULT_TRANSPOSE,(PetscErrorCodeFn*)MatMultTranspose_Shell));
     34: #if defined(PETSC_USE_COMPLEX)
     35:   PetscCall(MatShellSetOperation(*M,MATOP_MULT_HERMITIAN_TRANSPOSE,(PetscErrorCodeFn*)MatMultHermitianTranspose_Shell));
     36: #endif
     37:   PetscCall(MatShellSetOperation(*M,MATOP_GET_DIAGONAL,(PetscErrorCodeFn*)MatGetDiagonal_Shell));
     38:   PetscCall(MatShellSetOperation(*M,MATOP_DUPLICATE,(PetscErrorCodeFn*)MatDuplicate_Shell));
     39:   PetscFunctionReturn(PETSC_SUCCESS);
     40: }

     42: int main(int argc,char **argv)
     43: {
     44:   Mat            A,S,mat[1];
     45:   ST             st;
     46:   Vec            v,w;
     47:   STType         type;
     48:   KSP            ksp;
     49:   PC             pc;
     50:   PetscScalar    sigma;
     51:   PetscInt       n=10,i,Istart,Iend;

     53:   PetscFunctionBeginUser;
     54:   PetscCall(SlepcInitialize(&argc,&argv,NULL,help));
     55:   PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
     56:   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n1-D Laplacian with shell matrices, n=%" PetscInt_FMT "\n\n",n));

     58:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     59:      Compute the operator matrix for the 1-D Laplacian
     60:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

     62:   PetscCall(MatCreate(PETSC_COMM_WORLD,&A));
     63:   PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n));
     64:   PetscCall(MatSetFromOptions(A));

     66:   PetscCall(MatGetOwnershipRange(A,&Istart,&Iend));
     67:   for (i=Istart;i<Iend;i++) {
     68:     if (i>0) PetscCall(MatSetValue(A,i,i-1,-1.0,INSERT_VALUES));
     69:     if (i<n-1) PetscCall(MatSetValue(A,i,i+1,-1.0,INSERT_VALUES));
     70:     PetscCall(MatSetValue(A,i,i,2.0,INSERT_VALUES));
     71:   }
     72:   PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY));
     73:   PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY));

     75:   /* create the shell version of A */
     76:   PetscCall(MyShellMatCreate(&A,&S));

     78:   /* work vectors */
     79:   PetscCall(MatCreateVecs(A,&v,&w));
     80:   PetscCall(VecSet(v,1.0));

     82:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     83:                 Create the spectral transformation object
     84:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

     86:   PetscCall(STCreate(PETSC_COMM_WORLD,&st));
     87:   mat[0] = S;
     88:   PetscCall(STSetMatrices(st,1,mat));
     89:   PetscCall(STSetTransform(st,PETSC_TRUE));
     90:   PetscCall(STSetFromOptions(st));

     92:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     93:               Apply the transformed operator for several ST's
     94:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

     96:   /* shift, sigma=0.0 */
     97:   PetscCall(STSetUp(st));
     98:   PetscCall(STGetType(st,&type));
     99:   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"ST type %s\n",type));
    100:   PetscCall(STApply(st,v,w));
    101:   PetscCall(VecView(w,NULL));
    102:   PetscCall(STApplyHermitianTranspose(st,v,w));
    103:   PetscCall(VecView(w,NULL));

    105:   /* shift, sigma=0.1 */
    106:   sigma = 0.1;
    107:   PetscCall(STSetShift(st,sigma));
    108:   PetscCall(STGetShift(st,&sigma));
    109:   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"With shift=%g\n",(double)PetscRealPart(sigma)));
    110:   PetscCall(STApply(st,v,w));
    111:   PetscCall(VecView(w,NULL));

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

    127:   /* sinvert, sigma=-0.5 */
    128:   sigma = -0.5;
    129:   PetscCall(STSetShift(st,sigma));
    130:   PetscCall(STGetShift(st,&sigma));
    131:   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"With shift=%g\n",(double)PetscRealPart(sigma)));
    132:   PetscCall(STApply(st,v,w));
    133:   PetscCall(VecView(w,NULL));

    135:   PetscCall(STDestroy(&st));
    136:   PetscCall(MatDestroy(&A));
    137:   PetscCall(MatDestroy(&S));
    138:   PetscCall(VecDestroy(&v));
    139:   PetscCall(VecDestroy(&w));
    140:   PetscCall(SlepcFinalize());
    141:   return 0;
    142: }

    144: static PetscErrorCode MatMult_Shell(Mat S,Vec x,Vec y)
    145: {
    146:   Mat               *A;

    148:   PetscFunctionBeginUser;
    149:   PetscCall(MatShellGetContext(S,&A));
    150:   PetscCall(MatMult(*A,x,y));
    151:   PetscFunctionReturn(PETSC_SUCCESS);
    152: }

    154: static PetscErrorCode MatMultTranspose_Shell(Mat S,Vec x,Vec y)
    155: {
    156:   Mat               *A;

    158:   PetscFunctionBeginUser;
    159:   PetscCall(MatShellGetContext(S,&A));
    160:   PetscCall(MatMultTranspose(*A,x,y));
    161:   PetscFunctionReturn(PETSC_SUCCESS);
    162: }

    164: #if defined(PETSC_USE_COMPLEX)
    165: static PetscErrorCode MatMultHermitianTranspose_Shell(Mat S,Vec x,Vec y)
    166: {
    167:   Mat               *A;

    169:   PetscFunctionBeginUser;
    170:   PetscCall(MatShellGetContext(S,&A));
    171:   PetscCall(MatMultHermitianTranspose(*A,x,y));
    172:   PetscFunctionReturn(PETSC_SUCCESS);
    173: }
    174: #endif

    176: static PetscErrorCode MatGetDiagonal_Shell(Mat S,Vec diag)
    177: {
    178:   Mat               *A;

    180:   PetscFunctionBeginUser;
    181:   PetscCall(MatShellGetContext(S,&A));
    182:   PetscCall(MatGetDiagonal(*A,diag));
    183:   PetscFunctionReturn(PETSC_SUCCESS);
    184: }

    186: static PetscErrorCode MatDuplicate_Shell(Mat S,MatDuplicateOption op,Mat *M)
    187: {
    188:   Mat            *A;

    190:   PetscFunctionBeginUser;
    191:   PetscCall(MatShellGetContext(S,&A));
    192:   PetscCall(MyShellMatCreate(A,M));
    193:   PetscFunctionReturn(PETSC_SUCCESS);
    194: }

    196: /*TEST

    198:    test:
    199:       suffix: 1
    200:       args: -st_matmode {{inplace shell}}
    201:       requires: !single

    203: TEST*/