Actual source code: ex77.c

petsc-master 2020-07-05
Report Typos and Errors
  1:  #include <petsc.h>

  3: static char help[] = "Solves a linear system with a block of right-hand sides using KSPHPDDM.\n\n";

  5: int main(int argc,char **args)
  6: {
  7:   Mat                X,B;         /* computed solutions and RHS */
  8:   Vec                cx,cb;       /* columns of X and B */
  9:   Mat                A,KA = NULL; /* linear system matrix */
 10:   KSP                ksp;         /* linear solver context */
 11:   PC                 pc;          /* preconditioner context */
 12:   Mat                F;           /* factored matrix from the preconditioner context */
 13:   const PetscScalar  *b;
 14:   PetscScalar        *x,*S = NULL,*T = NULL;
 15:   PetscReal          norm;
 16:   PetscInt           m,M,N = 5,i,j;
 17:   const char         *deft = MATAIJ;
 18:   PetscViewer        viewer;
 19:   char               dir[PETSC_MAX_PATH_LEN],name[256],type[256];
 20:   PetscBool          breakdown = PETSC_FALSE,flg;
 21:   KSPConvergedReason reason;
 22:   PetscErrorCode     ierr;

 24:   PetscInitialize(&argc,&args,NULL,help);if (ierr) return ierr;
 25:   PetscStrcpy(dir,".");
 26:   PetscOptionsGetString(NULL,NULL,"-load_dir",dir,sizeof(dir),NULL);
 27:   PetscOptionsGetInt(NULL,NULL,"-N",&N,NULL);
 28:   PetscOptionsGetBool(NULL,NULL,"-breakdown",&breakdown,NULL);
 29:   MatCreate(PETSC_COMM_WORLD,&A);
 30:   KSPCreate(PETSC_COMM_WORLD,&ksp);
 31:   KSPSetOperators(ksp,A,A);
 32:   PetscSNPrintf(name,sizeof(name),"%s/A_400.dat",dir);
 33:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,name,FILE_MODE_READ,&viewer);
 34:   MatLoad(A,viewer);
 35:   PetscViewerDestroy(&viewer);
 36:   PetscOptionsBegin(PETSC_COMM_WORLD,"","","");
 37:   PetscOptionsFList("-mat_type","Matrix type","MatSetType",MatList,deft,type,256,&flg);
 38:   PetscOptionsEnd();
 39:   if (flg) {
 40:     PetscStrcmp(type,MATKAIJ,&flg);
 41:     if (!flg) {
 42:       MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);
 43:       MatConvert(A,type,MAT_INPLACE_MATRIX,&A);
 44:     }
 45:     else {
 46:       PetscCalloc2(N*N,&S,N*N,&T);
 47:       for (i=0; i<N; i++) {
 48:         for (j=0; j<N; j++) {
 49:           S[i*(N+1)] = 1e+6; /* really easy problem used for testing */
 50:           T[i*(N+1)] = 1e-2;
 51:         }
 52:       }
 53:       MatCreateKAIJ(A,N,N,S,T,&KA);
 54:     }
 55:   }
 56:   MatGetLocalSize(A,&m,NULL);
 57:   MatCreateDense(PETSC_COMM_WORLD,m,PETSC_DECIDE,PETSC_DECIDE,N,NULL,&B);
 58:   MatCreateDense(PETSC_COMM_WORLD,m,PETSC_DECIDE,PETSC_DECIDE,N,NULL,&X);
 59:   if (!breakdown) {
 60:     MatSetRandom(B,NULL);
 61:   }
 62:   KSPSetFromOptions(ksp);
 63:   if (!flg) {
 64:     if (!breakdown) {
 65:       KSPMatSolve(ksp,B,X);
 66:       KSPGetMatSolveBlockSize(ksp,&M);
 67:       if (M != PETSC_DECIDE) {
 68:         KSPSetMatSolveBlockSize(ksp,PETSC_DECIDE);
 69:         MatZeroEntries(X);
 70:         KSPMatSolve(ksp,B,X);
 71:       }
 72:       KSPGetPC(ksp,&pc);
 73:       PetscObjectTypeCompare((PetscObject)pc,PCLU,&flg);
 74:       if (flg) {
 75:         PCFactorGetMatrix(pc,&F);
 76:         MatMatSolve(F,B,B);
 77:         MatAYPX(B,-1.0,X,SAME_NONZERO_PATTERN);
 78:         MatNorm(B,NORM_INFINITY,&norm);
 79:         if (norm > 100*PETSC_MACHINE_EPSILON) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_PLIB,"KSPMatSolve() and MatMatSolve() difference has nonzero norm %g",(double)norm);
 80:       }
 81:     } else {
 82:       MatZeroEntries(B);
 83:       KSPMatSolve(ksp,B,X);
 84:       KSPGetConvergedReason(ksp,&reason);
 85:       if (reason != KSP_CONVERGED_HAPPY_BREAKDOWN) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_PLIB,"KSPConvergedReason() %s != KSP_CONVERGED_HAPPY_BREAKDOWN",KSPConvergedReasons[reason]);
 86:       MatDenseGetArrayWrite(B,&x);
 87:       for (i=0; i<m*N; ++i) x[i] = 1.0;
 88:       MatDenseRestoreArrayWrite(B,&x);
 89:       KSPMatSolve(ksp,B,X);
 90:       KSPGetConvergedReason(ksp,&reason);
 91:       if (reason != KSP_DIVERGED_BREAKDOWN) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_PLIB,"KSPConvergedReason() %s != KSP_DIVERGED_BREAKDOWN",KSPConvergedReasons[reason]);
 92:     }
 93:   } else {
 94:     KSPSetOperators(ksp,KA,KA);
 95:     MatGetSize(KA,&M,NULL);
 96:     /* from column- to row-major to be consistent with MatKAIJ format */
 97:     MatTranspose(B,MAT_INPLACE_MATRIX,&B);
 98:     MatDenseGetArrayRead(B,&b);
 99:     MatDenseGetArray(X,&x);
100:     VecCreateMPIWithArray(PETSC_COMM_WORLD,1,m*N,M,b,&cb);
101:     VecCreateMPIWithArray(PETSC_COMM_WORLD,1,m*N,M,x,&cx);
102:     KSPSolve(ksp,cb,cx);
103:     VecDestroy(&cx);
104:     VecDestroy(&cb);
105:     MatDenseRestoreArray(X,&x);
106:     MatDenseRestoreArrayRead(B,&b);
107:   }
108:   MatDestroy(&X);
109:   MatDestroy(&B);
110:   PetscFree2(S,T);
111:   MatDestroy(&KA);
112:   MatDestroy(&A);
113:   KSPDestroy(&ksp);
114:   PetscFinalize();
115:   return ierr;
116: }

