Actual source code: ex37.c
2: static char help[] = "Test MatGetMultiProcBlock() \n\
3: Reads a PETSc matrix and vector from a file and solves a linear system.\n\n";
4: /*
5: Example:
6: mpiexec -n 4 ./ex37 -f <input_file> -nsubcomm 2 -subcomm_view -subcomm_type <1 or 2>
7: */
9: #include <petscksp.h>
13: int main(int argc,char **args)
14: {
15: KSP subksp;
16: Mat A,subA;
17: Vec x,b,u,subb,subx,subu;
18: PetscViewer fd;
19: char file[PETSC_MAX_PATH_LEN];
20: PetscBool flg;
22: PetscInt i,m,n,its;
23: PetscReal norm;
24: PetscMPIInt rank,size;
25: MPI_Comm comm,subcomm;
26: PetscSubcomm psubcomm;
27: PetscInt nsubcomm=1,id;
28: PetscScalar *barray,*xarray,*uarray,*array,one=1.0;
29: PetscInt type=1;
31: PetscInitialize(&argc,&args,(char *)0,help);
32: /* Load the matrix */
33: PetscOptionsGetString(PETSC_NULL,"-f",file,PETSC_MAX_PATH_LEN,&flg);
34: if (!flg) SETERRQ(PETSC_COMM_WORLD,1,"Must indicate binary file with the -f option");
35: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);
36:
37: /* Load the matrix; then destroy the viewer.*/
38: MatCreate(PETSC_COMM_WORLD,&A);
39: MatLoad(A,fd);
40: PetscViewerDestroy(&fd);
42: PetscObjectGetComm((PetscObject)A,&comm);
43: MPI_Comm_size(comm,&size);
44: MPI_Comm_rank(comm,&rank);
45:
46: /* Create rhs vector b */
47: MatGetLocalSize(A,&m,PETSC_NULL);
48: VecCreate(PETSC_COMM_WORLD,&b);
49: VecSetSizes(b,m,PETSC_DECIDE);
50: VecSetFromOptions(b);
51: VecSet(b,one);
53: VecDuplicate(b,&x);
54: VecDuplicate(b,&u);
55: VecSet(x,0.0);
57: /* Test MatGetMultiProcBlock() */
58: PetscOptionsGetInt(PETSC_NULL,"-nsubcomm",&nsubcomm,PETSC_NULL);
59: PetscOptionsGetInt(PETSC_NULL,"-subcomm_type",&type,PETSC_NULL);
61: PetscSubcommCreate(comm,&psubcomm);
62: PetscSubcommSetNumber(psubcomm,nsubcomm);
63: if (type == PETSC_SUBCOMM_GENERAL){/* user provides color, subrank and duprank */
64: PetscMPIInt color,subrank,duprank,subsize;
65: duprank = size-1 - rank;
66: subsize = size/nsubcomm;
67: if (subsize*nsubcomm != size) SETERRQ2(comm,PETSC_ERR_SUP,"This example requires nsubcomm %D divides nproc %D",nsubcomm,size);
68: color = duprank/subsize;
69: subrank = duprank - color*subsize;
70: PetscSubcommSetTypeGeneral(psubcomm,color,subrank,duprank);
71: } else if (type == PETSC_SUBCOMM_CONTIGUOUS){
72: PetscSubcommSetType(psubcomm,PETSC_SUBCOMM_CONTIGUOUS);
73: } else if(type == PETSC_SUBCOMM_INTERLACED){
74: PetscSubcommSetType(psubcomm,PETSC_SUBCOMM_INTERLACED);
75: } else {
76: SETERRQ1(psubcomm->parent,PETSC_ERR_SUP,"PetscSubcommType %D is not supported yet",type);
77: }
78: subcomm = psubcomm->comm;
80: PetscOptionsHasName(PETSC_NULL, "-subcomm_view", &flg);
81: if (flg){
82: PetscMPIInt subsize,subrank,duprank;
83: MPI_Comm_size((MPI_Comm)subcomm,&subsize);
84: MPI_Comm_rank((MPI_Comm)subcomm,&subrank);
85: MPI_Comm_rank((MPI_Comm)psubcomm->dupparent,&duprank);
87: PetscSynchronizedPrintf(comm,"[%D], color %D, sub-size %D, sub-rank %D, duprank %D\n",rank,psubcomm->color,subsize,subrank,duprank);
88: PetscSynchronizedFlush(comm);
89: }
91: /* Create subA */
92: MatGetMultiProcBlock(A,subcomm,&subA);
94: /* Create sub vectors without arrays. Place b's and x's local arrays into subb and subx */
95: MatGetLocalSize(subA,&m,&n);
96: VecCreateMPIWithArray(subcomm,m,PETSC_DECIDE,PETSC_NULL,&subb);
97: VecCreateMPIWithArray(subcomm,n,PETSC_DECIDE,PETSC_NULL,&subx);
98: VecCreateMPIWithArray(subcomm,n,PETSC_DECIDE,PETSC_NULL,&subu);
100: VecGetArray(b,&barray);
101: VecGetArray(x,&xarray);
102: VecGetArray(u,&uarray);
103: VecPlaceArray(subb,barray);
104: VecPlaceArray(subx,xarray);
105: VecPlaceArray(subu,uarray);
107: /* Create linear solvers associated with subA */
108: KSPCreate(subcomm,&subksp);
109: KSPSetOperators(subksp,subA,subA,SAME_NONZERO_PATTERN);
110: KSPSetFromOptions(subksp);
112: /* Solve sub systems */
113: KSPSolve(subksp,subb,subx);
114: KSPGetIterationNumber(subksp,&its);
116: /* check residual */
117: MatMult(subA,subx,subu);
118: VecAXPY(subu,-1.0,subb);
119: VecNorm(u,NORM_2,&norm);
120: if (norm > 1.e-4 && !rank){
121: PetscPrintf(PETSC_COMM_WORLD,"[%D] Number of iterations = %3D\n",rank,its);
122: printf("Error: Residual norm of each block |subb - subA*subx |= %A\n",norm);
123: }
124: VecResetArray(subb);
125: VecResetArray(subx);
126: VecResetArray(subu);
127:
128: PetscOptionsGetInt(PETSC_NULL,"-subvec_view",&id,&flg);
129: if (flg && rank == id){
130: PetscPrintf(PETSC_COMM_SELF,"[%D] subb:\n", rank);
131: VecGetArray(subb,&array);
132: for (i=0; i<m; i++) printf("%G\n",PetscRealPart(array[i]));
133: VecRestoreArray(subb,&array);
134: PetscPrintf(PETSC_COMM_SELF,"[%D] subx:\n", rank);
135: VecGetArray(subx,&array);
136: for (i=0; i<m; i++) printf("%G\n",PetscRealPart(array[i]));
137: VecRestoreArray(subx,&array);
138: }
140: VecRestoreArray(x,&xarray);
141: VecRestoreArray(b,&barray);
142: VecRestoreArray(u,&uarray);
143: MatDestroy(&subA);
144: VecDestroy(&subb);
145: VecDestroy(&subx);
146: VecDestroy(&subu);
147: KSPDestroy(&subksp);
148: PetscSubcommDestroy(&psubcomm);
149: MatDestroy(&A); VecDestroy(&b);
150: VecDestroy(&u); VecDestroy(&x);
151:
152: PetscFinalize();
153: return 0;
154: }