Actual source code: ex125.c
petsc-3.8.4 2018-03-24
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: /* Create dense matrix C and X; C holds true solution with identical colums */
40: nrhs = 2;
41: PetscOptionsGetInt(NULL,NULL,"-nrhs",&nrhs,NULL);
42: if (!rank) PetscPrintf(PETSC_COMM_SELF,"ex125: nrhs %D\n",nrhs);
43: MatCreate(PETSC_COMM_WORLD,&C);
44: MatSetSizes(C,m,PETSC_DECIDE,PETSC_DECIDE,nrhs);
45: MatSetType(C,MATDENSE);
46: MatSetFromOptions(C);
47: MatSetUp(C);
48:
49: PetscRandomCreate(PETSC_COMM_WORLD,&rand);
50: PetscRandomSetFromOptions(rand);
51: /* #define DEBUGEX */
52: #if defined(DEBUGEX)
53: {
54: PetscInt row,j,M,cols[nrhs];
55: PetscScalar vals[nrhs];
56: MatGetSize(A,&M,NULL);
57: if (!rank) {
58: for (j=0; j<nrhs; j++) cols[j] = j;
59: for (row = 0; row < M; row++){
60: for (j=0; j<nrhs; j++) vals[j] = row;
61: MatSetValues(C,1,&row,nrhs,cols,vals,INSERT_VALUES);
62: }
63: }
64: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
65: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
66: }
67: #else
68: MatSetRandom(C,rand);
69: #endif
70: MatDuplicate(C,MAT_DO_NOT_COPY_VALUES,&X);
72: /* Create vectors */
73: VecCreate(PETSC_COMM_WORLD,&x);
74: VecSetSizes(x,n,PETSC_DECIDE);
75: VecSetFromOptions(x);
76: VecDuplicate(x,&b);
77: VecDuplicate(x,&u); /* save the true solution */
79: /* Test LU Factorization */
80: MatGetOrdering(A,MATORDERINGND,&perm,&iperm);
81: /*ISView(perm,PETSC_VIEWER_STDOUT_WORLD);*/
82: /*ISView(perm,PETSC_VIEWER_STDOUT_SELF);*/
84: PetscOptionsGetInt(NULL,NULL,"-mat_solver_package",&ipack,NULL);
85: switch (ipack) {
86: #if defined(PETSC_HAVE_SUPERLU)
87: case 0:
88: if (!rank) printf(" SUPERLU LU:\n");
89: MatGetFactor(A,MATSOLVERSUPERLU,MAT_FACTOR_LU,&F);
90: break;
91: #endif
92: #if defined(PETSC_HAVE_SUPERLU_DIST)
93: case 1:
94: if (!rank) printf(" SUPERLU_DIST LU:\n");
95: MatGetFactor(A,MATSOLVERSUPERLU_DIST,MAT_FACTOR_LU,&F);
96: break;
97: #endif
98: #if defined(PETSC_HAVE_MUMPS)
99: case 2:
100: if (!rank) printf(" MUMPS LU:\n");
101: MatGetFactor(A,MATSOLVERMUMPS,MAT_FACTOR_LU,&F);
102: {
103: /* test mumps options */
104: PetscInt icntl;
105: PetscReal cntl;
106:
107: icntl = 2; /* sequential matrix ordering */
108: MatMumpsSetIcntl(F,7,icntl);
110: cntl = 1.e-6; /* threshhold for row pivot detection */
111: MatMumpsSetIcntl(F,24,1);
112: MatMumpsSetCntl(F,3,cntl);
113: }
114: break;
115: #endif
116: #if defined(PETSC_HAVE_MKL_PARDISO)
117: case 3:
118: if (!rank) printf(" MKL_PARDISO LU:\n");
119: MatGetFactor(A,MATSOLVERMKL_PARDISO,MAT_FACTOR_LU,&F);
120: break;
121: #endif
122: default:
123: if (!rank) printf(" PETSC LU:\n");
124: MatGetFactor(A,MATSOLVERPETSC,MAT_FACTOR_LU,&F);
125: }
127: MatFactorInfoInitialize(&info);
128: info.fill = 5.0;
129: info.shifttype = (PetscReal) MAT_SHIFT_NONE;
130: MatLUFactorSymbolic(F,A,perm,iperm,&info);
132: for (nfact = 0; nfact < 2; nfact++) {
133: if (!rank) PetscPrintf(PETSC_COMM_SELF," %D-the LU numfactorization \n",nfact);
134: MatLUFactorNumeric(F,A,&info);
136: #if defined(PETSC_HAVE_SUPERLU_DIST)
137: { /* Test MatSuperluDistGetDiagU()
138: -- input: matrix factor F; output: main diagonal of matrix U on all processes */
139: PetscInt M;
140: PetscScalar *diag;
142: MatGetSize(F,&M,NULL);
143: PetscMalloc1(M,&diag);
144: MatSuperluDistGetDiagU(F,diag);
145: PetscFree(diag);
146: }
147: #endif
149: /* Test MatMatSolve() */
150: if (testMatMatSolve) {
151: if (!nfact) {
152: MatMatMult(A,C,MAT_INITIAL_MATRIX,2.0,&RHS);
153: } else {
154: MatMatMult(A,C,MAT_REUSE_MATRIX,2.0,&RHS);
155: }
156: for (nsolve = 0; nsolve < 2; nsolve++) {
157: if (!rank) PetscPrintf(PETSC_COMM_SELF," %D-the MatMatSolve \n",nsolve);
158: MatMatSolve(F,RHS,X);
160: /* Check the error */
161: MatAXPY(X,-1.0,C,SAME_NONZERO_PATTERN);
162: MatNorm(X,NORM_FROBENIUS,&norm);
163: if (norm > tol) {
164: if (!rank) {
165: PetscPrintf(PETSC_COMM_SELF,"%D-the MatMatSolve: Norm of error %g, nsolve %D\n",nsolve,norm,nsolve);
166: }
167: }
168: }
169: if (ipack == 2 && size == 1) {
170: Mat spRHS,spRHST,RHST;
172: MatTranspose(RHS,MAT_INITIAL_MATRIX,&RHST);
173: MatConvert(RHST,MATAIJ,MAT_INITIAL_MATRIX,&spRHST);
174: MatCreateTranspose(spRHST,&spRHS);
175: for (nsolve = 0; nsolve < 2; nsolve++) {
176: if (!rank) PetscPrintf(PETSC_COMM_SELF," %D-the sparse MatMatSolve \n",nsolve);
177: MatMatSolve(F,spRHS,X);
179: /* Check the error */
180: MatAXPY(X,-1.0,C,SAME_NONZERO_PATTERN);
181: MatNorm(X,NORM_FROBENIUS,&norm);
182: if (norm > tol) {
183: if (!rank) {
184: PetscPrintf(PETSC_COMM_SELF,"%D-the sparse MatMatSolve: Norm of error %g, nsolve %D\n",nsolve,norm,nsolve);
185: }
186: }
187: }
188: MatDestroy(&spRHST);
189: MatDestroy(&spRHS);
190: MatDestroy(&RHST);
191: }
192: }
194: /* Test MatSolve() */
195: if (testMatSolve) {
196: for (nsolve = 0; nsolve < 2; nsolve++) {
197: VecGetArray(x,&array);
198: for (i=0; i<m; i++) {
199: PetscRandomGetValue(rand,&rval);
200: array[i] = rval;
201: }
202: VecRestoreArray(x,&array);
203: VecCopy(x,u);
204: MatMult(A,x,b);
206: if (!rank) PetscPrintf(PETSC_COMM_SELF," %D-the MatSolve \n",nsolve);
207: MatSolve(F,b,x);
209: /* Check the error */
210: VecAXPY(u,-1.0,x); /* u <- (-1.0)x + u */
211: VecNorm(u,NORM_2,&norm);
212: if (norm > tol) {
213: PetscReal resi;
214: MatMult(A,x,u); /* u = A*x */
215: VecAXPY(u,-1.0,b); /* u <- (-1.0)b + u */
216: VecNorm(u,NORM_2,&resi);
217: if (!rank) {
218: PetscPrintf(PETSC_COMM_SELF,"MatSolve: Norm of error %g, resi %g, LU numfact %D\n",norm,resi,nfact);
219: }
220: }
221: }
222: }
223: }
225: /* Free data structures */
226: MatDestroy(&A);
227: MatDestroy(&C);
228: MatDestroy(&F);
229: MatDestroy(&X);
230: if (testMatMatSolve) {
231: MatDestroy(&RHS);
232: }
234: PetscRandomDestroy(&rand);
235: ISDestroy(&perm);
236: ISDestroy(&iperm);
237: VecDestroy(&x);
238: VecDestroy(&b);
239: VecDestroy(&u);
240: PetscFinalize();
241: return ierr;
242: }