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:   const PetscScalar *b;
 12:   PetscScalar       *x,*S = NULL,*T = NULL;
 13:   PetscInt          m,n,M,N = 5,i,j;
 14:   const char        *deft = MATAIJ;
 15:   PetscViewer       viewer;
 16:   char              dir[PETSC_MAX_PATH_LEN],name[256],type[256];
 17:   PetscBool         flg;
 18:   PetscErrorCode    ierr;

 20:   PetscInitialize(&argc,&args,NULL,help);if (ierr) return ierr;
 21:   PetscStrcpy(dir,".");
 22:   PetscOptionsGetString(NULL,NULL,"-load_dir",dir,sizeof(dir),NULL);
 23:   PetscOptionsGetInt(NULL,NULL,"-N",&N,NULL);
 24:   MatCreate(PETSC_COMM_WORLD,&A);
 25:   KSPCreate(PETSC_COMM_WORLD,&ksp);
 26:   KSPSetOperators(ksp,A,A);
 27:   PetscSNPrintf(name,sizeof(name),"%s/A_400.dat",dir);
 28:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,name,FILE_MODE_READ,&viewer);
 29:   MatLoad(A,viewer);
 30:   PetscViewerDestroy(&viewer);
 31:   PetscOptionsBegin(PETSC_COMM_WORLD,"","","");
 32:   PetscOptionsFList("-mat_type","Matrix type","MatSetType",MatList,deft,type,256,&flg);
 33:   PetscOptionsEnd();
 34:   if (flg) {
 35:     PetscStrcmp(type,MATKAIJ,&flg);
 36:     if (!flg) {
 37:       MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);
 38:       MatConvert(A,type,MAT_INPLACE_MATRIX,&A);
 39:     }
 40:     else {
 41:       PetscCalloc2(N*N,&S,N*N,&T);
 42:       for (i=0; i<N; i++) {
 43:         for (j=0; j<N; j++) {
 44:           S[i*(N+1)] = 1e+6; /* really easy problem used for testing */
 45:           T[i*(N+1)] = 1e-2;
 46:         }
 47:       }
 48:       MatCreateKAIJ(A,N,N,S,T,&KA);
 49:     }
 50:   }
 51:   MatGetLocalSize(A,&m,NULL);
 52:   MatCreateDense(PETSC_COMM_WORLD,m,PETSC_DECIDE,PETSC_DECIDE,N,NULL,&B);
 53:   MatCreateDense(PETSC_COMM_WORLD,m,PETSC_DECIDE,PETSC_DECIDE,N,NULL,&X);
 54:   MatSetRandom(B,NULL);
 55:   KSPSetFromOptions(ksp);
 56:   if (!flg) {
 57:     KSPSetUp(ksp);
 58:     PetscObjectTypeCompare((PetscObject)ksp,KSPHPDDM,&flg);
 59: #if defined(PETSC_HAVE_HPDDM)
 60:     if (flg) {
 61:       KSPHPDDMMatSolve(ksp,B,X);
 62:     } else
 63: #endif
 64:     {
 65:       MatGetSize(A,&M,NULL);
 66:       for (n=0; n<N; n++) {
 67:         MatDenseGetColumn(B,n,(PetscScalar**)&b);
 68:         MatDenseGetColumn(X,n,&x);
 69:         VecCreateMPIWithArray(PETSC_COMM_WORLD,1,m,M,b,&cb);
 70:         VecCreateMPIWithArray(PETSC_COMM_WORLD,1,m,M,x,&cx);
 71:         KSPSolve(ksp,cb,cx);
 72:         VecDestroy(&cx);
 73:         VecDestroy(&cb);
 74:         MatDenseRestoreColumn(X,&x);
 75:         MatDenseRestoreColumn(B,(PetscScalar**)&b);
 76:       }
 77:     }
 78:   } else {
 79:     KSPSetOperators(ksp,KA,KA);
 80:     MatGetSize(KA,&M,NULL);
 81:     /* from column- to row-major to be consistent with MatKAIJ format */
 82:     MatTranspose(B,MAT_INPLACE_MATRIX,&B);
 83:     MatDenseGetArrayRead(B,&b);
 84:     MatDenseGetArray(X,&x);
 85:     VecCreateMPIWithArray(PETSC_COMM_WORLD,1,m*N,M,b,&cb);
 86:     VecCreateMPIWithArray(PETSC_COMM_WORLD,1,m*N,M,x,&cx);
 87:     KSPSolve(ksp,cb,cx);
 88:     VecDestroy(&cx);
 89:     VecDestroy(&cb);
 90:     MatDenseRestoreArray(X,&x);
 91:     MatDenseRestoreArrayRead(B,&b);
 92:   }
 93:   MatDestroy(&X);
 94:   MatDestroy(&B);
 95:   PetscFree2(S,T);
 96:   MatDestroy(&KA);
 97:   MatDestroy(&A);
 98:   KSPDestroy(&ksp);
 99:   PetscFinalize();
100:   return ierr;
101: }

103: /*TEST

105:    testset:
106:       nsize: 2
107:       requires: hpddm datafilespath double !complex !define(PETSC_USE_64BIT_INDICES)
108:       args: -ksp_converged_reason -ksp_max_it 500 -load_dir ${DATAFILESPATH}/matrices/hpddm/GCRODR -mat_type {{aij sbaij}shared output}
109:       test:
110:          suffix: 1
111:          args:
112:       test:
113:          suffix: 2
114:          args: -ksp_type hpddm -pc_type asm -ksp_hpddm_type {{gmres bgmres}separate_output}
115:       test:
116:          suffix: 3
117:          args: -ksp_type hpddm -ksp_hpddm_recycle 10 -ksp_hpddm_type {{gcrodr bgcrodr}separate_output}

119:    test:
120:       nsize: 2
121:       requires: hpddm datafilespath double !complex !define(PETSC_USE_64BIT_INDICES)
122:       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}

124: TEST*/