Actual source code: ex214.c
petsc-3.9.4 2018-09-11
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;
13: Vec u,x,b;
14: PetscInt m,n,nrhs,M,N,i,Istart,Iend,Ii,j,J;
15: PetscScalar v;
16: PetscReal norm,tol=PETSC_SQRT_MACHINE_EPSILON;
17: PetscRandom rand;
18: PetscBool displ=PETSC_FALSE;
19: char solver[256];
20: #endif
22: PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
23: MPI_Comm_size(PETSC_COMM_WORLD,&size);
24: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
26: #if !defined(PETSC_HAVE_MUMPS)
27: if (!rank) {PetscPrintf(PETSC_COMM_SELF,"This example requires MUMPS, exit...\n");}
28: PetscFinalize();
29: return ierr;
30: #else
32: PetscOptionsGetBool(NULL,NULL,"-displ",&displ,NULL);
34: /* Create matrix A */
35: m = 4; n = 4;
36: PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL);
37: PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
39: MatCreate(PETSC_COMM_WORLD,&A);
40: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m*n,m*n);
41: MatSetFromOptions(A);
42: MatMPIAIJSetPreallocation(A,5,NULL,5,NULL);
43: MatSeqAIJSetPreallocation(A,5,NULL);
45: MatGetOwnershipRange(A,&Istart,&Iend);
46: for (Ii=Istart; Ii<Iend; Ii++) {
47: v = -1.0; i = Ii/n; j = Ii - i*n;
48: if (i>0) {J = Ii - n; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);}
49: if (i<m-1) {J = Ii + n; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);}
50: if (j>0) {J = Ii - 1; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);}
51: if (j<n-1) {J = Ii + 1; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);}
52: v = 4.0; MatSetValues(A,1,&Ii,1,&Ii,&v,ADD_VALUES);
53: }
54: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
55: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
57: MatGetLocalSize(A,&m,&n);
58: MatGetSize(A,&M,&N);
59: if (m != n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ, "This example is not intended for rectangular matrices (%d, %d)", m, n);
61: /* Create dense matrix C and X; C holds true solution with identical colums */
62: nrhs = N;
63: PetscOptionsGetInt(NULL,NULL,"-nrhs",&nrhs,NULL);
64: MatCreate(PETSC_COMM_WORLD,&C);
65: MatSetSizes(C,m,PETSC_DECIDE,PETSC_DECIDE,nrhs);
66: MatSetType(C,MATDENSE);
67: MatSetFromOptions(C);
68: MatSetUp(C);
70: PetscRandomCreate(PETSC_COMM_WORLD,&rand);
71: PetscRandomSetFromOptions(rand);
72: MatSetRandom(C,rand);
73: MatDuplicate(C,MAT_DO_NOT_COPY_VALUES,&X);
75: /* Create vectors */
76: VecCreate(PETSC_COMM_WORLD,&x);
77: VecSetSizes(x,n,PETSC_DECIDE);
78: VecSetFromOptions(x);
79: VecDuplicate(x,&b);
80: VecDuplicate(x,&u); /* save the true solution */
82: PetscStrcpy(solver,MATSOLVERMUMPS);
83: if (!rank && displ) {PetscPrintf(PETSC_COMM_SELF,"Solving with %s: nrhs %D, size mat %D x %D\n",solver,nrhs,M,N);}
85: /* Test LU Factorization */
86: MatGetFactor(A,solver,MAT_FACTOR_LU,&F);
87: MatLUFactorSymbolic(F,A,NULL,NULL,NULL);
89: VecSetRandom(x,rand);
90: MatLUFactorNumeric(F,A,NULL);
92: VecSetRandom(x,rand);
93: VecCopy(x,u);
95: /* (1) Test MatMatSolve(): dense RHS = A*C, C: true solutions */
96: MatMatMult(A,C,MAT_INITIAL_MATRIX,2.0,&RHS);
97: MatMatSolve(F,RHS,X);
99: /* Check the error */
100: MatAXPY(X,-1.0,C,SAME_NONZERO_PATTERN);
101: MatNorm(X,NORM_FROBENIUS,&norm);
102: if (norm > tol) {
103: PetscPrintf(PETSC_COMM_SELF,"(1) MatMatSolve: Norm of error %g\n",norm);
104: }
106: /* (2) Test MatMatSolve() for inv(A) with dense RHS:
107: RHS = [e[0],...,e[nrhs-1], dense X holds first nrhs columns of inv(A) */
108: MatZeroEntries(RHS);
109: for (i=0; i<nrhs; i++) {
110: v = 1.0;
111: MatSetValues(RHS,1,&i,1,&i,&v,INSERT_VALUES);
112: }
113: MatAssemblyBegin(RHS,MAT_FINAL_ASSEMBLY);
114: MatAssemblyEnd(RHS,MAT_FINAL_ASSEMBLY);
116: MatMatSolve(F,RHS,X);
117: if (displ) {
118: if (!rank) {PetscPrintf(PETSC_COMM_SELF," (2) first %D columns of inv(A) with dense RHS:\n",nrhs);}
119: MatView(X,PETSC_VIEWER_STDOUT_WORLD);
120: }
122: /* Check the residual */
123: MatMatMult(A,X,MAT_INITIAL_MATRIX,2.0,&AX);
124: MatAXPY(AX,-1.0,RHS,SAME_NONZERO_PATTERN);
125: MatNorm(AX,NORM_INFINITY,&norm);
126: if (norm > tol) {
127: PetscPrintf(PETSC_COMM_SELF,"(2) MatMatSolve: Norm of residual %g\n",norm);
128: }
130: if (size == 1) {
131: /* (3) Test MatMatSolve() for inv(A) with sparse RHS:
132: spRHS = [e[0],...,e[nrhs-1], dense X holds first nrhs columns of inv(A) */
133: Mat RHST,spRHST,spRHS;
135: MatTranspose(RHS,MAT_INITIAL_MATRIX,&RHST);
136: MatConvert(RHST,MATAIJ,MAT_INITIAL_MATRIX,&spRHST);
138: /* MUMPS requres spRHS in compressed column format, which PETSc does not support.
139: PETSc can create a new matrix object that shares same data structure with A and behaves like A^T */
140: MatCreateTranspose(spRHST,&spRHS);
141: MatMatSolve(F,spRHS,X);
142: if (displ) {
143: if (!rank) {PetscPrintf(PETSC_COMM_SELF," (3) first %D columns of inv(A) with sparse RHS:\n",nrhs);}
144: MatView(X,PETSC_VIEWER_STDOUT_SELF);
145: }
147: /* Check the residual */
148: MatMatMult(A,X,MAT_REUSE_MATRIX,2.0,&AX);
149: MatAXPY(AX,-1.0,RHS,SAME_NONZERO_PATTERN);
150: MatNorm(AX,NORM_INFINITY,&norm);
151: if (norm > tol) {
152: PetscPrintf(PETSC_COMM_SELF,"(3) MatMatSolve: Norm of residual %g\n",norm);
153: }
155: MatDestroy(&spRHST);
156: MatDestroy(&spRHS);
157: MatDestroy(&RHST);
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: if (nrhs == N) { /* mumps requires nrhs = n */
163: /* Create spRHS on proc[0] */
164: Mat spRHS = NULL,spRHST;
165: if (!rank) {
166: /* Create spRHST = spRHS^T in compressed row format (aij) and set its nonzero structure */
167: MatCreateSeqAIJ(PETSC_COMM_SELF,N,N,2,NULL,&spRHST);
168: v = 0.0;
169: for (i=0; i<N; i++) {
170: MatSetValues(spRHST,1,&i,1,&i,&v,INSERT_VALUES);
171: }
172: for (i=1; i<N; i++) {
173: j = i - 1;
174: MatSetValues(spRHST,1,&i,1,&j,&v,INSERT_VALUES);
175: }
176: MatAssemblyBegin(spRHST,MAT_FINAL_ASSEMBLY);
177: MatAssemblyEnd(spRHST,MAT_FINAL_ASSEMBLY);
179: /* Create spRHS = spRHST^T. Two matrices share internal matrix data structure */
180: MatCreateTranspose(spRHST,&spRHS);
181: }
183: MatMumpsGetInverse(F,spRHS);
185: if (!rank) {
186: if (displ) {
187: /* spRHST and spRHS share same internal data structure */
188: Mat spRHSTT;
189: MatTranspose(spRHST,MAT_INITIAL_MATRIX,&spRHSTT);
190: PetscPrintf(PETSC_COMM_SELF,"\nSelected entries of inv(A):\n");
191: MatView(spRHSTT,0);
192: MatDestroy(&spRHSTT);
193: }
195: MatDestroy(&spRHST);
196: MatDestroy(&spRHS);
197: }
198: }
200: /* Free data structures */
201: MatDestroy(&AX);
202: MatDestroy(&A);
203: MatDestroy(&C);
204: MatDestroy(&F);
205: MatDestroy(&X);
206: MatDestroy(&RHS);
207: PetscRandomDestroy(&rand);
208: VecDestroy(&x);
209: VecDestroy(&b);
210: VecDestroy(&u);
211: PetscFinalize();
212: return ierr;
213: #endif
214: }
216: /*TEST
219: test:
220: requires: mumps double !complex !define(PETSC_USE_64BIT_INDICES)
222: test:
223: suffix: 2
224: requires: mumps double !complex !define(PETSC_USE_64BIT_INDICES)
225: nsize: 2
227: TEST*/