Actual source code: ex214.c
petsc-3.13.6 2020-09-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);
100: /* Check the error */
101: MatAXPY(X,-1.0,C,SAME_NONZERO_PATTERN);
102: MatNorm(X,NORM_FROBENIUS,&norm);
103: if (norm > tol) {
104: PetscPrintf(PETSC_COMM_SELF,"(1) MatMatSolve: Norm of error %g\n",norm);
105: }
107: /* Test X=RHS */
108: MatMatSolve(F,RHS,RHS);
109: /* Check the error */
110: MatAXPY(RHS,-1.0,C,SAME_NONZERO_PATTERN);
111: MatNorm(RHS,NORM_FROBENIUS,&norm);
112: if (norm > tol) {
113: PetscPrintf(PETSC_COMM_SELF,"(1.1) MatMatSolve: Norm of error %g\n",norm);
114: }
116: /* (2) Test MatMatSolve() for inv(A) with dense RHS:
117: RHS = [e[0],...,e[nrhs-1]], dense X holds first nrhs columns of inv(A) */
118: /* -------------------------------------------------------------------- */
119: MatZeroEntries(RHS);
120: for (i=0; i<nrhs; i++) {
121: v = 1.0;
122: MatSetValues(RHS,1,&i,1,&i,&v,INSERT_VALUES);
123: }
124: MatAssemblyBegin(RHS,MAT_FINAL_ASSEMBLY);
125: MatAssemblyEnd(RHS,MAT_FINAL_ASSEMBLY);
127: MatMatSolve(F,RHS,X);
128: if (displ) {
129: if (!rank) {PetscPrintf(PETSC_COMM_SELF," \n(2) first %D columns of inv(A) with dense RHS:\n",nrhs);}
130: MatView(X,PETSC_VIEWER_STDOUT_WORLD);
131: }
133: /* Check the residual */
134: MatMatMult(A,X,MAT_INITIAL_MATRIX,2.0,&AX);
135: MatAXPY(AX,-1.0,RHS,SAME_NONZERO_PATTERN);
136: MatNorm(AX,NORM_INFINITY,&norm);
137: if (norm > tol) {
138: PetscPrintf(PETSC_COMM_SELF,"(2) MatMatSolve: Norm of residual %g\n",norm);
139: }
140: MatZeroEntries(X);
142: /* (3) Test MatMatTransposeSolve() for inv(A) with sparse RHS stored in the host:
143: spRHST = [e[0],...,e[nrhs-1]]^T, dense X holds first nrhs columns of inv(A) */
144: /* --------------------------------------------------------------------------- */
145: /* Create spRHST: PETSc does not support compressed column format which is required by MUMPS for sparse RHS matrix,
146: thus user must create spRHST=spRHS^T and call MatMatTransposeSolve() */
147: MatCreate(PETSC_COMM_WORLD,&spRHST);
148: if (!rank) {
149: /* MUMPS requirs RHS be centralized on the host! */
150: MatSetSizes(spRHST,nrhs,M,PETSC_DECIDE,PETSC_DECIDE);
151: } else {
152: MatSetSizes(spRHST,0,0,PETSC_DECIDE,PETSC_DECIDE);
153: }
154: MatSetType(spRHST,MATAIJ);
155: MatSetFromOptions(spRHST);
156: MatSetUp(spRHST);
157: if (!rank) {
158: v = 1.0;
159: for (i=0; i<nrhs; i++) {
160: MatSetValues(spRHST,1,&i,1,&i,&v,INSERT_VALUES);
161: }
162: }
163: MatAssemblyBegin(spRHST,MAT_FINAL_ASSEMBLY);
164: MatAssemblyEnd(spRHST,MAT_FINAL_ASSEMBLY);
166: MatMatTransposeSolve(F,spRHST,X);
168: if (displ) {
169: if (!rank) {PetscPrintf(PETSC_COMM_SELF," \n(3) first %D columns of inv(A) with sparse RHS:\n",nrhs);}
170: MatView(X,PETSC_VIEWER_STDOUT_WORLD);
171: }
173: /* Check the residual: R = A*X - RHS */
174: MatMatMult(A,X,MAT_REUSE_MATRIX,2.0,&AX);
176: MatAXPY(AX,-1.0,RHS,SAME_NONZERO_PATTERN);
177: MatNorm(AX,NORM_INFINITY,&norm);
178: if (norm > tol) {
179: PetscPrintf(PETSC_COMM_SELF,"(3) MatMatSolve: Norm of residual %g\n",norm);
180: }
182: /* (4) Test MatMatSolve() for inv(A) with selected entries:
183: input: spRHS gives selected indices; output: spRHS holds selected entries of inv(A) */
184: /* --------------------------------------------------------------------------------- */
185: if (nrhs == N) { /* mumps requires nrhs = n */
186: /* Create spRHS on proc[0] */
187: Mat spRHS = NULL;
189: /* Create spRHS = spRHST^T. Two matrices share internal matrix data structure */
190: MatCreateTranspose(spRHST,&spRHS);
191: MatMumpsGetInverse(F,spRHS);
192: MatDestroy(&spRHS);
194: MatMumpsGetInverseTranspose(F,spRHST);
195: if (displ) {
196: PetscPrintf(PETSC_COMM_WORLD,"\nSelected entries of inv(A^T):\n");
197: MatView(spRHST,PETSC_VIEWER_STDOUT_WORLD);
198: }
199: MatDestroy(&spRHS);
200: }
201: MatDestroy(&AX);
202: MatDestroy(&F);
203: MatDestroy(&RHS);
204: MatDestroy(&spRHST);
205: }
207: /* Free data structures */
208: MatDestroy(&A);
209: MatDestroy(&C);
210: MatDestroy(&X);
211: PetscRandomDestroy(&rand);
212: PetscFinalize();
213: return ierr;
214: #endif
215: }
217: /*TEST
220: test:
221: requires: mumps double !complex
223: test:
224: suffix: 2
225: requires: mumps double !complex
226: nsize: 2
228: TEST*/