Actual source code: ex125.c
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;
10: PetscMPIInt size;
11: PetscInt m,n,nfact,nsolve,nrhs,ipack=0;
12: PetscReal norm,tol=1.e-10;
13: IS perm,iperm;
14: MatFactorInfo info;
15: PetscRandom rand;
16: PetscBool flg,testMatSolve=PETSC_TRUE,testMatMatSolve=PETSC_TRUE;
17: PetscBool chol=PETSC_FALSE,view=PETSC_FALSE,matsolvexx = PETSC_FALSE;
18: #if defined(PETSC_HAVE_MUMPS)
19: PetscBool test_mumps_opts=PETSC_FALSE;
20: #endif
21: PetscViewer fd; /* viewer */
22: char file[PETSC_MAX_PATH_LEN]; /* input file name */
24: PetscInitialize(&argc,&args,(char*)0,help);
25: MPI_Comm_size(PETSC_COMM_WORLD, &size);
27: /* Determine file from which we read the matrix A */
28: PetscOptionsGetString(NULL,NULL,"-f",file,sizeof(file),&flg);
29: if (flg) { /* Load matrix A */
30: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);
31: MatCreate(PETSC_COMM_WORLD,&A);
32: MatSetFromOptions(A);
33: MatLoad(A,fd);
34: PetscViewerDestroy(&fd);
35: } else {
36: n = 13;
37: PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
38: MatCreate(PETSC_COMM_WORLD,&A);
39: MatSetType(A,MATAIJ);
40: MatSetFromOptions(A);
41: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n);
42: MatSetUp(A);
43: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
44: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
45: MatShift(A,1.0);
46: }
47: MatGetLocalSize(A,&m,&n);
50: /* if A is symmetric, set its flag -- required by MatGetInertia() */
51: MatIsSymmetric(A,0.0,&flg);
53: MatViewFromOptions(A,NULL,"-A_view");
55: /* Create dense matrix C and X; C holds true solution with identical columns */
56: nrhs = 2;
57: PetscOptionsGetInt(NULL,NULL,"-nrhs",&nrhs,NULL);
58: PetscPrintf(PETSC_COMM_WORLD,"ex125: nrhs %" PetscInt_FMT "\n",nrhs);
59: MatCreate(PETSC_COMM_WORLD,&C);
60: MatSetOptionsPrefix(C,"rhs_");
61: MatSetSizes(C,m,PETSC_DECIDE,PETSC_DECIDE,nrhs);
62: MatSetType(C,MATDENSE);
63: MatSetFromOptions(C);
64: MatSetUp(C);
66: PetscOptionsGetBool(NULL,NULL,"-view_factor",&view,NULL);
67: PetscOptionsGetBool(NULL,NULL,"-test_matmatsolve",&testMatMatSolve,NULL);
68: PetscOptionsGetBool(NULL,NULL,"-cholesky",&chol,NULL);
69: #if defined(PETSC_HAVE_MUMPS)
70: PetscOptionsGetBool(NULL,NULL,"-test_mumps_opts",&test_mumps_opts,NULL);
71: #endif
73: PetscRandomCreate(PETSC_COMM_WORLD,&rand);
74: PetscRandomSetFromOptions(rand);
75: MatSetRandom(C,rand);
76: MatDuplicate(C,MAT_DO_NOT_COPY_VALUES,&X);
78: /* Create vectors */
79: MatCreateVecs(A,&x,&b);
80: VecDuplicate(x,&u); /* save the true solution */
82: /* Test Factorization */
83: MatGetOrdering(A,MATORDERINGND,&perm,&iperm);
85: PetscOptionsGetInt(NULL,NULL,"-mat_solver_type",&ipack,NULL);
86: switch (ipack) {
87: #if defined(PETSC_HAVE_SUPERLU)
88: case 0:
90: PetscPrintf(PETSC_COMM_WORLD," SUPERLU LU:\n");
91: MatGetFactor(A,MATSOLVERSUPERLU,MAT_FACTOR_LU,&F);
92: matsolvexx = PETSC_TRUE;
93: break;
94: #endif
95: #if defined(PETSC_HAVE_SUPERLU_DIST)
96: case 1:
98: PetscPrintf(PETSC_COMM_WORLD," SUPERLU_DIST LU:\n");
99: MatGetFactor(A,MATSOLVERSUPERLU_DIST,MAT_FACTOR_LU,&F);
100: matsolvexx = PETSC_TRUE;
101: break;
102: #endif
103: #if defined(PETSC_HAVE_MUMPS)
104: case 2:
105: if (chol) {
106: PetscPrintf(PETSC_COMM_WORLD," MUMPS CHOLESKY:\n");
107: MatGetFactor(A,MATSOLVERMUMPS,MAT_FACTOR_CHOLESKY,&F);
108: } else {
109: PetscPrintf(PETSC_COMM_WORLD," MUMPS LU:\n");
110: MatGetFactor(A,MATSOLVERMUMPS,MAT_FACTOR_LU,&F);
111: }
112: matsolvexx = PETSC_TRUE;
113: if (test_mumps_opts) {
114: /* test mumps options */
115: PetscInt icntl;
116: PetscReal cntl;
118: icntl = 2; /* sequential matrix ordering */
119: MatMumpsSetIcntl(F,7,icntl);
121: cntl = 1.e-6; /* threshold for row pivot detection */
122: MatMumpsSetIcntl(F,24,1);
123: MatMumpsSetCntl(F,3,cntl);
124: }
125: break;
126: #endif
127: #if defined(PETSC_HAVE_MKL_PARDISO)
128: case 3:
129: if (chol) {
130: PetscPrintf(PETSC_COMM_WORLD," MKL_PARDISO CHOLESKY:\n");
131: MatGetFactor(A,MATSOLVERMKL_PARDISO,MAT_FACTOR_CHOLESKY,&F);
132: } else {
133: PetscPrintf(PETSC_COMM_WORLD," MKL_PARDISO LU:\n");
134: MatGetFactor(A,MATSOLVERMKL_PARDISO,MAT_FACTOR_LU,&F);
135: }
136: break;
137: #endif
138: #if defined(PETSC_HAVE_CUDA)
139: case 4:
140: if (chol) {
141: PetscPrintf(PETSC_COMM_WORLD," CUSPARSE CHOLESKY:\n");
142: MatGetFactor(A,MATSOLVERCUSPARSE,MAT_FACTOR_CHOLESKY,&F);
143: } else {
144: PetscPrintf(PETSC_COMM_WORLD," CUSPARSE LU:\n");
145: MatGetFactor(A,MATSOLVERCUSPARSE,MAT_FACTOR_LU,&F);
146: }
147: break;
148: #endif
149: default:
150: if (chol) {
151: PetscPrintf(PETSC_COMM_WORLD," PETSC CHOLESKY:\n");
152: MatGetFactor(A,MATSOLVERPETSC,MAT_FACTOR_CHOLESKY,&F);
153: } else {
154: PetscPrintf(PETSC_COMM_WORLD," PETSC LU:\n");
155: MatGetFactor(A,MATSOLVERPETSC,MAT_FACTOR_LU,&F);
156: }
157: matsolvexx = PETSC_TRUE;
158: }
160: MatFactorInfoInitialize(&info);
161: info.fill = 5.0;
162: info.shifttype = (PetscReal) MAT_SHIFT_NONE;
163: if (chol) {
164: MatCholeskyFactorSymbolic(F,A,perm,&info);
165: } else {
166: MatLUFactorSymbolic(F,A,perm,iperm,&info);
167: }
169: for (nfact = 0; nfact < 2; nfact++) {
170: if (chol) {
171: PetscPrintf(PETSC_COMM_WORLD," %" PetscInt_FMT "-the CHOLESKY numfactorization \n",nfact);
172: MatCholeskyFactorNumeric(F,A,&info);
173: } else {
174: PetscPrintf(PETSC_COMM_WORLD," %" PetscInt_FMT "-the LU numfactorization \n",nfact);
175: MatLUFactorNumeric(F,A,&info);
176: }
177: if (view) {
178: PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO);
179: MatView(F,PETSC_VIEWER_STDOUT_WORLD);
180: PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);
181: view = PETSC_FALSE;
182: }
184: #if defined(PETSC_HAVE_SUPERLU_DIST)
185: if (ipack == 1) { /* Test MatSuperluDistGetDiagU()
186: -- input: matrix factor F; output: main diagonal of matrix U on all processes */
187: PetscInt M;
188: PetscScalar *diag;
189: #if !defined(PETSC_USE_COMPLEX)
190: PetscInt nneg,nzero,npos;
191: #endif
193: MatGetSize(F,&M,NULL);
194: PetscMalloc1(M,&diag);
195: MatSuperluDistGetDiagU(F,diag);
196: PetscFree(diag);
198: #if !defined(PETSC_USE_COMPLEX)
199: /* Test MatGetInertia() */
200: MatGetInertia(F,&nneg,&nzero,&npos);
201: PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_WORLD," MatInertia: nneg: %" PetscInt_FMT ", nzero: %" PetscInt_FMT ", npos: %" PetscInt_FMT "\n",nneg,nzero,npos);
202: #endif
203: }
204: #endif
206: #if defined(PETSC_HAVE_MUMPS)
207: /* mumps interface allows repeated call of MatCholeskyFactorSymbolic(), while the succession calls do nothing */
208: if (ipack == 2) {
209: if (chol) {
210: MatCholeskyFactorSymbolic(F,A,perm,&info);
211: MatCholeskyFactorNumeric(F,A,&info);
212: } else {
213: MatLUFactorSymbolic(F,A,perm,iperm,&info);
214: MatLUFactorNumeric(F,A,&info);
215: }
216: }
217: #endif
219: /* Test MatMatSolve() */
220: if (testMatMatSolve) {
221: if (!nfact) {
222: MatMatMult(A,C,MAT_INITIAL_MATRIX,2.0,&RHS);
223: } else {
224: MatMatMult(A,C,MAT_REUSE_MATRIX,2.0,&RHS);
225: }
226: for (nsolve = 0; nsolve < 2; nsolve++) {
227: PetscPrintf(PETSC_COMM_WORLD," %" PetscInt_FMT "-the MatMatSolve \n",nsolve);
228: MatMatSolve(F,RHS,X);
230: /* Check the error */
231: MatAXPY(X,-1.0,C,SAME_NONZERO_PATTERN);
232: MatNorm(X,NORM_FROBENIUS,&norm);
233: if (norm > tol) {
234: PetscPrintf(PETSC_COMM_WORLD,"%" PetscInt_FMT "-the MatMatSolve: Norm of error %g, nsolve %" PetscInt_FMT "\n",nsolve,(double)norm,nsolve);
235: }
236: }
237: if (matsolvexx) {
238: /* Test MatMatSolve(F,RHS,RHS), RHS is a dense matrix */
239: MatCopy(RHS,X,SAME_NONZERO_PATTERN);
240: MatMatSolve(F,X,X);
241: /* Check the error */
242: MatAXPY(X,-1.0,C,SAME_NONZERO_PATTERN);
243: MatNorm(X,NORM_FROBENIUS,&norm);
244: if (norm > tol) {
245: PetscPrintf(PETSC_COMM_WORLD,"MatMatSolve(F,RHS,RHS): Norm of error %g\n",(double)norm);
246: }
247: }
249: if (ipack == 2 && size == 1) {
250: Mat spRHS,spRHST,RHST;
252: MatTranspose(RHS,MAT_INITIAL_MATRIX,&RHST);
253: MatConvert(RHST,MATAIJ,MAT_INITIAL_MATRIX,&spRHST);
254: MatCreateTranspose(spRHST,&spRHS);
255: for (nsolve = 0; nsolve < 2; nsolve++) {
256: PetscPrintf(PETSC_COMM_WORLD," %" PetscInt_FMT "-the sparse MatMatSolve \n",nsolve);
257: MatMatSolve(F,spRHS,X);
259: /* Check the error */
260: MatAXPY(X,-1.0,C,SAME_NONZERO_PATTERN);
261: MatNorm(X,NORM_FROBENIUS,&norm);
262: if (norm > tol) {
263: PetscPrintf(PETSC_COMM_WORLD,"%" PetscInt_FMT "-the sparse MatMatSolve: Norm of error %g, nsolve %" PetscInt_FMT "\n",nsolve,(double)norm,nsolve);
264: }
265: }
266: MatDestroy(&spRHST);
267: MatDestroy(&spRHS);
268: MatDestroy(&RHST);
269: }
270: }
272: /* Test MatSolve() */
273: if (testMatSolve) {
274: for (nsolve = 0; nsolve < 2; nsolve++) {
275: VecSetRandom(x,rand);
276: VecCopy(x,u);
277: MatMult(A,x,b);
279: PetscPrintf(PETSC_COMM_WORLD," %" PetscInt_FMT "-the MatSolve \n",nsolve);
280: MatSolve(F,b,x);
282: /* Check the error */
283: VecAXPY(u,-1.0,x); /* u <- (-1.0)x + u */
284: VecNorm(u,NORM_2,&norm);
285: if (norm > tol) {
286: PetscReal resi;
287: MatMult(A,x,u); /* u = A*x */
288: VecAXPY(u,-1.0,b); /* u <- (-1.0)b + u */
289: VecNorm(u,NORM_2,&resi);
290: PetscPrintf(PETSC_COMM_WORLD,"MatSolve: Norm of error %g, resi %g, numfact %" PetscInt_FMT "\n",(double)norm,(double)resi,nfact);
291: }
292: }
293: }
294: }
296: /* Free data structures */
297: MatDestroy(&A);
298: MatDestroy(&C);
299: MatDestroy(&F);
300: MatDestroy(&X);
301: if (testMatMatSolve) {
302: MatDestroy(&RHS);
303: }
305: PetscRandomDestroy(&rand);
306: ISDestroy(&perm);
307: ISDestroy(&iperm);
308: VecDestroy(&x);
309: VecDestroy(&b);
310: VecDestroy(&u);
311: PetscFinalize();
312: return 0;
313: }
315: /*TEST
317: test:
318: requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
319: args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type 10
320: output_file: output/ex125.out
322: test:
323: suffix: 2
324: args: -mat_solver_type 10
325: output_file: output/ex125.out
327: test:
328: suffix: mkl_pardiso
329: requires: mkl_pardiso datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
330: args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type 3
332: test:
333: suffix: mkl_pardiso_2
334: requires: mkl_pardiso
335: args: -mat_solver_type 3
336: output_file: output/ex125_mkl_pardiso.out
338: test:
339: suffix: mumps
340: requires: mumps datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
341: args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type 2
342: output_file: output/ex125_mumps_seq.out
344: test:
345: suffix: mumps_2
346: nsize: 3
347: requires: mumps datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
348: args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type 2
349: output_file: output/ex125_mumps_par.out
351: test:
352: suffix: mumps_3
353: requires: mumps
354: args: -mat_solver_type 2
355: output_file: output/ex125_mumps_seq.out
357: test:
358: suffix: mumps_4
359: nsize: 3
360: requires: mumps
361: args: -mat_solver_type 2
362: output_file: output/ex125_mumps_par.out
364: test:
365: suffix: mumps_5
366: nsize: 3
367: requires: mumps
368: args: -mat_solver_type 2 -cholesky
369: output_file: output/ex125_mumps_par_cholesky.out
371: test:
372: suffix: superlu_dist
373: nsize: {{1 3}}
374: requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES) superlu_dist
375: args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type 1 -mat_superlu_dist_rowperm NOROWPERM
377: test:
378: suffix: superlu_dist_2
379: nsize: {{1 3}}
380: requires: superlu_dist !complex
381: args: -n 36 -mat_solver_type 1 -mat_superlu_dist_rowperm NOROWPERM
382: output_file: output/ex125_superlu_dist.out
384: test:
385: suffix: superlu_dist_complex
386: nsize: 3
387: requires: datafilespath superlu_dist complex double !defined(PETSC_USE_64BIT_INDICES)
388: args: -f ${DATAFILESPATH}/matrices/farzad_B_rhs -mat_solver_type 1
389: output_file: output/ex125_superlu_dist_complex.out
391: test:
392: suffix: superlu_dist_complex_2
393: nsize: 3
394: requires: superlu_dist complex
395: args: -mat_solver_type 1
396: output_file: output/ex125_superlu_dist_complex.out
398: test:
399: suffix: cusparse
400: requires: cuda datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
401: args: -mat_type aijcusparse -f ${DATAFILESPATH}/matrices/small -mat_solver_type 4 -cholesky {{0 1}separate output}
403: test:
404: suffix: cusparse_2
405: requires: cuda
406: args: -mat_type aijcusparse -mat_solver_type 4 -cholesky {{0 1}separate output}
408: TEST*/