Actual source code: ex125.c
petsc-3.9.4 2018-09-11
1:
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: PetscViewer fd; /* viewer */
21: char file[PETSC_MAX_PATH_LEN]; /* input file name */
23: PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
24: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
25: MPI_Comm_size(PETSC_COMM_WORLD, &size);
27: /* Determine file from which we read the matrix A */
28: PetscOptionsGetString(NULL,NULL,"-f",file,PETSC_MAX_PATH_LEN,&flg);
29: if (!flg) SETERRQ(PETSC_COMM_WORLD,1,"Must indicate binary file with the -f option");
31: /* Load matrix A */
32: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);
33: MatCreate(PETSC_COMM_WORLD,&A);
34: MatLoad(A,fd);
35: PetscViewerDestroy(&fd);
36: MatGetLocalSize(A,&m,&n);
37: if (m != n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ, "This example is not intended for rectangular matrices (%D, %D)", m, n);
39: /* if A is symmetric, set its flag -- required by MatGetInertia() */
40: MatIsSymmetric(A,0.0,&flg);
42: /* Create dense matrix C and X; C holds true solution with identical colums */
43: nrhs = 2;
44: PetscOptionsGetInt(NULL,NULL,"-nrhs",&nrhs,NULL);
45: PetscPrintf(PETSC_COMM_WORLD,"ex125: nrhs %D\n",nrhs);
46: MatCreate(PETSC_COMM_WORLD,&C);
47: MatSetSizes(C,m,PETSC_DECIDE,PETSC_DECIDE,nrhs);
48: MatSetType(C,MATDENSE);
49: MatSetFromOptions(C);
50: MatSetUp(C);
52: PetscRandomCreate(PETSC_COMM_WORLD,&rand);
53: PetscRandomSetFromOptions(rand);
54: /* #define DEBUGEX */
55: #if defined(DEBUGEX)
56: {
57: PetscInt row,j,M,cols[nrhs];
58: PetscScalar vals[nrhs];
59: MatGetSize(A,&M,NULL);
60: if (!rank) {
61: for (j=0; j<nrhs; j++) cols[j] = j;
62: for (row = 0; row < M; row++){
63: for (j=0; j<nrhs; j++) vals[j] = row;
64: MatSetValues(C,1,&row,nrhs,cols,vals,INSERT_VALUES);
65: }
66: }
67: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
68: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
69: }
70: #else
71: MatSetRandom(C,rand);
72: #endif
73: MatDuplicate(C,MAT_DO_NOT_COPY_VALUES,&X);
75: /* Create vectors */
76: VecCreate(PETSC_COMM_WORLD,&x);
77: VecSetSizes(x,n,PETSC_DECIDE);
78: VecSetFromOptions(x);
79: VecDuplicate(x,&b);
80: VecDuplicate(x,&u); /* save the true solution */
82: /* Test LU Factorization */
83: MatGetOrdering(A,MATORDERINGND,&perm,&iperm);
84: /*ISView(perm,PETSC_VIEWER_STDOUT_WORLD);*/
85: /*ISView(perm,PETSC_VIEWER_STDOUT_SELF);*/
87: PetscOptionsGetInt(NULL,NULL,"-mat_solver_type",&ipack,NULL);
88: switch (ipack) {
89: #if defined(PETSC_HAVE_SUPERLU)
90: case 0:
91: PetscPrintf(PETSC_COMM_WORLD," SUPERLU LU:\n");
92: MatGetFactor(A,MATSOLVERSUPERLU,MAT_FACTOR_LU,&F);
93: break;
94: #endif
95: #if defined(PETSC_HAVE_SUPERLU_DIST)
96: case 1:
97: PetscPrintf(PETSC_COMM_WORLD," SUPERLU_DIST LU:\n");
98: MatGetFactor(A,MATSOLVERSUPERLU_DIST,MAT_FACTOR_LU,&F);
99: break;
100: #endif
101: #if defined(PETSC_HAVE_MUMPS)
102: case 2:
103: PetscPrintf(PETSC_COMM_WORLD," MUMPS LU:\n");
104: MatGetFactor(A,MATSOLVERMUMPS,MAT_FACTOR_LU,&F);
105: {
106: /* test mumps options */
107: PetscInt icntl;
108: PetscReal cntl;
110: icntl = 2; /* sequential matrix ordering */
111: MatMumpsSetIcntl(F,7,icntl);
113: cntl = 1.e-6; /* threshhold for row pivot detection */
114: MatMumpsSetIcntl(F,24,1);
115: MatMumpsSetCntl(F,3,cntl);
116: }
117: break;
118: #endif
119: #if defined(PETSC_HAVE_MKL_PARDISO)
120: case 3:
121: PetscPrintf(PETSC_COMM_WORLD," MKL_PARDISO LU:\n");
122: MatGetFactor(A,MATSOLVERMKL_PARDISO,MAT_FACTOR_LU,&F);
123: break;
124: #endif
125: default:
126: PetscPrintf(PETSC_COMM_WORLD," PETSC LU:\n");
127: MatGetFactor(A,MATSOLVERPETSC,MAT_FACTOR_LU,&F);
128: }
130: MatFactorInfoInitialize(&info);
131: info.fill = 5.0;
132: info.shifttype = (PetscReal) MAT_SHIFT_NONE;
133: MatLUFactorSymbolic(F,A,perm,iperm,&info);
135: for (nfact = 0; nfact < 2; nfact++) {
136: PetscPrintf(PETSC_COMM_WORLD," %D-the LU numfactorization \n",nfact);
137: MatLUFactorNumeric(F,A,&info);
139: #if defined(PETSC_HAVE_SUPERLU_DIST)
140: if (ipack == 1) { /* Test MatSuperluDistGetDiagU()
141: -- input: matrix factor F; output: main diagonal of matrix U on all processes */
142: PetscInt M;
143: PetscScalar *diag;
144: #if !defined(PETSC_USE_COMPLEX)
145: PetscInt nneg,nzero,npos;
146: #endif
148: MatGetSize(F,&M,NULL);
149: PetscMalloc1(M,&diag);
150: MatSuperluDistGetDiagU(F,diag);
151: PetscFree(diag);
153: #if !defined(PETSC_USE_COMPLEX)
154: /* Test MatGetInertia() */
155: MatGetInertia(F,&nneg,&nzero,&npos);
156: if (!rank) {
157: PetscPrintf(PETSC_COMM_SELF," MatInertia: nneg: %D, nzero: %D, npos: %D\n",nneg,nzero,npos);
158: }
159: #endif
160: }
161: #endif
163: /* Test MatMatSolve() */
164: if (testMatMatSolve) {
165: if (!nfact) {
166: MatMatMult(A,C,MAT_INITIAL_MATRIX,2.0,&RHS);
167: } else {
168: MatMatMult(A,C,MAT_REUSE_MATRIX,2.0,&RHS);
169: }
170: for (nsolve = 0; nsolve < 2; nsolve++) {
171: PetscPrintf(PETSC_COMM_WORLD," %D-the MatMatSolve \n",nsolve);
172: MatMatSolve(F,RHS,X);
174: /* Check the error */
175: MatAXPY(X,-1.0,C,SAME_NONZERO_PATTERN);
176: MatNorm(X,NORM_FROBENIUS,&norm);
177: if (norm > tol) {
178: PetscPrintf(PETSC_COMM_WORLD,"%D-the MatMatSolve: Norm of error %g, nsolve %D\n",nsolve,norm,nsolve);
179: }
180: }
181: if (ipack == 2 && size == 1) {
182: Mat spRHS,spRHST,RHST;
184: MatTranspose(RHS,MAT_INITIAL_MATRIX,&RHST);
185: MatConvert(RHST,MATAIJ,MAT_INITIAL_MATRIX,&spRHST);
186: MatCreateTranspose(spRHST,&spRHS);
187: for (nsolve = 0; nsolve < 2; nsolve++) {
188: PetscPrintf(PETSC_COMM_WORLD," %D-the sparse MatMatSolve \n",nsolve);
189: MatMatSolve(F,spRHS,X);
191: /* Check the error */
192: MatAXPY(X,-1.0,C,SAME_NONZERO_PATTERN);
193: MatNorm(X,NORM_FROBENIUS,&norm);
194: if (norm > tol) {
195: PetscPrintf(PETSC_COMM_WORLD,"%D-the sparse MatMatSolve: Norm of error %g, nsolve %D\n",nsolve,norm,nsolve);
196: }
197: }
198: MatDestroy(&spRHST);
199: MatDestroy(&spRHS);
200: MatDestroy(&RHST);
201: }
202: }
204: /* Test MatSolve() */
205: if (testMatSolve) {
206: for (nsolve = 0; nsolve < 2; nsolve++) {
207: VecGetArray(x,&array);
208: for (i=0; i<m; i++) {
209: PetscRandomGetValue(rand,&rval);
210: array[i] = rval;
211: }
212: VecRestoreArray(x,&array);
213: VecCopy(x,u);
214: MatMult(A,x,b);
216: PetscPrintf(PETSC_COMM_WORLD," %D-the MatSolve \n",nsolve);
217: MatSolve(F,b,x);
219: /* Check the error */
220: VecAXPY(u,-1.0,x); /* u <- (-1.0)x + u */
221: VecNorm(u,NORM_2,&norm);
222: if (norm > tol) {
223: PetscReal resi;
224: MatMult(A,x,u); /* u = A*x */
225: VecAXPY(u,-1.0,b); /* u <- (-1.0)b + u */
226: VecNorm(u,NORM_2,&resi);
227: PetscPrintf(PETSC_COMM_WORLD,"MatSolve: Norm of error %g, resi %g, LU numfact %D\n",norm,resi,nfact);
228: }
229: }
230: }
231: }
233: /* Free data structures */
234: MatDestroy(&A);
235: MatDestroy(&C);
236: MatDestroy(&F);
237: MatDestroy(&X);
238: if (testMatMatSolve) {
239: MatDestroy(&RHS);
240: }
242: PetscRandomDestroy(&rand);
243: ISDestroy(&perm);
244: ISDestroy(&iperm);
245: VecDestroy(&x);
246: VecDestroy(&b);
247: VecDestroy(&u);
248: PetscFinalize();
249: return ierr;
250: }
253: /*TEST
255: test:
256: requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
257: args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type 10
258: output_file: output/ex125.out
260: test:
261: suffix: mkl_pardiso
262: requires: mkl_pardiso datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
263: args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type 3
265: test:
266: suffix: mumps
267: requires: mumps datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
268: args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type 2
269: output_file: output/ex125_mumps_seq.out
271: test:
272: suffix: mumps_2
273: nsize: 3
274: requires: mumps datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
275: args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type 2
276: output_file: output/ex125_mumps_par.out
278: test:
279: suffix: superlu_dist
280: requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES) superlu_dist
281: args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type 1 -mat_superlu_dist_rowperm NATURAL
283: test:
284: suffix: superlu_dist_2
285: nsize: 3
286: requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES) superlu_dist
287: args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type 1 -mat_superlu_dist_rowperm NATURAL
288: output_file: output/ex125_superlu_dist.out
290: test:
291: suffix: superlu_dist_complex
292: nsize: 3
293: requires: datafilespath superlu_dist complex double !define(PETSC_USE_64BIT_INDICES)
294: args: -f ${DATAFILESPATH}/matrices/farzad_B_rhs -mat_solver_type 1
295: output_file: output/ex125_superlu_dist_complex.out
297: TEST*/