Actual source code: ex1.c

petsc-3.13.4 2020-08-01
Report Typos and Errors

  2: static char help[] = "Tests solving linear system on 0 by 0 matrix, and KSPLSQR convergence test handling.\n\n";

  4:  #include <petscksp.h>

  6: static PetscErrorCode GetConvergenceTestName(PetscErrorCode (*converged)(KSP,PetscInt,PetscReal,KSPConvergedReason*,void*),char name[],size_t n)
  7: {

 11:   if (converged == KSPConvergedDefault) {
 12:     PetscStrncpy(name,"default",n);
 13:   } else if (converged == KSPConvergedSkip) {
 14:     PetscStrncpy(name,"skip",n);
 15:   } else if (converged == KSPLSQRConvergedDefault) {
 16:     PetscStrncpy(name,"lsqr",n);
 17:   } else {
 18:     PetscStrncpy(name,"other",n);
 19:   }
 20:   return(0);
 21: }

 23: int main(int argc,char **args)
 24: {
 25:   Mat            C;
 27:   PetscInt       N = 0;
 28:   Vec            u,b,x;
 29:   KSP            ksp;
 30:   PetscReal      norm;
 31:   PetscBool      flg=PETSC_FALSE;

 33:   PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;

 35:   /* create stiffness matrix */
 36:   MatCreate(PETSC_COMM_WORLD,&C);
 37:   MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,N,N);
 38:   MatSetFromOptions(C);
 39:   MatSetUp(C);
 40:   MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
 41:   MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);

 43:   /* create right hand side and solution */
 44:   VecCreate(PETSC_COMM_WORLD,&u);
 45:   VecSetSizes(u,PETSC_DECIDE,N);
 46:   VecSetFromOptions(u);
 47:   VecDuplicate(u,&b);
 48:   VecDuplicate(u,&x);
 49:   VecSet(u,0.0);
 50:   VecSet(b,0.0);

 52:   VecAssemblyBegin(b);
 53:   VecAssemblyEnd(b);

 55:   /* solve linear system */
 56:   KSPCreate(PETSC_COMM_WORLD,&ksp);
 57:   KSPSetOperators(ksp,C,C);
 58:   KSPSetFromOptions(ksp);
 59:   KSPSolve(ksp,b,u);

 61:   /* test proper handling of convergence test by KSPLSQR */
 62:   PetscOptionsGetBool(NULL,NULL,"-test_lsqr",&flg,NULL);
 63:   if (flg) {
 64:     char                  *type;
 65:     char                  convtestname[16];
 66:     PetscBool             islsqr;
 67:     PetscErrorCode        (*converged)(KSP,PetscInt,PetscReal,KSPConvergedReason*,void*);
 68:     PetscErrorCode        (*converged1)(KSP,PetscInt,PetscReal,KSPConvergedReason*,void*);
 69:     PetscErrorCode        (*destroy)(void*),(*destroy1)(void*);
 70:     void                  *ctx,*ctx1;

 72:     {
 73:       const char *typeP;
 74:       KSPGetType(ksp,&typeP);
 75:       PetscStrallocpy(typeP,&type);
 76:     }
 77:     PetscStrcmp(type,KSPLSQR,&islsqr);
 78:     KSPGetConvergenceTest(ksp,&converged,&ctx,&destroy);
 79:     GetConvergenceTestName(converged,convtestname,16);
 80:     PetscPrintf(PETSC_COMM_WORLD,"convergence test: %s\n",convtestname);
 81:     KSPSetType(ksp,KSPLSQR);
 82:     KSPGetConvergenceTest(ksp,&converged1,&ctx1,&destroy1);
 83:     if (converged1 != KSPLSQRConvergedDefault)  SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"convergence test should be KSPLSQRConvergedDefault");
 84:     if (destroy1 != KSPConvergedDefaultDestroy) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"convergence test destroy function should be KSPConvergedDefaultDestroy");
 85:     if (islsqr) {
 86:       if (converged1 != converged) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"convergence test should be kept");
 87:       if (destroy1 != destroy)     SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"convergence test destroy function should be kept");
 88:       if (ctx1 != ctx)             SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"convergence test context should be kept");
 89:     }
 90:     GetConvergenceTestName(converged1,convtestname,16);
 91:     KSPViewFromOptions(ksp,NULL,"-ksp1_view");
 92:     PetscPrintf(PETSC_COMM_WORLD,"convergence test: %s\n",convtestname);
 93:     KSPSetType(ksp,type);
 94:     KSPGetConvergenceTest(ksp,&converged1,&ctx1,&destroy1);
 95:     if (converged1 != converged) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"convergence test not reverted properly");
 96:     if (destroy1 != destroy)     SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"convergence test destroy function not reverted properly");
 97:     if (ctx1 != ctx)             SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"convergence test context not reverted properly");
 98:     GetConvergenceTestName(converged1,convtestname,16);
 99:     KSPViewFromOptions(ksp,NULL,"-ksp2_view");
100:     PetscPrintf(PETSC_COMM_WORLD,"convergence test: %s\n",convtestname);
101:     PetscFree(type);
102:   }

104:   MatMult(C,u,x);
105:   VecAXPY(x,-1.0,b);
106:   VecNorm(x,NORM_2,&norm);

108:   KSPDestroy(&ksp);
109:   VecDestroy(&u);
110:   VecDestroy(&x);
111:   VecDestroy(&b);
112:   MatDestroy(&C);
113:   PetscFinalize();
114:   return ierr;
115: }

117: /*TEST

119:     test:
120:       args:  -pc_type jacobi -ksp_monitor_short -ksp_gmres_cgs_refinement_type refine_always 

122:     test:
123:       suffix: 2
124:       nsize: 2
125:       args: -pc_type jacobi -ksp_monitor_short -ksp_gmres_cgs_refinement_type refine_always 

127:     test:
128:       suffix: 3
129:       args: -pc_type sor -pc_sor_symmetric -ksp_monitor_short -ksp_gmres_cgs_refinement_type refine_always

131:     test:
132:       suffix: 5
133:       args: -pc_type eisenstat -ksp_monitor_short -ksp_gmres_cgs_refinement_type refine_always

135:     testset:
136:       args: -test_lsqr -ksp{,1,2}_view -pc_type jacobi
137:       filter: grep -E "(^  type:|preconditioning|norm type|convergence test:)"
138:       test:
139:         suffix: lsqr_0
140:         args: -ksp_convergence_test {{default skip}separate output}
141:       test:
142:         suffix: lsqr_1
143:         args: -ksp_type cg -ksp_convergence_test {{default skip}separate output}
144:       test:
145:         suffix: lsqr_2
146:         args: -ksp_type lsqr


149: TEST*/