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