Actual source code: ex214.c
petsc-3.12.5 2020-03-29
2: static char help[] = "Tests MatMatSolve() and MatMatTransposeSolve() 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,test;
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: for (test=0; test<2; test++) {
78: if (test == 0) {
79: /* Test LU Factorization */
80: PetscPrintf(PETSC_COMM_WORLD,"test LU factorization\n");
81: MatGetFactor(A,solver,MAT_FACTOR_LU,&F);
82: MatLUFactorSymbolic(F,A,NULL,NULL,NULL);
83: MatLUFactorNumeric(F,A,NULL);
84: } else {
85: /* Test Cholesky Factorization */
86: PetscBool flg;
87: MatIsSymmetric(A,0.0,&flg);
88: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"A must be symmetric for Cholesky factorization");
90: PetscPrintf(PETSC_COMM_WORLD,"test Cholesky factorization\n");
91: MatGetFactor(A,solver,MAT_FACTOR_CHOLESKY,&F);
92: MatCholeskyFactorSymbolic(F,A,NULL,NULL);
93: MatCholeskyFactorNumeric(F,A,NULL);
94: }
96: /* (1) Test MatMatSolve(): dense RHS = A*C, C: true solutions */
97: /* ---------------------------------------------------------- */
98: MatMatMult(A,C,MAT_INITIAL_MATRIX,2.0,&RHS);
99: MatMatSolve(F,RHS,X);
101: /* Check the error */
102: MatAXPY(X,-1.0,C,SAME_NONZERO_PATTERN);
103: MatNorm(X,NORM_FROBENIUS,&norm);
104: if (norm > tol) {
105: PetscPrintf(PETSC_COMM_SELF,"(1) MatMatSolve: Norm of error %g\n",norm);
106: }
108: /* (2) Test MatMatSolve() for inv(A) with dense RHS:
109: RHS = [e[0],...,e[nrhs-1]], dense X holds first nrhs columns of inv(A) */
110: /* -------------------------------------------------------------------- */
111: MatZeroEntries(RHS);
112: for (i=0; i<nrhs; i++) {
113: v = 1.0;
114: MatSetValues(RHS,1,&i,1,&i,&v,INSERT_VALUES);
115: }
116: MatAssemblyBegin(RHS,MAT_FINAL_ASSEMBLY);
117: MatAssemblyEnd(RHS,MAT_FINAL_ASSEMBLY);
119: MatMatSolve(F,RHS,X);
120: if (displ) {
121: if (!rank) {PetscPrintf(PETSC_COMM_SELF," \n(2) first %D columns of inv(A) with dense RHS:\n",nrhs);}
122: MatView(X,PETSC_VIEWER_STDOUT_WORLD);
123: }
125: /* Check the residual */
126: MatMatMult(A,X,MAT_INITIAL_MATRIX,2.0,&AX);
127: MatAXPY(AX,-1.0,RHS,SAME_NONZERO_PATTERN);
128: MatNorm(AX,NORM_INFINITY,&norm);
129: if (norm > tol) {
130: PetscPrintf(PETSC_COMM_SELF,"(2) MatMatSolve: Norm of residual %g\n",norm);
131: }
132: MatZeroEntries(X);
134: /* (3) Test MatMatTransposeSolve() for inv(A) with sparse RHS stored in the host:
135: spRHST = [e[0],...,e[nrhs-1]]^T, dense X holds first nrhs columns of inv(A) */
136: /* --------------------------------------------------------------------------- */
137: /* Create spRHST: PETSc does not support compressed column format which is required by MUMPS for sparse RHS matrix,
138: thus user must create spRHST=spRHS^T and call MatMatTransposeSolve() */
139: MatCreate(PETSC_COMM_WORLD,&spRHST);
140: if (!rank) {
141: /* MUMPS requirs RHS be centralized on the host! */
142: MatSetSizes(spRHST,nrhs,M,PETSC_DECIDE,PETSC_DECIDE);
143: } else {
144: MatSetSizes(spRHST,0,0,PETSC_DECIDE,PETSC_DECIDE);
145: }
146: MatSetType(spRHST,MATAIJ);
147: MatSetFromOptions(spRHST);
148: MatSetUp(spRHST);
149: if (!rank) {
150: v = 1.0;
151: for (i=0; i<nrhs; i++) {
152: MatSetValues(spRHST,1,&i,1,&i,&v,INSERT_VALUES);
153: }
154: }
155: MatAssemblyBegin(spRHST,MAT_FINAL_ASSEMBLY);
156: MatAssemblyEnd(spRHST,MAT_FINAL_ASSEMBLY);
158: MatMatTransposeSolve(F,spRHST,X);
160: if (displ) {
161: if (!rank) {PetscPrintf(PETSC_COMM_SELF," \n(3) first %D columns of inv(A) with sparse RHS:\n",nrhs);}
162: MatView(X,PETSC_VIEWER_STDOUT_WORLD);
163: }
165: /* Check the residual: R = A*X - RHS */
166: MatMatMult(A,X,MAT_REUSE_MATRIX,2.0,&AX);
168: MatAXPY(AX,-1.0,RHS,SAME_NONZERO_PATTERN);
169: MatNorm(AX,NORM_INFINITY,&norm);
170: if (norm > tol) {
171: PetscPrintf(PETSC_COMM_SELF,"(3) MatMatSolve: Norm of residual %g\n",norm);
172: }
174: /* (4) Test MatMatSolve() for inv(A) with selected entries:
175: input: spRHS gives selected indices; output: spRHS holds selected entries of inv(A) */
176: /* --------------------------------------------------------------------------------- */
177: if (nrhs == N) { /* mumps requires nrhs = n */
178: /* Create spRHS on proc[0] */
179: Mat spRHS = NULL;
181: /* Create spRHS = spRHST^T. Two matrices share internal matrix data structure */
182: MatCreateTranspose(spRHST,&spRHS);
183: MatMumpsGetInverse(F,spRHS);
184: MatDestroy(&spRHS);
186: MatMumpsGetInverseTranspose(F,spRHST);
187: if (displ) {
188: PetscPrintf(PETSC_COMM_WORLD,"\nSelected entries of inv(A^T):\n");
189: MatView(spRHST,PETSC_VIEWER_STDOUT_WORLD);
190: }
191: MatDestroy(&spRHS);
192: }
193: MatDestroy(&AX);
194: MatDestroy(&F);
195: MatDestroy(&RHS);
196: MatDestroy(&spRHST);
197: }
199: /* Free data structures */
200: MatDestroy(&A);
201: MatDestroy(&C);
202: MatDestroy(&X);
203: PetscRandomDestroy(&rand);
204: PetscFinalize();
205: return ierr;
206: #endif
207: }
209: /*TEST
212: test:
213: requires: mumps double !complex !define(PETSC_USE_64BIT_INDICES)
215: test:
216: suffix: 2
217: requires: mumps double !complex !define(PETSC_USE_64BIT_INDICES)
218: nsize: 2
220: TEST*/