1 #include <petsc.h> 2 3 static char help[] = "Solves a series of linear systems using KSPHPDDM.\n\n"; 4 5 int main(int argc,char **args) 6 { 7 Vec x,b; /* computed solution and RHS */ 8 Mat A; /* linear system matrix */ 9 KSP ksp; /* linear solver context */ 10 #if defined(PETSC_HAVE_HPDDM) 11 Mat U; /* deflation space */ 12 PetscBool flg; 13 #endif 14 PetscInt i,j,nmat = 10; 15 PetscViewer viewer; 16 PetscBool reset = PETSC_FALSE; 17 char name[256]; 18 char dir[PETSC_MAX_PATH_LEN]; 19 PetscErrorCode ierr; 20 21 ierr = PetscInitialize(&argc,&args,NULL,help);if (ierr) return ierr; 22 ierr = PetscStrcpy(dir,".");CHKERRQ(ierr); 23 ierr = PetscOptionsGetString(NULL,NULL,"-load_dir",dir,sizeof(dir),NULL);CHKERRQ(ierr); 24 ierr = PetscOptionsGetInt(NULL,NULL,"-nmat",&nmat,NULL);CHKERRQ(ierr); 25 ierr = PetscOptionsGetBool(NULL,NULL,"-reset",&reset,NULL);CHKERRQ(ierr); 26 ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr); 27 ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr); 28 ierr = KSPSetOperators(ksp,A,A);CHKERRQ(ierr); 29 for (i=0; i<nmat; i++) { 30 j = i+400; 31 ierr = PetscSNPrintf(name,sizeof(name),"%s/A_%d.dat",dir,j);CHKERRQ(ierr); 32 ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,name,FILE_MODE_READ,&viewer);CHKERRQ(ierr); 33 ierr = MatLoad(A,viewer);CHKERRQ(ierr); 34 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 35 if (i == 0) { 36 ierr = MatCreateVecs(A,&x,&b);CHKERRQ(ierr); 37 } 38 ierr = PetscSNPrintf(name,sizeof(name),"%s/rhs_%d.dat",dir,j);CHKERRQ(ierr); 39 ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,name,FILE_MODE_READ,&viewer);CHKERRQ(ierr); 40 ierr = VecLoad(b,viewer);CHKERRQ(ierr); 41 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 42 ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr); 43 ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr); 44 #if defined(PETSC_HAVE_HPDDM) 45 ierr = PetscObjectTypeCompare((PetscObject)ksp,KSPHPDDM,&flg);CHKERRQ(ierr); 46 if (flg && reset) { 47 ierr = KSPHPDDMGetDeflationSpace(ksp,&U);CHKERRQ(ierr); 48 ierr = KSPReset(ksp);CHKERRQ(ierr); 49 ierr = KSPSetOperators(ksp,A,A);CHKERRQ(ierr); 50 ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr); 51 ierr = KSPSetUp(ksp);CHKERRQ(ierr); 52 if (U) { 53 ierr = KSPHPDDMSetDeflationSpace(ksp,U);CHKERRQ(ierr); 54 ierr = MatDestroy(&U);CHKERRQ(ierr); 55 } 56 } 57 #endif 58 } 59 ierr = VecDestroy(&x);CHKERRQ(ierr); 60 ierr = VecDestroy(&b);CHKERRQ(ierr); 61 ierr = MatDestroy(&A);CHKERRQ(ierr); 62 ierr = KSPDestroy(&ksp);CHKERRQ(ierr); 63 ierr = PetscFinalize(); 64 return ierr; 65 } 66 67 /*TEST 68 69 test: 70 suffix: 1 71 nsize: 1 72 requires: hpddm datafilespath double !complex !define(PETSC_USE_64BIT_INDICES) 73 args: -nmat 1 -pc_type none -ksp_converged_reason -ksp_type {{gmres hpddm}shared ouput} -ksp_max_it 1000 -ksp_gmres_restart 1000 -ksp_rtol 1e-10 -ksp_hpddm_type {{gmres bgmres}shared output} -options_left no -load_dir ${DATAFILESPATH}/matrices/hpddm/GCRODR 74 75 test: 76 requires: hpddm datafilespath double !complex !define(PETSC_USE_64BIT_INDICES) 77 suffix: 1_icc 78 nsize: 1 79 args: -nmat 1 -pc_type icc -ksp_converged_reason -ksp_type {{gmres hpddm}shared ouput} -ksp_max_it 1000 -ksp_gmres_restart 1000 -ksp_rtol 1e-10 -load_dir ${DATAFILESPATH}/matrices/hpddm/GCRODR 80 81 testset: 82 requires: hpddm datafilespath double !complex !define(PETSC_USE_64BIT_INDICES) 83 args: -nmat 3 -pc_type none -ksp_converged_reason -ksp_type hpddm -ksp_max_it 1000 -ksp_gmres_restart 40 -ksp_rtol 1e-10 -ksp_hpddm_type {{gcrodr bgcrodr}shared output} -ksp_hpddm_recycle 20 -load_dir ${DATAFILESPATH}/matrices/hpddm/GCRODR 84 test: 85 nsize: 1 86 suffix: 2_seq 87 output_file: output/ex75_2.out 88 test: 89 nsize: 2 90 suffix: 2_par 91 output_file: output/ex75_2.out 92 93 test: 94 requires: hpddm datafilespath double !complex !define(PETSC_USE_64BIT_INDICES) 95 suffix: 2_icc 96 nsize: 1 97 args: -nmat 3 -pc_type icc -ksp_converged_reason -ksp_type hpddm -ksp_max_it 1000 -ksp_gmres_restart 40 -ksp_rtol 1e-10 -ksp_hpddm_type gcrodr -ksp_hpddm_recycle 20 -reset {{false true}shared output} -load_dir ${DATAFILESPATH}/matrices/hpddm/GCRODR 98 99 TEST*/ 100