Actual source code: ex125.c
petsc-3.13.6 2020-09-29
1: static char help[] = "Tests MatSolve() and MatMatSolve() (interface to superlu_dist, mumps and mkl_pardiso).\n\
2: Example: mpiexec -n <np> ./ex125 -f <matrix binary file> -nrhs 4 \n\n";
4: #include <petscmat.h>
6: int main(int argc,char **args)
7: {
8: Mat A,RHS,C,F,X;
9: Vec u,x,b;
11: PetscMPIInt size;
12: PetscInt m,n,nfact,nsolve,nrhs,ipack=0;
13: PetscReal norm,tol=1.e-10;
14: IS perm,iperm;
15: MatFactorInfo info;
16: PetscRandom rand;
17: PetscBool flg,testMatSolve=PETSC_TRUE,testMatMatSolve=PETSC_TRUE;
18: PetscBool chol=PETSC_FALSE,view=PETSC_FALSE,matsolvexx = PETSC_FALSE;
19: #if defined(PETSC_HAVE_MUMPS)
20: PetscBool test_mumps_opts=PETSC_FALSE;
21: #endif
22: PetscViewer fd; /* viewer */
23: char file[PETSC_MAX_PATH_LEN]; /* input file name */
25: PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
26: MPI_Comm_size(PETSC_COMM_WORLD, &size);
28: /* Determine file from which we read the matrix A */
29: PetscOptionsGetString(NULL,NULL,"-f",file,PETSC_MAX_PATH_LEN,&flg);
30: if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"Must indicate binary file with the -f option");
32: /* Load matrix A */
33: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);
34: MatCreate(PETSC_COMM_WORLD,&A);
35: MatSetFromOptions(A);
36: MatLoad(A,fd);
37: PetscViewerDestroy(&fd);
38: MatGetLocalSize(A,&m,&n);
39: if (m != n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ, "This example is not intended for rectangular matrices (%D, %D)", m, n);
41: /* if A is symmetric, set its flag -- required by MatGetInertia() */
42: MatIsSymmetric(A,0.0,&flg);
44: MatViewFromOptions(A,NULL,"-A_view");
46: /* Create dense matrix C and X; C holds true solution with identical colums */
47: nrhs = 2;
48: PetscOptionsGetInt(NULL,NULL,"-nrhs",&nrhs,NULL);
49: PetscPrintf(PETSC_COMM_WORLD,"ex125: nrhs %D\n",nrhs);
50: MatCreate(PETSC_COMM_WORLD,&C);
51: MatSetOptionsPrefix(C,"rhs_");
52: MatSetSizes(C,m,PETSC_DECIDE,PETSC_DECIDE,nrhs);
53: MatSetType(C,MATDENSE);
54: MatSetFromOptions(C);
55: MatSetUp(C);
57: PetscOptionsGetBool(NULL,NULL,"-view_factor",&view,NULL);
58: PetscOptionsGetBool(NULL,NULL,"-test_matmatsolve",&testMatMatSolve,NULL);
59: PetscOptionsGetBool(NULL,NULL,"-cholesky",&chol,NULL);
60: #if defined(PETSC_HAVE_MUMPS)
61: PetscOptionsGetBool(NULL,NULL,"-test_mumps_opts",&test_mumps_opts,NULL);
62: #endif
64: PetscRandomCreate(PETSC_COMM_WORLD,&rand);
65: PetscRandomSetFromOptions(rand);
66: MatSetRandom(C,rand);
67: MatDuplicate(C,MAT_DO_NOT_COPY_VALUES,&X);
69: /* Create vectors */
70: VecCreate(PETSC_COMM_WORLD,&x);
71: VecSetSizes(x,n,PETSC_DECIDE);
72: VecSetFromOptions(x);
73: VecDuplicate(x,&b);
74: VecDuplicate(x,&u); /* save the true solution */
76: /* Test Factorization */
77: MatGetOrdering(A,MATORDERINGND,&perm,&iperm);
79: PetscOptionsGetInt(NULL,NULL,"-mat_solver_type",&ipack,NULL);
80: switch (ipack) {
81: #if defined(PETSC_HAVE_SUPERLU)
82: case 0:
83: if (chol) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"SuperLU does not provide Cholesky!");
84: PetscPrintf(PETSC_COMM_WORLD," SUPERLU LU:\n");
85: MatGetFactor(A,MATSOLVERSUPERLU,MAT_FACTOR_LU,&F);
86: matsolvexx = PETSC_TRUE;
87: break;
88: #endif
89: #if defined(PETSC_HAVE_SUPERLU_DIST)
90: case 1:
91: if (chol) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"SuperLU does not provide Cholesky!");
92: PetscPrintf(PETSC_COMM_WORLD," SUPERLU_DIST LU:\n");
93: MatGetFactor(A,MATSOLVERSUPERLU_DIST,MAT_FACTOR_LU,&F);
94: matsolvexx = PETSC_TRUE;
95: break;
96: #endif
97: #if defined(PETSC_HAVE_MUMPS)
98: case 2:
99: if (chol) {
100: PetscPrintf(PETSC_COMM_WORLD," MUMPS CHOLESKY:\n");
101: MatGetFactor(A,MATSOLVERMUMPS,MAT_FACTOR_CHOLESKY,&F);
102: } else {
103: PetscPrintf(PETSC_COMM_WORLD," MUMPS LU:\n");
104: MatGetFactor(A,MATSOLVERMUMPS,MAT_FACTOR_LU,&F);
105: }
106: matsolvexx = PETSC_TRUE;
107: if (test_mumps_opts) {
108: /* test mumps options */
109: PetscInt icntl;
110: PetscReal cntl;
112: icntl = 2; /* sequential matrix ordering */
113: MatMumpsSetIcntl(F,7,icntl);
115: cntl = 1.e-6; /* threshold for row pivot detection */
116: MatMumpsSetIcntl(F,24,1);
117: MatMumpsSetCntl(F,3,cntl);
118: }
119: break;
120: #endif
121: #if defined(PETSC_HAVE_MKL_PARDISO)
122: case 3:
123: if (chol) {
124: PetscPrintf(PETSC_COMM_WORLD," MKL_PARDISO CHOLESKY:\n");
125: MatGetFactor(A,MATSOLVERMKL_PARDISO,MAT_FACTOR_CHOLESKY,&F);
126: } else {
127: PetscPrintf(PETSC_COMM_WORLD," MKL_PARDISO LU:\n");
128: MatGetFactor(A,MATSOLVERMKL_PARDISO,MAT_FACTOR_LU,&F);
129: }
130: break;
131: #endif
132: default:
133: if (chol) {
134: PetscPrintf(PETSC_COMM_WORLD," PETSC CHOLESKY:\n");
135: MatGetFactor(A,MATSOLVERPETSC,MAT_FACTOR_CHOLESKY,&F);
136: } else {
137: PetscPrintf(PETSC_COMM_WORLD," PETSC LU:\n");
138: MatGetFactor(A,MATSOLVERPETSC,MAT_FACTOR_LU,&F);
139: }
140: matsolvexx = PETSC_TRUE;
141: }
143: MatFactorInfoInitialize(&info);
144: info.fill = 5.0;
145: info.shifttype = (PetscReal) MAT_SHIFT_NONE;
146: if (chol) {
147: MatCholeskyFactorSymbolic(F,A,perm,&info);
148: } else {
149: MatLUFactorSymbolic(F,A,perm,iperm,&info);
150: }
152: for (nfact = 0; nfact < 2; nfact++) {
153: if (chol) {
154: PetscPrintf(PETSC_COMM_WORLD," %D-the CHOLESKY numfactorization \n",nfact);
155: MatCholeskyFactorNumeric(F,A,&info);
156: } else {
157: PetscPrintf(PETSC_COMM_WORLD," %D-the LU numfactorization \n",nfact);
158: MatLUFactorNumeric(F,A,&info);
159: }
160: if (view) {
161: PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO);
162: MatView(F,PETSC_VIEWER_STDOUT_WORLD);
163: PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);
164: view = PETSC_FALSE;
165: }
167: #if defined(PETSC_HAVE_SUPERLU_DIST)
168: if (ipack == 1) { /* Test MatSuperluDistGetDiagU()
169: -- input: matrix factor F; output: main diagonal of matrix U on all processes */
170: PetscInt M;
171: PetscScalar *diag;
172: #if !defined(PETSC_USE_COMPLEX)
173: PetscInt nneg,nzero,npos;
174: #endif
176: MatGetSize(F,&M,NULL);
177: PetscMalloc1(M,&diag);
178: MatSuperluDistGetDiagU(F,diag);
179: PetscFree(diag);
181: #if !defined(PETSC_USE_COMPLEX)
182: /* Test MatGetInertia() */
183: MatGetInertia(F,&nneg,&nzero,&npos);
184: PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_WORLD," MatInertia: nneg: %D, nzero: %D, npos: %D\n",nneg,nzero,npos);
185: #endif
186: }
187: #endif
189: /* Test MatMatSolve() */
190: if (testMatMatSolve) {
191: if (!nfact) {
192: MatMatMult(A,C,MAT_INITIAL_MATRIX,2.0,&RHS);
193: } else {
194: MatMatMult(A,C,MAT_REUSE_MATRIX,2.0,&RHS);
195: }
196: for (nsolve = 0; nsolve < 2; nsolve++) {
197: PetscPrintf(PETSC_COMM_WORLD," %D-the MatMatSolve \n",nsolve);
198: MatMatSolve(F,RHS,X);
200: /* Check the error */
201: MatAXPY(X,-1.0,C,SAME_NONZERO_PATTERN);
202: MatNorm(X,NORM_FROBENIUS,&norm);
203: if (norm > tol) {
204: PetscPrintf(PETSC_COMM_WORLD,"%D-the MatMatSolve: Norm of error %g, nsolve %D\n",nsolve,(double)norm,nsolve);
205: }
206: }
207: if (matsolvexx) {
208: /* Test MatMatSolve(F,RHS,RHS), RHS is a dense matrix */
209: MatCopy(RHS,X,SAME_NONZERO_PATTERN);
210: MatMatSolve(F,X,X);
211: /* Check the error */
212: MatAXPY(X,-1.0,C,SAME_NONZERO_PATTERN);
213: MatNorm(X,NORM_FROBENIUS,&norm);
214: if (norm > tol) {
215: PetscPrintf(PETSC_COMM_WORLD,"MatMatSolve(F,RHS,RHS): Norm of error %g\n",(double)norm);
216: }
217: }
219: if (ipack == 2 && size == 1) {
220: Mat spRHS,spRHST,RHST;
222: MatTranspose(RHS,MAT_INITIAL_MATRIX,&RHST);
223: MatConvert(RHST,MATAIJ,MAT_INITIAL_MATRIX,&spRHST);
224: MatCreateTranspose(spRHST,&spRHS);
225: for (nsolve = 0; nsolve < 2; nsolve++) {
226: PetscPrintf(PETSC_COMM_WORLD," %D-the sparse MatMatSolve \n",nsolve);
227: MatMatSolve(F,spRHS,X);
229: /* Check the error */
230: MatAXPY(X,-1.0,C,SAME_NONZERO_PATTERN);
231: MatNorm(X,NORM_FROBENIUS,&norm);
232: if (norm > tol) {
233: PetscPrintf(PETSC_COMM_WORLD,"%D-the sparse MatMatSolve: Norm of error %g, nsolve %D\n",nsolve,(double)norm,nsolve);
234: }
235: }
236: MatDestroy(&spRHST);
237: MatDestroy(&spRHS);
238: MatDestroy(&RHST);
239: }
240: }
242: /* Test MatSolve() */
243: if (testMatSolve) {
244: for (nsolve = 0; nsolve < 2; nsolve++) {
245: VecSetRandom(x,rand);
246: VecCopy(x,u);
247: MatMult(A,x,b);
249: PetscPrintf(PETSC_COMM_WORLD," %D-the MatSolve \n",nsolve);
250: MatSolve(F,b,x);
252: /* Check the error */
253: VecAXPY(u,-1.0,x); /* u <- (-1.0)x + u */
254: VecNorm(u,NORM_2,&norm);
255: if (norm > tol) {
256: PetscReal resi;
257: MatMult(A,x,u); /* u = A*x */
258: VecAXPY(u,-1.0,b); /* u <- (-1.0)b + u */
259: VecNorm(u,NORM_2,&resi);
260: PetscPrintf(PETSC_COMM_WORLD,"MatSolve: Norm of error %g, resi %g, numfact %D\n",(double)norm,(double)resi,nfact);
261: }
262: }
263: }
264: }
266: /* Free data structures */
267: MatDestroy(&A);
268: MatDestroy(&C);
269: MatDestroy(&F);
270: MatDestroy(&X);
271: if (testMatMatSolve) {
272: MatDestroy(&RHS);
273: }
275: PetscRandomDestroy(&rand);
276: ISDestroy(&perm);
277: ISDestroy(&iperm);
278: VecDestroy(&x);
279: VecDestroy(&b);
280: VecDestroy(&u);
281: PetscFinalize();
282: return ierr;
283: }
286: /*TEST
288: test:
289: requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
290: args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type 10
291: output_file: output/ex125.out
293: test:
294: suffix: mkl_pardiso
295: requires: mkl_pardiso datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
296: args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type 3
298: test:
299: suffix: mumps
300: requires: mumps datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
301: args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type 2
302: output_file: output/ex125_mumps_seq.out
304: test:
305: suffix: mumps_2
306: nsize: 3
307: requires: mumps datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
308: args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type 2
309: output_file: output/ex125_mumps_par.out
311: test:
312: suffix: superlu_dist
313: requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES) superlu_dist
314: args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type 1 -mat_superlu_dist_rowperm NOROWPERM
316: test:
317: suffix: superlu_dist_2
318: nsize: 3
319: requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES) superlu_dist
320: args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type 1 -mat_superlu_dist_rowperm NOROWPERM
321: output_file: output/ex125_superlu_dist.out
323: test:
324: suffix: superlu_dist_complex
325: nsize: 3
326: requires: datafilespath superlu_dist complex double !define(PETSC_USE_64BIT_INDICES)
327: args: -f ${DATAFILESPATH}/matrices/farzad_B_rhs -mat_solver_type 1
328: output_file: output/ex125_superlu_dist_complex.out
330: TEST*/