Actual source code: ex214.c
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: {
9: PetscMPIInt size,rank;
10: #if defined(PETSC_HAVE_MUMPS)
11: Mat A,RHS,C,F,X,AX,spRHST;
12: PetscInt m,n,nrhs,M,N,i,Istart,Iend,Ii,j,J,test;
13: PetscScalar v;
14: PetscReal norm,tol=PETSC_SQRT_MACHINE_EPSILON;
15: PetscRandom rand;
16: PetscBool displ=PETSC_FALSE;
17: char solver[256];
18: #endif
20: PetscInitialize(&argc,&args,(char*)0,help);
21: MPI_Comm_size(PETSC_COMM_WORLD,&size);
22: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
24: #if !defined(PETSC_HAVE_MUMPS)
25: if (rank == 0) PetscPrintf(PETSC_COMM_SELF,"This example requires MUMPS, exit...\n");
26: PetscFinalize();
27: return 0;
28: #else
30: PetscOptionsGetBool(NULL,NULL,"-displ",&displ,NULL);
32: /* Create matrix A */
33: m = 4; n = 4;
34: PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL);
35: PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
37: MatCreate(PETSC_COMM_WORLD,&A);
38: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m*n,m*n);
39: MatSetFromOptions(A);
40: MatMPIAIJSetPreallocation(A,5,NULL,5,NULL);
41: MatSeqAIJSetPreallocation(A,5,NULL);
43: MatGetOwnershipRange(A,&Istart,&Iend);
44: for (Ii=Istart; Ii<Iend; Ii++) {
45: v = -1.0; i = Ii/n; j = Ii - i*n;
46: if (i>0) {J = Ii - n; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);}
47: if (i<m-1) {J = Ii + n; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);}
48: if (j>0) {J = Ii - 1; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);}
49: if (j<n-1) {J = Ii + 1; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);}
50: v = 4.0; MatSetValues(A,1,&Ii,1,&Ii,&v,ADD_VALUES);
51: }
52: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
53: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
55: MatGetLocalSize(A,&m,&n);
56: MatGetSize(A,&M,&N);
59: /* Create dense matrix C and X; C holds true solution with identical columns */
60: nrhs = N;
61: PetscOptionsGetInt(NULL,NULL,"-nrhs",&nrhs,NULL);
62: MatCreate(PETSC_COMM_WORLD,&C);
63: MatSetSizes(C,m,PETSC_DECIDE,PETSC_DECIDE,nrhs);
64: MatSetType(C,MATDENSE);
65: MatSetFromOptions(C);
66: MatSetUp(C);
68: PetscRandomCreate(PETSC_COMM_WORLD,&rand);
69: PetscRandomSetFromOptions(rand);
70: MatSetRandom(C,rand);
71: MatDuplicate(C,MAT_DO_NOT_COPY_VALUES,&X);
73: PetscStrcpy(solver,MATSOLVERMUMPS);
74: if (rank == 0 && displ) PetscPrintf(PETSC_COMM_SELF,"Solving with %s: nrhs %" PetscInt_FMT ", size mat %" PetscInt_FMT " x %" PetscInt_FMT "\n",solver,nrhs,M,N);
76: for (test=0; test<2; test++) {
77: if (test == 0) {
78: /* Test LU Factorization */
79: PetscPrintf(PETSC_COMM_WORLD,"test LU factorization\n");
80: MatGetFactor(A,solver,MAT_FACTOR_LU,&F);
81: MatLUFactorSymbolic(F,A,NULL,NULL,NULL);
82: MatLUFactorNumeric(F,A,NULL);
83: } else {
84: /* Test Cholesky Factorization */
85: PetscBool flg;
86: MatIsSymmetric(A,0.0,&flg);
89: PetscPrintf(PETSC_COMM_WORLD,"test Cholesky factorization\n");
90: MatGetFactor(A,solver,MAT_FACTOR_CHOLESKY,&F);
91: MatCholeskyFactorSymbolic(F,A,NULL,NULL);
92: MatCholeskyFactorNumeric(F,A,NULL);
93: }
95: /* (1) Test MatMatSolve(): dense RHS = A*C, C: true solutions */
96: /* ---------------------------------------------------------- */
97: MatMatMult(A,C,MAT_INITIAL_MATRIX,2.0,&RHS);
98: 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: /* Test X=RHS */
107: MatMatSolve(F,RHS,RHS);
108: /* Check the error */
109: MatAXPY(RHS,-1.0,C,SAME_NONZERO_PATTERN);
110: MatNorm(RHS,NORM_FROBENIUS,&norm);
111: if (norm > tol) {
112: PetscPrintf(PETSC_COMM_SELF,"(1.1) MatMatSolve: Norm of error %g\n",norm);
113: }
115: /* (2) Test MatMatSolve() for inv(A) with dense RHS:
116: RHS = [e[0],...,e[nrhs-1]], dense X holds first nrhs columns of inv(A) */
117: /* -------------------------------------------------------------------- */
118: MatZeroEntries(RHS);
119: for (i=0; i<nrhs; i++) {
120: v = 1.0;
121: MatSetValues(RHS,1,&i,1,&i,&v,INSERT_VALUES);
122: }
123: MatAssemblyBegin(RHS,MAT_FINAL_ASSEMBLY);
124: MatAssemblyEnd(RHS,MAT_FINAL_ASSEMBLY);
126: MatMatSolve(F,RHS,X);
127: if (displ) {
128: if (rank == 0) PetscPrintf(PETSC_COMM_SELF," \n(2) first %" PetscInt_FMT " columns of inv(A) with dense RHS:\n",nrhs);
129: MatView(X,PETSC_VIEWER_STDOUT_WORLD);
130: }
132: /* Check the residual */
133: MatMatMult(A,X,MAT_INITIAL_MATRIX,2.0,&AX);
134: MatAXPY(AX,-1.0,RHS,SAME_NONZERO_PATTERN);
135: MatNorm(AX,NORM_INFINITY,&norm);
136: if (norm > tol) {
137: PetscPrintf(PETSC_COMM_SELF,"(2) MatMatSolve: Norm of residual %g\n",norm);
138: }
139: MatZeroEntries(X);
141: /* (3) Test MatMatTransposeSolve() for inv(A) with sparse RHS stored in the host:
142: spRHST = [e[0],...,e[nrhs-1]]^T, dense X holds first nrhs columns of inv(A) */
143: /* --------------------------------------------------------------------------- */
144: /* Create spRHST: PETSc does not support compressed column format which is required by MUMPS for sparse RHS matrix,
145: thus user must create spRHST=spRHS^T and call MatMatTransposeSolve() */
146: MatCreate(PETSC_COMM_WORLD,&spRHST);
147: if (rank == 0) {
148: /* MUMPS requirs RHS be centralized on the host! */
149: MatSetSizes(spRHST,nrhs,M,PETSC_DECIDE,PETSC_DECIDE);
150: } else {
151: MatSetSizes(spRHST,0,0,PETSC_DECIDE,PETSC_DECIDE);
152: }
153: MatSetType(spRHST,MATAIJ);
154: MatSetFromOptions(spRHST);
155: MatSetUp(spRHST);
156: if (rank == 0) {
157: v = 1.0;
158: for (i=0; i<nrhs; i++) {
159: MatSetValues(spRHST,1,&i,1,&i,&v,INSERT_VALUES);
160: }
161: }
162: MatAssemblyBegin(spRHST,MAT_FINAL_ASSEMBLY);
163: MatAssemblyEnd(spRHST,MAT_FINAL_ASSEMBLY);
165: MatMatTransposeSolve(F,spRHST,X);
167: if (displ) {
168: if (rank == 0) PetscPrintf(PETSC_COMM_SELF," \n(3) first %" PetscInt_FMT " columns of inv(A) with sparse RHS:\n",nrhs);
169: MatView(X,PETSC_VIEWER_STDOUT_WORLD);
170: }
172: /* Check the residual: R = A*X - RHS */
173: MatMatMult(A,X,MAT_REUSE_MATRIX,2.0,&AX);
175: MatAXPY(AX,-1.0,RHS,SAME_NONZERO_PATTERN);
176: MatNorm(AX,NORM_INFINITY,&norm);
177: if (norm > tol) {
178: PetscPrintf(PETSC_COMM_SELF,"(3) MatMatSolve: Norm of residual %g\n",norm);
179: }
181: /* (4) Test MatMatSolve() for inv(A) with selected entries:
182: input: spRHS gives selected indices; output: spRHS holds selected entries of inv(A) */
183: /* --------------------------------------------------------------------------------- */
184: if (nrhs == N) { /* mumps requires nrhs = n */
185: /* Create spRHS on proc[0] */
186: Mat spRHS = NULL;
188: /* Create spRHS = spRHST^T. Two matrices share internal matrix data structure */
189: MatCreateTranspose(spRHST,&spRHS);
190: MatMumpsGetInverse(F,spRHS);
191: MatDestroy(&spRHS);
193: MatMumpsGetInverseTranspose(F,spRHST);
194: if (displ) {
195: PetscPrintf(PETSC_COMM_WORLD,"\nSelected entries of inv(A^T):\n");
196: MatView(spRHST,PETSC_VIEWER_STDOUT_WORLD);
197: }
198: MatDestroy(&spRHS);
199: }
200: MatDestroy(&AX);
201: MatDestroy(&F);
202: MatDestroy(&RHS);
203: MatDestroy(&spRHST);
204: }
206: /* Free data structures */
207: MatDestroy(&A);
208: MatDestroy(&C);
209: MatDestroy(&X);
210: PetscRandomDestroy(&rand);
211: PetscFinalize();
212: return 0;
213: #endif
214: }
216: /*TEST
218: test:
219: requires: mumps double !complex
221: test:
222: suffix: 2
223: requires: mumps double !complex
224: nsize: 2
226: TEST*/