Actual source code: ex214.c
petsc-3.10.5 2019-03-28
2: static char help[] = "Tests MatMatSolve() for computing inv(A) with MUMPS.\n\
3: Example: mpiexec -n <np> ./ex214 -displ \n\n";
5: #include <petscmat.h>
7: int main(int argc,char **args)
8: {
10: PetscMPIInt size,rank;
11: #if defined(PETSC_HAVE_MUMPS)
12: Mat A,RHS,C,F,X,AX,spRHST;
13: PetscInt m,n,nrhs,M,N,i,Istart,Iend,Ii,j,J;
14: PetscScalar v;
15: PetscReal norm,tol=PETSC_SQRT_MACHINE_EPSILON;
16: PetscRandom rand;
17: PetscBool displ=PETSC_FALSE;
18: char solver[256];
19: #endif
21: PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
22: MPI_Comm_size(PETSC_COMM_WORLD,&size);
23: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
25: #if !defined(PETSC_HAVE_MUMPS)
26: if (!rank) {PetscPrintf(PETSC_COMM_SELF,"This example requires MUMPS, exit...\n");}
27: PetscFinalize();
28: return ierr;
29: #else
31: PetscOptionsGetBool(NULL,NULL,"-displ",&displ,NULL);
33: /* Create matrix A */
34: m = 4; n = 4;
35: PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL);
36: PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
38: MatCreate(PETSC_COMM_WORLD,&A);
39: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m*n,m*n);
40: MatSetFromOptions(A);
41: MatMPIAIJSetPreallocation(A,5,NULL,5,NULL);
42: MatSeqAIJSetPreallocation(A,5,NULL);
44: MatGetOwnershipRange(A,&Istart,&Iend);
45: for (Ii=Istart; Ii<Iend; Ii++) {
46: v = -1.0; i = Ii/n; j = Ii - i*n;
47: if (i>0) {J = Ii - n; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);}
48: if (i<m-1) {J = Ii + n; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);}
49: if (j>0) {J = Ii - 1; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);}
50: if (j<n-1) {J = Ii + 1; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);}
51: v = 4.0; MatSetValues(A,1,&Ii,1,&Ii,&v,ADD_VALUES);
52: }
53: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
54: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
56: MatGetLocalSize(A,&m,&n);
57: MatGetSize(A,&M,&N);
58: if (m != n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ, "This example is not intended for rectangular matrices (%d, %d)", m, n);
60: /* Create dense matrix C and X; C holds true solution with identical colums */
61: nrhs = N;
62: PetscOptionsGetInt(NULL,NULL,"-nrhs",&nrhs,NULL);
63: MatCreate(PETSC_COMM_WORLD,&C);
64: MatSetSizes(C,m,PETSC_DECIDE,PETSC_DECIDE,nrhs);
65: MatSetType(C,MATDENSE);
66: MatSetFromOptions(C);
67: MatSetUp(C);
69: PetscRandomCreate(PETSC_COMM_WORLD,&rand);
70: PetscRandomSetFromOptions(rand);
71: MatSetRandom(C,rand);
72: MatDuplicate(C,MAT_DO_NOT_COPY_VALUES,&X);
74: PetscStrcpy(solver,MATSOLVERMUMPS);
75: if (!rank && displ) {PetscPrintf(PETSC_COMM_SELF,"Solving with %s: nrhs %D, size mat %D x %D\n",solver,nrhs,M,N);}
77: /* Test LU Factorization */
78: MatGetFactor(A,solver,MAT_FACTOR_LU,&F);
79: MatLUFactorSymbolic(F,A,NULL,NULL,NULL);
80: MatLUFactorNumeric(F,A,NULL);
82: /* (1) Test MatMatSolve(): dense RHS = A*C, C: true solutions */
83: /* ---------------------------------------------------------- */
84: MatMatMult(A,C,MAT_INITIAL_MATRIX,2.0,&RHS);
85: MatMatSolve(F,RHS,X);
87: /* Check the error */
88: MatAXPY(X,-1.0,C,SAME_NONZERO_PATTERN);
89: MatNorm(X,NORM_FROBENIUS,&norm);
90: if (norm > tol) {
91: PetscPrintf(PETSC_COMM_SELF,"(1) MatMatSolve: Norm of error %g\n",norm);
92: }
94: /* (2) Test MatMatSolve() for inv(A) with dense RHS:
95: RHS = [e[0],...,e[nrhs-1]], dense X holds first nrhs columns of inv(A) */
96: /* -------------------------------------------------------------------- */
97: MatZeroEntries(RHS);
98: for (i=0; i<nrhs; i++) {
99: v = 1.0;
100: MatSetValues(RHS,1,&i,1,&i,&v,INSERT_VALUES);
101: }
102: MatAssemblyBegin(RHS,MAT_FINAL_ASSEMBLY);
103: MatAssemblyEnd(RHS,MAT_FINAL_ASSEMBLY);
105: MatMatSolve(F,RHS,X);
106: if (displ) {
107: if (!rank) {PetscPrintf(PETSC_COMM_SELF," \n(2) first %D columns of inv(A) with dense RHS:\n",nrhs);}
108: MatView(X,PETSC_VIEWER_STDOUT_WORLD);
109: }
111: /* Check the residual */
112: MatMatMult(A,X,MAT_INITIAL_MATRIX,2.0,&AX);
113: MatAXPY(AX,-1.0,RHS,SAME_NONZERO_PATTERN);
114: MatNorm(AX,NORM_INFINITY,&norm);
115: if (norm > tol) {
116: PetscPrintf(PETSC_COMM_SELF,"(2) MatMatSolve: Norm of residual %g\n",norm);
117: }
118: MatZeroEntries(X);
120: /* (3) Test MatMatTransposeSolve() for inv(A) with sparse RHS stored in the host:
121: spRHST = [e[0],...,e[nrhs-1]]^T, dense X holds first nrhs columns of inv(A) */
122: /* --------------------------------------------------------------------------- */
123: /* Create spRHST: PETSc does not support compressed column format which is required by MUMPS for sparse RHS matrix,
124: thus user must create spRHST=spRHS^T and call MatMatTransposeSolve() */
125: MatCreate(PETSC_COMM_WORLD,&spRHST);
126: if (!rank) {
127: /* MUMPS requirs RHS be centralized on the host! */
128: MatSetSizes(spRHST,nrhs,M,PETSC_DECIDE,PETSC_DECIDE);
129: } else {
130: MatSetSizes(spRHST,0,0,PETSC_DECIDE,PETSC_DECIDE);
131: }
132: MatSetType(spRHST,MATAIJ);
133: MatSetFromOptions(spRHST);
134: MatSetUp(spRHST);
135: if (!rank) {
136: v = 1.0;
137: for (i=0; i<nrhs; i++) {
138: MatSetValues(spRHST,1,&i,1,&i,&v,INSERT_VALUES);
139: }
140: }
141: MatAssemblyBegin(spRHST,MAT_FINAL_ASSEMBLY);
142: MatAssemblyEnd(spRHST,MAT_FINAL_ASSEMBLY);
144: MatMatTransposeSolve(F,spRHST,X);
146: if (displ) {
147: if (!rank) {PetscPrintf(PETSC_COMM_SELF," \n(3) first %D columns of inv(A) with sparse RHS:\n",nrhs);}
148: MatView(X,PETSC_VIEWER_STDOUT_WORLD);
149: }
151: /* Check the residual */
152: MatMatMult(A,X,MAT_REUSE_MATRIX,2.0,&AX);
154: MatAXPY(AX,-1.0,RHS,SAME_NONZERO_PATTERN);
155: MatNorm(AX,NORM_INFINITY,&norm);
156: if (norm > tol) {
157: PetscPrintf(PETSC_COMM_SELF,"(3) MatMatSolve: Norm of residual %g\n",norm);
158: }
160: /* (4) Test MatMatSolve() for inv(A) with selected entries:
161: input: spRHS gives selected indices; output: spRHS holds selected entries of inv(A) */
162: /* --------------------------------------------------------------------------------- */
163: if (nrhs == N) { /* mumps requires nrhs = n */
164: /* Create spRHS on proc[0] */
165: Mat spRHS = NULL;
167: /* Create spRHS = spRHST^T. Two matrices share internal matrix data structure */
168: MatCreateTranspose(spRHST,&spRHS);
169: MatMumpsGetInverse(F,spRHS);
170: MatDestroy(&spRHS);
172: MatMumpsGetInverseTranspose(F,spRHST);
173: if (displ) {
174: PetscPrintf(PETSC_COMM_SELF,"\nSelected entries of inv(A^T):\n");
175: MatView(spRHST,PETSC_VIEWER_STDOUT_WORLD);
176: }
177: MatDestroy(&spRHS);
178: }
180: /* Free data structures */
181: MatDestroy(&AX);
182: MatDestroy(&A);
183: MatDestroy(&C);
184: MatDestroy(&F);
185: MatDestroy(&X);
186: MatDestroy(&RHS);
187: MatDestroy(&spRHST);
188: PetscRandomDestroy(&rand);
189: PetscFinalize();
190: return ierr;
191: #endif
192: }
194: /*TEST
197: test:
198: requires: mumps double !complex !define(PETSC_USE_64BIT_INDICES)
200: test:
201: suffix: 2
202: requires: mumps double !complex !define(PETSC_USE_64BIT_INDICES)
203: nsize: 2
205: TEST*/