Actual source code: ex125.c
petsc-3.14.6 2021-03-30
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,sizeof(file),&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: MatCreateVecs(A,&x,&b);
71: VecDuplicate(x,&u); /* save the true solution */
73: /* Test Factorization */
74: MatGetOrdering(A,MATORDERINGND,&perm,&iperm);
76: PetscOptionsGetInt(NULL,NULL,"-mat_solver_type",&ipack,NULL);
77: switch (ipack) {
78: #if defined(PETSC_HAVE_SUPERLU)
79: case 0:
80: if (chol) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"SuperLU does not provide Cholesky!");
81: PetscPrintf(PETSC_COMM_WORLD," SUPERLU LU:\n");
82: MatGetFactor(A,MATSOLVERSUPERLU,MAT_FACTOR_LU,&F);
83: matsolvexx = PETSC_TRUE;
84: break;
85: #endif
86: #if defined(PETSC_HAVE_SUPERLU_DIST)
87: case 1:
88: if (chol) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"SuperLU does not provide Cholesky!");
89: PetscPrintf(PETSC_COMM_WORLD," SUPERLU_DIST LU:\n");
90: MatGetFactor(A,MATSOLVERSUPERLU_DIST,MAT_FACTOR_LU,&F);
91: matsolvexx = PETSC_TRUE;
92: break;
93: #endif
94: #if defined(PETSC_HAVE_MUMPS)
95: case 2:
96: if (chol) {
97: PetscPrintf(PETSC_COMM_WORLD," MUMPS CHOLESKY:\n");
98: MatGetFactor(A,MATSOLVERMUMPS,MAT_FACTOR_CHOLESKY,&F);
99: } else {
100: PetscPrintf(PETSC_COMM_WORLD," MUMPS LU:\n");
101: MatGetFactor(A,MATSOLVERMUMPS,MAT_FACTOR_LU,&F);
102: }
103: matsolvexx = PETSC_TRUE;
104: if (test_mumps_opts) {
105: /* test mumps options */
106: PetscInt icntl;
107: PetscReal cntl;
109: icntl = 2; /* sequential matrix ordering */
110: MatMumpsSetIcntl(F,7,icntl);
112: cntl = 1.e-6; /* threshold for row pivot detection */
113: MatMumpsSetIcntl(F,24,1);
114: MatMumpsSetCntl(F,3,cntl);
115: }
116: break;
117: #endif
118: #if defined(PETSC_HAVE_MKL_PARDISO)
119: case 3:
120: if (chol) {
121: PetscPrintf(PETSC_COMM_WORLD," MKL_PARDISO CHOLESKY:\n");
122: MatGetFactor(A,MATSOLVERMKL_PARDISO,MAT_FACTOR_CHOLESKY,&F);
123: } else {
124: PetscPrintf(PETSC_COMM_WORLD," MKL_PARDISO LU:\n");
125: MatGetFactor(A,MATSOLVERMKL_PARDISO,MAT_FACTOR_LU,&F);
126: }
127: break;
128: #endif
129: #if defined(PETSC_HAVE_CUDA)
130: case 4:
131: if (chol) {
132: PetscPrintf(PETSC_COMM_WORLD," CUSPARSE CHOLESKY:\n");
133: MatGetFactor(A,MATSOLVERCUSPARSE,MAT_FACTOR_CHOLESKY,&F);
134: } else {
135: PetscPrintf(PETSC_COMM_WORLD," CUSPARSE LU:\n");
136: MatGetFactor(A,MATSOLVERCUSPARSE,MAT_FACTOR_LU,&F);
137: }
138: break;
139: #endif
140: default:
141: if (chol) {
142: PetscPrintf(PETSC_COMM_WORLD," PETSC CHOLESKY:\n");
143: MatGetFactor(A,MATSOLVERPETSC,MAT_FACTOR_CHOLESKY,&F);
144: } else {
145: PetscPrintf(PETSC_COMM_WORLD," PETSC LU:\n");
146: MatGetFactor(A,MATSOLVERPETSC,MAT_FACTOR_LU,&F);
147: }
148: matsolvexx = PETSC_TRUE;
149: }
151: MatFactorInfoInitialize(&info);
152: info.fill = 5.0;
153: info.shifttype = (PetscReal) MAT_SHIFT_NONE;
154: if (chol) {
155: MatCholeskyFactorSymbolic(F,A,perm,&info);
156: } else {
157: MatLUFactorSymbolic(F,A,perm,iperm,&info);
158: }
160: for (nfact = 0; nfact < 2; nfact++) {
161: if (chol) {
162: PetscPrintf(PETSC_COMM_WORLD," %D-the CHOLESKY numfactorization \n",nfact);
163: MatCholeskyFactorNumeric(F,A,&info);
164: } else {
165: PetscPrintf(PETSC_COMM_WORLD," %D-the LU numfactorization \n",nfact);
166: MatLUFactorNumeric(F,A,&info);
167: }
168: if (view) {
169: PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO);
170: MatView(F,PETSC_VIEWER_STDOUT_WORLD);
171: PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);
172: view = PETSC_FALSE;
173: }
175: #if defined(PETSC_HAVE_SUPERLU_DIST)
176: if (ipack == 1) { /* Test MatSuperluDistGetDiagU()
177: -- input: matrix factor F; output: main diagonal of matrix U on all processes */
178: PetscInt M;
179: PetscScalar *diag;
180: #if !defined(PETSC_USE_COMPLEX)
181: PetscInt nneg,nzero,npos;
182: #endif
184: MatGetSize(F,&M,NULL);
185: PetscMalloc1(M,&diag);
186: MatSuperluDistGetDiagU(F,diag);
187: PetscFree(diag);
189: #if !defined(PETSC_USE_COMPLEX)
190: /* Test MatGetInertia() */
191: MatGetInertia(F,&nneg,&nzero,&npos);
192: PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_WORLD," MatInertia: nneg: %D, nzero: %D, npos: %D\n",nneg,nzero,npos);
193: #endif
194: }
195: #endif
197: /* Test MatMatSolve() */
198: if (testMatMatSolve) {
199: if (!nfact) {
200: MatMatMult(A,C,MAT_INITIAL_MATRIX,2.0,&RHS);
201: } else {
202: MatMatMult(A,C,MAT_REUSE_MATRIX,2.0,&RHS);
203: }
204: for (nsolve = 0; nsolve < 2; nsolve++) {
205: PetscPrintf(PETSC_COMM_WORLD," %D-the MatMatSolve \n",nsolve);
206: MatMatSolve(F,RHS,X);
208: /* Check the error */
209: MatAXPY(X,-1.0,C,SAME_NONZERO_PATTERN);
210: MatNorm(X,NORM_FROBENIUS,&norm);
211: if (norm > tol) {
212: PetscPrintf(PETSC_COMM_WORLD,"%D-the MatMatSolve: Norm of error %g, nsolve %D\n",nsolve,(double)norm,nsolve);
213: }
214: }
215: if (matsolvexx) {
216: /* Test MatMatSolve(F,RHS,RHS), RHS is a dense matrix */
217: MatCopy(RHS,X,SAME_NONZERO_PATTERN);
218: MatMatSolve(F,X,X);
219: /* Check the error */
220: MatAXPY(X,-1.0,C,SAME_NONZERO_PATTERN);
221: MatNorm(X,NORM_FROBENIUS,&norm);
222: if (norm > tol) {
223: PetscPrintf(PETSC_COMM_WORLD,"MatMatSolve(F,RHS,RHS): Norm of error %g\n",(double)norm);
224: }
225: }
227: if (ipack == 2 && size == 1) {
228: Mat spRHS,spRHST,RHST;
230: MatTranspose(RHS,MAT_INITIAL_MATRIX,&RHST);
231: MatConvert(RHST,MATAIJ,MAT_INITIAL_MATRIX,&spRHST);
232: MatCreateTranspose(spRHST,&spRHS);
233: for (nsolve = 0; nsolve < 2; nsolve++) {
234: PetscPrintf(PETSC_COMM_WORLD," %D-the sparse MatMatSolve \n",nsolve);
235: MatMatSolve(F,spRHS,X);
237: /* Check the error */
238: MatAXPY(X,-1.0,C,SAME_NONZERO_PATTERN);
239: MatNorm(X,NORM_FROBENIUS,&norm);
240: if (norm > tol) {
241: PetscPrintf(PETSC_COMM_WORLD,"%D-the sparse MatMatSolve: Norm of error %g, nsolve %D\n",nsolve,(double)norm,nsolve);
242: }
243: }
244: MatDestroy(&spRHST);
245: MatDestroy(&spRHS);
246: MatDestroy(&RHST);
247: }
248: }
250: /* Test MatSolve() */
251: if (testMatSolve) {
252: for (nsolve = 0; nsolve < 2; nsolve++) {
253: VecSetRandom(x,rand);
254: VecCopy(x,u);
255: MatMult(A,x,b);
257: PetscPrintf(PETSC_COMM_WORLD," %D-the MatSolve \n",nsolve);
258: MatSolve(F,b,x);
260: /* Check the error */
261: VecAXPY(u,-1.0,x); /* u <- (-1.0)x + u */
262: VecNorm(u,NORM_2,&norm);
263: if (norm > tol) {
264: PetscReal resi;
265: MatMult(A,x,u); /* u = A*x */
266: VecAXPY(u,-1.0,b); /* u <- (-1.0)b + u */
267: VecNorm(u,NORM_2,&resi);
268: PetscPrintf(PETSC_COMM_WORLD,"MatSolve: Norm of error %g, resi %g, numfact %D\n",(double)norm,(double)resi,nfact);
269: }
270: }
271: }
272: }
274: /* Free data structures */
275: MatDestroy(&A);
276: MatDestroy(&C);
277: MatDestroy(&F);
278: MatDestroy(&X);
279: if (testMatMatSolve) {
280: MatDestroy(&RHS);
281: }
283: PetscRandomDestroy(&rand);
284: ISDestroy(&perm);
285: ISDestroy(&iperm);
286: VecDestroy(&x);
287: VecDestroy(&b);
288: VecDestroy(&u);
289: PetscFinalize();
290: return ierr;
291: }
294: /*TEST
296: test:
297: requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
298: args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type 10
299: output_file: output/ex125.out
301: test:
302: suffix: mkl_pardiso
303: requires: mkl_pardiso datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
304: args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type 3
306: test:
307: suffix: mumps
308: requires: mumps datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
309: args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type 2
310: output_file: output/ex125_mumps_seq.out
312: test:
313: suffix: mumps_2
314: nsize: 3
315: requires: mumps datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
316: args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type 2
317: output_file: output/ex125_mumps_par.out
319: test:
320: suffix: superlu_dist
321: requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES) superlu_dist
322: args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type 1 -mat_superlu_dist_rowperm NOROWPERM
324: test:
325: suffix: superlu_dist_2
326: nsize: 3
327: requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES) superlu_dist
328: args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type 1 -mat_superlu_dist_rowperm NOROWPERM
329: output_file: output/ex125_superlu_dist.out
331: test:
332: suffix: superlu_dist_complex
333: nsize: 3
334: requires: datafilespath superlu_dist complex double !define(PETSC_USE_64BIT_INDICES)
335: args: -f ${DATAFILESPATH}/matrices/farzad_B_rhs -mat_solver_type 1
336: output_file: output/ex125_superlu_dist_complex.out
338: test:
339: suffix: cusparse
340: requires: cuda datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
341: args: -mat_type aijcusparse -f ${DATAFILESPATH}/matrices/small -mat_solver_type 4 -cholesky {{0 1}separate output}
343: TEST*/