118: /*TEST

120:    testset:
121:       nsize: 2
122:       requires: datafilespath double !complex !define(PETSC_USE_64BIT_INDICES)
123:       args: -ksp_converged_reason -ksp_max_it 500 -load_dir ${DATAFILESPATH}/matrices/hpddm/GCRODR -mat_type {{aij sbaij}shared output}
124:       test:
125:          suffix: 1
126:          args:
127:       test:
128:          suffix: 2
129:          requires: hpddm
130:          args: -ksp_type hpddm -pc_type asm -ksp_hpddm_type {{gmres bgmres}separate output}
131:       test:
132:          suffix: 3
133:          requires: hpddm
134:          args: -ksp_type hpddm -ksp_hpddm_recycle 5 -ksp_hpddm_type {{gcrodr bgcrodr}separate output}
135:       test:
136:          nsize: 4
137:          suffix: 4
138:          requires: hpddm
139:          args: -ksp_rtol 1e-4 -ksp_type hpddm -ksp_hpddm_recycle 5 -ksp_hpddm_type bgcrodr -ksp_view_final_residual -N 12 -ksp_matsolve_block_size 5

141:    test:
142:       nsize: 1
143:       requires: hpddm datafilespath double !complex !define(PETSC_USE_64BIT_INDICES)
144:       suffix: preonly
145:       args: -N 6 -load_dir ${DATAFILESPATH}/matrices/hpddm/GCRODR -pc_type lu -ksp_type hpddm -ksp_hpddm_type preonly

147:    # to avoid breakdown failures, use -ksp_hpddm_deflation_tol, cf. KSPHPDDM documentation
148:    test:
149:       nsize: 1
150:       requires: hpddm datafilespath double !complex !define(PETSC_USE_64BIT_INDICES)
151:       suffix: breakdown
152:       output_file: output/ex77_preonly.out
153:       args: -N 3 -load_dir ${DATAFILESPATH}/matrices/hpddm/GCRODR -pc_type none -ksp_type hpddm -ksp_hpddm_type {{bcg bgmres bgcrodr bfbcg}shared output} -breakdown

155:    test:
156:       nsize: 2
157:       requires: hpddm datafilespath double !complex !define(PETSC_USE_64BIT_INDICES)
158:       args: -N 12 -ksp_converged_reason -ksp_max_it 500 -load_dir ${DATAFILESPATH}/matrices/hpddm/GCRODR -mat_type kaij -pc_type pbjacobi -ksp_type hpddm -ksp_hpddm_type {{gmres bgmres}separate output}

160:    test:
161:       nsize: 4
162:       requires: hpddm datafilespath double !complex !define(PETSC_USE_64BIT_INDICES) slepc
163:       suffix: 4_slepc
164:       output_file: output/ex77_4.out
165:       filter: sed "/^ksp_hpddm_recycle_ Linear eigensolve converged/d"
166:       args: -ksp_converged_reason -ksp_max_it 500 -load_dir ${DATAFILESPATH}/matrices/hpddm/GCRODR -ksp_rtol 1e-4 -ksp_type hpddm -ksp_hpddm_recycle 5 -ksp_hpddm_type bgcrodr -ksp_view_final_residual -N 12 -ksp_matsolve_block_size 5 -ksp_hpddm_recycle_redistribute 2 -ksp_hpddm_recycle_mat_type {{aij dense}shared output} -ksp_hpddm_recycle_eps_converged_reason -ksp_hpddm_recycle_st_pc_type redundant

168:    test:
169:       nsize: 4
170:       requires: hpddm datafilespath double !complex !define(PETSC_USE_64BIT_INDICES) slepc elemental
171:       suffix: 4_elemental
172:       output_file: output/ex77_4.out
173:       filter: sed "/^ksp_hpddm_recycle_ Linear eigensolve converged/d"
174:       args: -ksp_converged_reason -ksp_max_it 500 -load_dir ${DATAFILESPATH}/matrices/hpddm/GCRODR -ksp_rtol 1e-4 -ksp_type hpddm -ksp_hpddm_recycle 5 -ksp_hpddm_type bgcrodr -ksp_view_final_residual -N 12 -ksp_matsolve_block_size 5 -ksp_hpddm_recycle_redistribute 2 -ksp_hpddm_recycle_mat_type elemental -ksp_hpddm_recycle_eps_converged_reason

176: TEST*/